Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

emutil.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "all_imageio.h"
00037 #include "portable_fileio.h"
00038 #include "emcache.h"
00039 #include "emdata.h"
00040 #include "ctf.h"
00041 #include "emassert.h"
00042 #include "exception.h"
00043 #include "hdf_filecache.h"
00044 
00045 #include <boost/shared_ptr.hpp>
00046 using boost::shared_ptr;
00047 
00048 #ifdef  HDFIO_CACHE
00049         #include <boost/filesystem.hpp>
00050         #include <boost/version.hpp>
00051 #endif  //HDFIO_CACHE
00052 
00053 #ifdef WIN32
00054         #include <windows.h>
00055         #define MAXPATHLEN (MAX_PATH*4)
00056 #else
00057         #include <sys/param.h>
00058 #endif  //WIN32
00059 
00060 //#ifdef EMAN2_USING_CUDA_MALLOC
00061 //#include "cuda/cuda_util.h"
00062 //#endif
00063 
00064 using namespace EMAN;
00065 
00066 static const int ATTR_NAME_LEN = 128;
00067 
00068 EMUtil::ImageType EMUtil::get_image_ext_type(const string & file_ext)
00069 {
00070         ENTERFUNC;
00071         static bool initialized = false;
00072         static map < string, ImageType > imagetypes;
00073 
00074         if (!initialized) {
00075                 imagetypes["rec"] = IMAGE_MRC;
00076                 imagetypes["mrc"] = IMAGE_MRC;
00077                 imagetypes["MRC"] = IMAGE_MRC;
00078                 imagetypes["ali"] = IMAGE_MRC;
00079                 imagetypes["st"] = IMAGE_MRC;           // IMOD stack file
00080 
00081                 imagetypes["tnf"] = IMAGE_MRC;
00082                 imagetypes["TNF"] = IMAGE_MRC;
00083 
00084                 imagetypes["ccp4"] = IMAGE_MRC;
00085                 imagetypes["map"] = IMAGE_MRC;
00086 
00087                 imagetypes["dm3"] = IMAGE_DM3;
00088                 imagetypes["DM3"] = IMAGE_DM3;
00089 
00090                 imagetypes["dm4"] = IMAGE_DM4;
00091                 imagetypes["DM4"] = IMAGE_DM4;
00092 
00093                 imagetypes["spi"] = IMAGE_SPIDER;
00094                 imagetypes["SPI"] = IMAGE_SPIDER;
00095 
00096                 imagetypes["spider"] = IMAGE_SPIDER;
00097                 imagetypes["SPIDER"] = IMAGE_SPIDER;
00098 
00099                 imagetypes["spidersingle"] = IMAGE_SINGLE_SPIDER;
00100                 imagetypes["SPIDERSINGLE"] = IMAGE_SINGLE_SPIDER;
00101 
00102                 imagetypes["singlespider"] = IMAGE_SINGLE_SPIDER;
00103                 imagetypes["SINGLESPIDER"] = IMAGE_SINGLE_SPIDER;
00104 
00105                 imagetypes["img"] = IMAGE_IMAGIC;
00106                 imagetypes["IMG"] = IMAGE_IMAGIC;
00107 
00108                 imagetypes["hed"] = IMAGE_IMAGIC;
00109                 imagetypes["HED"] = IMAGE_IMAGIC;
00110 
00111                 imagetypes["imagic"] = IMAGE_IMAGIC;
00112                 imagetypes["IMAGIC"] = IMAGE_IMAGIC;
00113 
00114                 imagetypes["pgm"] = IMAGE_PGM;
00115                 imagetypes["PGM"] = IMAGE_PGM;
00116 
00117                 imagetypes["lst"] = IMAGE_LST;
00118                 imagetypes["LST"] = IMAGE_LST;
00119 
00120                 imagetypes["lsx"] = IMAGE_LSTFAST;      // but .lst or another extension would also be ok
00121                 imagetypes["LSX"] = IMAGE_LSTFAST;
00122 
00123                 imagetypes["pif"] = IMAGE_PIF;
00124                 imagetypes["PIF"] = IMAGE_PIF;
00125 
00126                 imagetypes["png"] = IMAGE_PNG;
00127                 imagetypes["PNG"] = IMAGE_PNG;
00128 
00129                 imagetypes["h5"] = IMAGE_HDF;
00130                 imagetypes["H5"] = IMAGE_HDF;
00131 
00132                 imagetypes["hd5"] = IMAGE_HDF;
00133                 imagetypes["HD5"] = IMAGE_HDF;
00134 
00135                 imagetypes["hdf"] = IMAGE_HDF;
00136                 imagetypes["HDF"] = IMAGE_HDF;
00137 
00138                 imagetypes["tif"] = IMAGE_TIFF;
00139                 imagetypes["TIF"] = IMAGE_TIFF;
00140 
00141                 imagetypes["tiff"] = IMAGE_TIFF;
00142                 imagetypes["TIFF"] = IMAGE_TIFF;
00143 
00144                 imagetypes["fts"] = IMAGE_FITS;
00145                 imagetypes["FTS"] = IMAGE_FITS;
00146 
00147                 imagetypes["vtk"] = IMAGE_VTK;
00148                 imagetypes["VTK"] = IMAGE_VTK;
00149 
00150                 imagetypes["hdr"] = IMAGE_SAL;
00151                 imagetypes["HDR"] = IMAGE_SAL;
00152 
00153                 imagetypes["sal"] = IMAGE_SAL;
00154                 imagetypes["SAL"] = IMAGE_SAL;
00155 
00156                 imagetypes["map"] = IMAGE_ICOS;
00157                 imagetypes["MAP"] = IMAGE_ICOS;
00158 
00159                 imagetypes["icos"] = IMAGE_ICOS;
00160                 imagetypes["ICOS"] = IMAGE_ICOS;
00161 
00162                 imagetypes["am"] = IMAGE_AMIRA;
00163                 imagetypes["AM"] = IMAGE_AMIRA;
00164 
00165                 imagetypes["amira"] = IMAGE_AMIRA;
00166                 imagetypes["AMIRA"] = IMAGE_AMIRA;
00167 
00168                 imagetypes["emim"] = IMAGE_EMIM;
00169                 imagetypes["EMIM"] = IMAGE_EMIM;
00170 
00171                 imagetypes["xplor"] = IMAGE_XPLOR;
00172                 imagetypes["XPLOR"] = IMAGE_XPLOR;
00173 
00174                 imagetypes["em"] = IMAGE_EM;
00175                 imagetypes["EM"] = IMAGE_EM;
00176 
00177                 imagetypes["dm2"] = IMAGE_GATAN2;
00178                 imagetypes["DM2"] = IMAGE_GATAN2;
00179 
00180                 imagetypes["v4l"] = IMAGE_V4L;
00181                 imagetypes["V4L"] = IMAGE_V4L;
00182 
00183                 imagetypes["jpg"] = IMAGE_JPEG;
00184                 imagetypes["JPG"] = IMAGE_JPEG;
00185                 imagetypes["jpeg"] = IMAGE_JPEG;
00186                 imagetypes["JPEG"] = IMAGE_JPEG;
00187 
00188                 imagetypes["df3"] = IMAGE_DF3;
00189                 imagetypes["DF3"] = IMAGE_DF3;
00190 
00191                 imagetypes["Omap"] = IMAGE_OMAP;
00192                 imagetypes["omap"] = IMAGE_OMAP;
00193                 imagetypes["OMAP"] = IMAGE_OMAP;
00194                 imagetypes["BRIX"] = IMAGE_OMAP;
00195                 imagetypes["brix"] = IMAGE_OMAP;
00196                 imagetypes["DSN6"] = IMAGE_OMAP;
00197 
00198                 imagetypes["situs"] = IMAGE_SITUS;
00199                 imagetypes["SITUS"] = IMAGE_SITUS;
00200 
00201                 imagetypes["ser"] = IMAGE_SER;
00202                 imagetypes["SER"] = IMAGE_SER;
00203 
00204                 initialized = true;
00205         }
00206 
00207         ImageType result = IMAGE_UNKNOWN;
00208 
00209         if (imagetypes.find(file_ext) != imagetypes.end()) {
00210                 result = imagetypes[file_ext];
00211         }
00212 
00213         EXITFUNC;
00214         return result;
00215 }
00216 
00217 
00218 
00219 bool EMUtil::is_valid_filename(const string & filename) {
00220         ImageType type = get_image_ext_type(Util::get_filename_ext(filename));
00221         return (type != IMAGE_UNKNOWN);
00222 }
00223 
00224 EMUtil::ImageType EMUtil::fast_get_image_type(const string & filename,
00225                                                                                           const void *first_block,
00226                                                                                           off_t file_size)
00227 {
00228         ENTERFUNC;
00229         Assert(filename != "");
00230         Assert(first_block != 0);
00231         Assert(file_size > 0);
00232 
00233 #ifdef ENABLE_V4L2
00234         if (filename.compare(0,5,"/dev/")==0) return IMAGE_V4L;
00235 #endif
00236 
00237         string ext = Util::get_filename_ext(filename);
00238         if (ext == "") {
00239                 return IMAGE_UNKNOWN;
00240         }
00241         ImageType image_type = get_image_ext_type(ext);
00242 
00243         switch (image_type) {
00244         case IMAGE_MRC:
00245                 if (MrcIO::is_valid(first_block, file_size)) {
00246                         return IMAGE_MRC;
00247                 }
00248                 break;
00249         case IMAGE_DM3:
00250                 if (DM3IO::is_valid(first_block)) {
00251                         return IMAGE_DM3;
00252                 }
00253                 break;
00254         case IMAGE_DM4:
00255                 if (DM4IO::is_valid(first_block)) {
00256                         return IMAGE_DM4;
00257                 }
00258                 break;
00259 #ifdef EM_HDF5
00260         case IMAGE_HDF:
00261                 if (HdfIO2::is_valid(first_block)) {
00262                         return IMAGE_HDF;
00263                 }
00264                 break;
00265 #endif
00266         case IMAGE_LST:
00267                 if (LstIO::is_valid(first_block)) {
00268                         return IMAGE_LST;
00269                 }
00270                 break;
00271         case IMAGE_LSTFAST:
00272                 if (LstFastIO::is_valid(first_block)) {
00273                         return IMAGE_LSTFAST;
00274                 }
00275                 break;
00276 #ifdef EM_TIFF
00277         case IMAGE_TIFF:
00278                 if (TiffIO::is_valid(first_block)) {
00279                         return IMAGE_TIFF;
00280                 }
00281                 break;
00282 #endif
00283         case IMAGE_SPIDER:
00284                 if (SpiderIO::is_valid(first_block)) {
00285                         return IMAGE_SPIDER;
00286                 }
00287                 break;
00288         case IMAGE_SINGLE_SPIDER:
00289                 if (SingleSpiderIO::is_valid(first_block)) {
00290                         return IMAGE_SINGLE_SPIDER;
00291                 }
00292                 break;
00293         case IMAGE_PIF:
00294                 if (PifIO::is_valid(first_block)) {
00295                         return IMAGE_PIF;
00296                 }
00297                 break;
00298 #ifdef EM_PNG
00299         case IMAGE_PNG:
00300                 if (PngIO::is_valid(first_block)) {
00301                         return IMAGE_PNG;
00302                 }
00303                 break;
00304 #endif
00305         case IMAGE_VTK:
00306                 if (VtkIO::is_valid(first_block)) {
00307                         return IMAGE_VTK;
00308                 }
00309                 break;
00310         case IMAGE_PGM:
00311                 if (PgmIO::is_valid(first_block)) {
00312                         return IMAGE_PGM;
00313                 }
00314                 break;
00315         case IMAGE_ICOS:
00316                 if (IcosIO::is_valid(first_block)) {
00317                         return IMAGE_ICOS;
00318                 }
00319                 break;
00320         case IMAGE_SAL:
00321                 if (SalIO::is_valid(first_block)) {
00322                         return IMAGE_SAL;
00323                 }
00324                 break;
00325         case IMAGE_AMIRA:
00326                 if (AmiraIO::is_valid(first_block)) {
00327                         return IMAGE_AMIRA;
00328                 }
00329                 break;
00330         case IMAGE_XPLOR:
00331                 if (XplorIO::is_valid(first_block)) {
00332                         return IMAGE_XPLOR;
00333                 }
00334                 break;
00335         case IMAGE_GATAN2:
00336                 if (Gatan2IO::is_valid(first_block)) {
00337                         return IMAGE_GATAN2;
00338                 }
00339                 break;
00340         case IMAGE_EM:
00341                 if (EmIO::is_valid(first_block, file_size)) {
00342                         return IMAGE_EM;
00343                 }
00344                 break;
00345         case IMAGE_DF3:
00346                 if (EmIO::is_valid(first_block, file_size)) {
00347                         return IMAGE_DF3;
00348                 }
00349                 break;
00350         case IMAGE_OMAP:
00351                 if (OmapIO::is_valid(first_block, file_size)) {
00352                         return IMAGE_OMAP;
00353                 }
00354                 break;
00355         case IMAGE_SITUS:
00356                 if (SitusIO::is_valid(first_block)) {
00357                         return IMAGE_SITUS;
00358                 }
00359                 break;
00360         case IMAGE_SER:
00361                 if (SerIO::is_valid(first_block)) {
00362                         return IMAGE_SER;
00363                 }
00364                 break;
00365         case IMAGE_IMAGIC:
00366                 if (ImagicIO::is_valid(first_block)) {
00367                         return IMAGE_IMAGIC;
00368                 }
00369                 break;
00370         default:
00371                 return IMAGE_UNKNOWN;
00372         }
00373         EXITFUNC;
00374         return IMAGE_UNKNOWN;
00375 }
00376 
00377 
00378 EMUtil::ImageType EMUtil::get_image_type(const string & in_filename)
00379 {
00380         ENTERFUNC;
00381         Assert(in_filename != "");
00382 
00383 #ifdef ENABLE_V4L2
00384         if (in_filename.compare(0,5,"/dev/")==0) return IMAGE_V4L;
00385 #endif
00386 
00387         string filename = in_filename;
00388 
00389         string old_ext = Util::get_filename_ext(filename);
00390         if (old_ext == ImagicIO::IMG_EXT) {
00391                 filename = Util::change_filename_ext(filename, ImagicIO::HED_EXT);
00392         }
00393         else if(old_ext == "hdf") {
00394                 return IMAGE_HDF;
00395         }
00396 
00397         FILE *in = fopen(filename.c_str(), "rb");
00398         if (!in) {
00399                 throw FileAccessException(filename);
00400         }
00401 
00402         char first_block[1024];
00403         size_t n = fread(first_block, sizeof(char), sizeof(first_block), in);
00404         portable_fseek(in, 0, SEEK_END);
00405         off_t file_size = portable_ftell(in);
00406 
00407         if (n == 0) {
00408 //              This produces annoying console messages         
00409 //              LOGERR("file '%s' is an empty file", filename.c_str());
00410                 fclose(in);
00411                 return IMAGE_UNKNOWN;
00412         }
00413         fclose(in);
00414 
00415         ImageType image_type = fast_get_image_type(filename, first_block, file_size);
00416         if (image_type != IMAGE_UNKNOWN) {
00417                 return image_type;
00418         }
00419 
00420         if (SpiderIO::is_valid(first_block)) {
00421                 image_type = IMAGE_SPIDER;
00422         }
00423         else if (SingleSpiderIO::is_valid(first_block)) {
00424                 image_type = IMAGE_SINGLE_SPIDER;
00425         }
00426         else if (MrcIO::is_valid(first_block, file_size)) {
00427                 image_type = IMAGE_MRC;
00428         }
00429         else if (DM3IO::is_valid(first_block)) {
00430                 image_type = IMAGE_DM3;
00431         }
00432         else if (DM4IO::is_valid(first_block)) {
00433                 image_type = IMAGE_DM4;
00434         }
00435 #ifdef EM_HDF5
00436         else if (HdfIO2::is_valid(first_block)) {
00437                 image_type = IMAGE_HDF;
00438         }
00439 #endif
00440         else if (LstIO::is_valid(first_block)) {
00441                 image_type = IMAGE_LST;
00442         }
00443         else if (LstFastIO::is_valid(first_block)) {
00444                 image_type = IMAGE_LSTFAST;
00445         }
00446 #ifdef EM_TIFF
00447         else if (TiffIO::is_valid(first_block)) {
00448                 image_type = IMAGE_TIFF;
00449         }
00450 #endif
00451         else if (PifIO::is_valid(first_block)) {
00452                 image_type = IMAGE_PIF;
00453         }
00454 #ifdef EM_PNG
00455         else if (PngIO::is_valid(first_block)) {
00456                 image_type = IMAGE_PNG;
00457         }
00458 #endif
00459         else if (VtkIO::is_valid(first_block)) {
00460                 image_type = IMAGE_VTK;
00461         }
00462         else if (PgmIO::is_valid(first_block)) {
00463                 image_type = IMAGE_PGM;
00464         }
00465         else if (IcosIO::is_valid(first_block)) {
00466                 image_type = IMAGE_ICOS;
00467         }
00468         else if (SalIO::is_valid(first_block)) {
00469                 image_type = IMAGE_SAL;
00470         }
00471         else if (AmiraIO::is_valid(first_block)) {
00472                 image_type = IMAGE_AMIRA;
00473         }
00474         else if (XplorIO::is_valid(first_block)) {
00475                 image_type = IMAGE_XPLOR;
00476         }
00477         else if (Gatan2IO::is_valid(first_block)) {
00478                 image_type = IMAGE_GATAN2;
00479         }
00480         else if (FitsIO::is_valid(first_block)) {
00481                 image_type = IMAGE_FITS;
00482         }
00483         else if (EmIO::is_valid(first_block, file_size)) {
00484                 image_type = IMAGE_EM;
00485         }
00486         else if(OmapIO::is_valid(first_block, file_size)) {
00487                 image_type = IMAGE_OMAP;
00488         }
00489         else if(SitusIO::is_valid(first_block)) {
00490                 image_type = IMAGE_SITUS;
00491         }
00492         else if(SerIO::is_valid(first_block)) {
00493                 image_type = IMAGE_SER;
00494         }
00495         else if (ImagicIO::is_valid(first_block)) {
00496                 image_type = IMAGE_IMAGIC;
00497         }
00498         else if (Df3IO::is_valid(first_block)) {
00499                 image_type = IMAGE_DF3;
00500         }
00501         else {
00502                 //LOGERR("I don't know this image's type: '%s'", filename.c_str());
00503                 throw ImageFormatException("invalid image type");
00504         }
00505 
00506         EXITFUNC;
00507         return image_type;
00508 }
00509 
00510 
00511 int EMUtil::get_image_count(const string & filename)
00512 {
00513         ENTERFUNC;
00514         Assert(filename != "");
00515 
00516         int nimg = 0;
00517         ImageIO *imageio = get_imageio(filename, ImageIO::READ_ONLY);
00518 
00519         if (imageio) {
00520                 nimg = imageio->get_nimg();
00521         }
00522 
00523 #ifndef IMAGEIO_CACHE
00524         if( imageio )
00525         {
00526 #ifdef HDFIO_CACHE
00527                 if(dynamic_cast<HdfIO2*>(imageio)==NULL && dynamic_cast<HdfIO*>(imageio)==NULL) {
00528 #endif  //HDFIO_CACHE
00529                         delete imageio;
00530                         imageio = 0;
00531 #ifdef HDFIO_CACHE
00532                 }
00533 #endif  //HDFIO_CACHE
00534         }
00535 #endif  //IMAGEIO_CACHE
00536 
00537         EXITFUNC;
00538         return nimg;
00539 }
00540 
00541 
00542 ImageIO *EMUtil::get_imageio(const string & filename, int rw,
00543                                                          ImageType image_type)
00544 {
00545         ENTERFUNC;
00546         Assert(filename != "");
00547         Assert(rw == ImageIO::READ_ONLY ||
00548                    rw == ImageIO::READ_WRITE ||
00549                    rw == ImageIO::WRITE_ONLY);
00550 
00551         ImageIO *imageio = 0;
00552 #ifdef IMAGEIO_CACHE
00553         imageio = GlobalCache::instance()->get_imageio(filename, rw);
00554         if (imageio) {
00555                 return imageio;
00556         }
00557 #endif
00558 
00559         ImageIO::IOMode rw_mode = static_cast < ImageIO::IOMode > (rw);
00560 
00561         if (image_type == IMAGE_UNKNOWN) {
00562                 if(rw == ImageIO::WRITE_ONLY || rw == ImageIO::READ_WRITE) {
00563                         throw ImageFormatException("writing to this image format not supported.");
00564                 }
00565 
00566                 image_type = get_image_type(filename);
00567         }
00568 
00569         switch (image_type) {
00570 #ifdef ENABLE_V4L2
00571         case IMAGE_V4L:
00572                 imageio = new V4L2IO(filename, rw_mode);
00573                 break;
00574 #endif
00575         case IMAGE_MRC:
00576                 imageio = new MrcIO(filename, rw_mode);
00577                 break;
00578         case IMAGE_IMAGIC:
00579                 imageio = new ImagicIO2(filename, rw_mode);
00580                 if (rw_mode==ImageIO::READ_ONLY && ((ImagicIO2 *)imageio)->init_test()==-1 ) {
00581                         delete imageio;
00582                         imageio = new ImagicIO(filename, rw_mode);
00583                 }
00584                 break;
00585         case IMAGE_DM3:
00586                 imageio = new DM3IO(filename, rw_mode);
00587                 break;
00588         case IMAGE_DM4:
00589                 imageio = new DM4IO(filename, rw_mode);
00590                 break;
00591 #ifdef EM_TIFF
00592         case IMAGE_TIFF:
00593                 imageio = new TiffIO(filename, rw_mode);
00594                 break;
00595 #endif
00596 #ifdef EM_HDF5
00597         #ifdef  HDFIO_CACHE
00598         case IMAGE_HDF:
00599         {
00600                 bool readonly;
00601                 if(rw == ImageIO::READ_ONLY) {
00602                         readonly = true;
00603                 }
00604                 else{
00605                         readonly = false;
00606                 }
00607 
00608                 boost::filesystem::path p(filename);
00609 
00610 #if BOOST_VERSION >= 104600 && BOOST_FILESYSTEM_VERSION >= 3
00611                 boost::filesystem::path full_p = boost::filesystem::absolute(p);
00612 #else
00613                 boost::filesystem::path full_p = boost::filesystem::complete(p);
00614 #endif
00615 
00616                 FileItem * fitem = HDFCache::instance()->get_file(full_p.string());
00617                 if(fitem && readonly==fitem->get_readonly()) {
00618                         imageio = fitem->get_imgio();
00619                 }
00620                 else {
00621                         imageio = new HdfIO2(filename, rw_mode);
00622                         if (((HdfIO2 *)imageio)->init_test()==-1) {
00623                                 delete imageio;
00624                                 imageio = new HdfIO(filename, rw_mode);
00625                         }
00626 
00627                         if (imageio !=NULL) {
00628                                 FileItem * fitem = new FileItem(full_p.string(), imageio, time(0), readonly);
00629                                 HDFCache::instance()->add_file(fitem);
00630                         }
00631                 }
00632         }
00633                 break;
00634         #else   //HDFIO_CACHE
00635         case IMAGE_HDF:
00636                 imageio = new HdfIO2(filename, rw_mode);
00637                 if (((HdfIO2 *)imageio)->init_test()==-1) {
00638                         delete imageio;
00639                         imageio = new HdfIO(filename, rw_mode);
00640                 }
00641                 break;
00642         #endif  //HDFIO_CACHE
00643 #endif  //EM_HDF5
00644         case IMAGE_LST:
00645                 imageio = new LstIO(filename, rw_mode);
00646                 break;
00647         case IMAGE_LSTFAST:
00648                 imageio = new LstFastIO(filename, rw_mode);
00649                 break;
00650         case IMAGE_PIF:
00651                 imageio = new PifIO(filename, rw_mode);
00652                 break;
00653         case IMAGE_VTK:
00654                 imageio = new VtkIO(filename, rw_mode);
00655                 break;
00656         case IMAGE_SPIDER:
00657                 imageio = new SpiderIO(filename, rw_mode);
00658                 break;
00659         case IMAGE_SINGLE_SPIDER:
00660                 imageio = new SingleSpiderIO(filename, rw_mode);
00661                 break;
00662         case IMAGE_PGM:
00663                 imageio = new PgmIO(filename, rw_mode);
00664                 break;
00665 #ifdef EM_JPEG
00666         case IMAGE_JPEG:
00667                 imageio = new JpegIO(filename,rw_mode);
00668                 break;
00669 #endif
00670         case IMAGE_ICOS:
00671                 imageio = new IcosIO(filename, rw_mode);
00672                 break;
00673 #ifdef EM_PNG
00674         case IMAGE_PNG:
00675                 imageio = new PngIO(filename, rw_mode);
00676                 break;
00677 #endif
00678         case IMAGE_SAL:
00679                 imageio = new SalIO(filename, rw_mode);
00680                 break;
00681         case IMAGE_AMIRA:
00682                 imageio = new AmiraIO(filename, rw_mode);
00683                 break;
00684         case IMAGE_GATAN2:
00685                 imageio = new Gatan2IO(filename, rw_mode);
00686                 break;
00687         case IMAGE_EM:
00688                 imageio = new EmIO(filename, rw_mode);
00689                 break;
00690         case IMAGE_XPLOR:
00691                 imageio = new XplorIO(filename, rw_mode);
00692                 break;
00693         case IMAGE_FITS:
00694                 imageio = new FitsIO(filename, rw_mode);
00695                 break;
00696         case IMAGE_DF3:
00697                 imageio = new Df3IO(filename, rw_mode);
00698                 break;
00699         case IMAGE_OMAP:
00700                 imageio = new OmapIO(filename, rw_mode);
00701                 break;
00702         case IMAGE_SITUS:
00703                 imageio = new SitusIO(filename, rw_mode);
00704                 break;
00705         case IMAGE_SER:
00706                 imageio = new SerIO(filename, rw_mode);
00707                 break;
00708         default:
00709                 break;
00710         }
00711 #ifdef IMAGEIO_CACHE
00712         GlobalCache::instance()->add_imageio(filename, rw, imageio);
00713 #endif
00714         EXITFUNC;
00715         return imageio;
00716 }
00717 
00718 
00719 
00720 const char *EMUtil::get_imagetype_name(ImageType t)
00721 {
00722         switch (t) {
00723         case IMAGE_V4L:
00724                 return "V4L2";
00725                 break;
00726         case IMAGE_MRC:
00727                 return "MRC";
00728                 break;
00729         case IMAGE_SPIDER:
00730                 return "SPIDER";
00731                 break;
00732         case IMAGE_SINGLE_SPIDER:
00733                 return "Single-SPIDER";
00734                 break;
00735         case IMAGE_IMAGIC:
00736                 return "IMAGIC";
00737                 break;
00738         case IMAGE_PGM:
00739                 return "PGM";
00740                 break;
00741         case IMAGE_LST:
00742                 return "LST";
00743                 break;
00744         case IMAGE_LSTFAST:
00745                 return "Fast LST";
00746                 break;
00747         case IMAGE_PIF:
00748                 return "PIF";
00749                 break;
00750         case IMAGE_PNG:
00751                 return "PNG";
00752                 break;
00753         case IMAGE_HDF:
00754                 return "HDF5";
00755                 break;
00756         case IMAGE_DM3:
00757                 return "GatanDM3";
00758                 break;
00759         case IMAGE_DM4:
00760                 return "GatanDM4";
00761                 break;
00762         case IMAGE_TIFF:
00763                 return "TIFF";
00764                 break;
00765         case IMAGE_VTK:
00766                 return "VTK";
00767                 break;
00768         case IMAGE_SAL:
00769                 return "HDR";
00770                 break;
00771         case IMAGE_ICOS:
00772                 return "ICOS_MAP";
00773                 break;
00774         case IMAGE_EMIM:
00775                 return "EMIM";
00776                 break;
00777         case IMAGE_GATAN2:
00778                 return "GatanDM2";
00779                 break;
00780         case IMAGE_JPEG:
00781                 return "JPEG";
00782                 break;
00783         case IMAGE_AMIRA:
00784                 return "AmiraMesh";
00785                 break;
00786         case IMAGE_XPLOR:
00787                 return "XPLOR";
00788                 break;
00789         case IMAGE_EM:
00790                 return "EM";
00791                 break;
00792         case IMAGE_FITS:
00793                 return "FITS";
00794                 break;
00795         case IMAGE_DF3:
00796                 return "DF3";
00797                 break;
00798         case IMAGE_OMAP:
00799                 return "OMAP";
00800                 break;
00801         case IMAGE_SITUS:
00802                 return "SITUS";
00803                 break;
00804         case IMAGE_SER:
00805                 return "SER";
00806                 break;
00807         case IMAGE_UNKNOWN:
00808                 return "unknown";
00809         }
00810         return "unknown";
00811 }
00812 
00813 const char *EMUtil::get_datatype_string(EMDataType type)
00814 {
00815         switch (type) {
00816         case EM_CHAR:
00817                 return "CHAR";
00818         case EM_UCHAR:
00819                 return "UNSIGNED CHAR";
00820         case EM_SHORT:
00821                 return "SHORT";
00822         case EM_USHORT:
00823                 return "UNSIGNED SHORT";
00824         case EM_INT:
00825                 return "INT";
00826         case EM_UINT:
00827                 return "UNSIGNED INT";
00828         case EM_FLOAT:
00829                 return "FLOAT";
00830         case EM_DOUBLE:
00831                 return "DOUBLE";
00832         case EM_SHORT_COMPLEX:
00833                 return "SHORT_COMPLEX";
00834         case EM_USHORT_COMPLEX:
00835                 return "USHORT_COMPLEX";
00836         case EM_FLOAT_COMPLEX:
00837                 return "FLOAT_COMPLEX";
00838         case EM_UNKNOWN:
00839                 return "UNKNOWN";
00840         }
00841         return "UNKNOWN";
00842 }
00843 
00844 void EMUtil::get_region_dims(const Region * area, int nx, int *area_x,
00845                                                          int ny, int *area_y, int nz, int *area_z)
00846 {
00847         Assert(area_x);
00848         Assert(area_y);
00849 
00850         if (!area) {
00851                 *area_x = nx;
00852                 *area_y = ny;
00853                 if (area_z) {
00854                         *area_z = nz;
00855                 }
00856         }
00857         else {
00858                 Vec3i size = area->get_size();
00859                 *area_x = size[0];
00860                 *area_y = size[1];
00861 
00862                 if (area_z) {
00863                         if (area->get_ndim() > 2 && nz > 1) {
00864                                 *area_z = size[2];
00865                         }
00866                         else {
00867                                 *area_z = 1;
00868                         }
00869                 }
00870 
00871         }
00872 }
00873 
00874 void EMUtil::get_region_origins(const Region * area, int *p_x0, int *p_y0, int *p_z0,
00875                                                                 int nz, int image_index)
00876 {
00877         Assert(p_x0);
00878         Assert(p_y0);
00879 
00880         if (area) {
00881                 *p_x0 = static_cast < int >(area->origin[0]);
00882                 *p_y0 = static_cast < int >(area->origin[1]);
00883 
00884                 if (p_z0 && nz > 1 && area->get_ndim() > 2) {
00885                         *p_z0 = static_cast < int >(area->origin[2]);
00886                 }
00887         }
00888         else {
00889                 *p_x0 = 0;
00890                 *p_y0 = 0;
00891                 if (p_z0) {
00892                         *p_z0 = nz > 1 ? 0 : image_index;
00893                 }
00894         }
00895 }
00896 
00897 
00898 void EMUtil::process_region_io(void *vdata, FILE * file,
00899                                                            int rw_mode, int image_index,
00900                                                            size_t mode_size, int nx, int ny, int nz,
00901                                                            const Region * area, bool need_flip,
00902                                                            ImageType imgtype, int pre_row, int post_row)
00903 {
00904         Assert(vdata != 0);
00905         Assert(file != 0);
00906         Assert(rw_mode == ImageIO::READ_ONLY ||
00907                    rw_mode == ImageIO::READ_WRITE ||
00908                    rw_mode == ImageIO::WRITE_ONLY);
00909 
00910         if (mode_size == 0) throw UnexpectedBehaviorException("The mode size was 0?");
00911 
00912         unsigned char * cdata = (unsigned char *)vdata;
00913 
00914         int dx0 = 0; // data x0
00915         int dy0 = 0; // data y0
00916         int dz0 = 0; // data z0
00917 
00918         int fx0 = 0; // file x0
00919         int fy0 = 0; // file y0
00920         int fz0 = nz > 1 ? 0 : image_index; // file z0
00921 
00922 
00923         int xlen = 0;
00924         int ylen = 0;
00925         int zlen = 0;
00926         get_region_dims(area, nx, &xlen, ny, &ylen, nz, &zlen);
00927 
00928         if (area) { // Accommodate for all boundary overlaps of the region
00929 
00930                 Vec3i origin = area->get_origin();
00931 
00932 
00933                 fx0 = origin[0]; dx0 = origin[0];
00934                 fy0 = origin[1]; dy0 = origin[1];
00935                 if (nz > 1 && area->get_ndim() > 2) {
00936                         fz0 = origin[2]; dz0 = origin[2];
00937                 }
00938 
00939                 if (need_flip) {
00940                         Vec3i size = area->get_size();
00941                         fy0 = ny-(origin[1]+size[1]);
00942                 }
00943 
00944                 if (fx0 < 0) {
00945                         dx0 *= -1;
00946                         xlen = xlen + fx0; // because there are less reads
00947                         fx0 = 0;
00948                 }else {
00949                         dx0 = 0;
00950                         //fx0 *= -1;
00951                 }
00952                 if (fy0 < 0) {
00953                         dy0 *= -1;
00954                         ylen = ylen + fy0; // because there are less reads
00955                         fy0 = 0;
00956                 }else {
00957                         if (need_flip){
00958                                 dy0*=-1;
00959                         }
00960                         else dy0 = 0;
00961                         //fy0 *= -1;
00962                 }
00963                 if (fz0 < 0) {
00964                         dz0 *= -1;
00965                         zlen = zlen + fz0; // because there are less reads
00966                         fz0 = 0;
00967                 }else {
00968                         dz0 = 0;
00969                         //fz0 *= -1;
00970                 }
00971 
00972                 if ((fx0 + xlen)> nx) xlen = nx-fx0;
00973                 if ((fy0 + ylen)> ny) ylen = ny-fy0;
00974                 if ((fz0 + zlen)> nz) zlen = nz-fz0;
00975                 if ( xlen <= 0 || ylen <= 0 || zlen <= 0 ) return; // This is fine the region was entirely outside the image
00976         }
00977 
00978         if ( xlen <= 0 ) {
00979                 cout << "Xlen was too small " << xlen << endl;
00980                 return;
00981         }
00982 
00983         Vec3i size;
00984         if (area != 0) size = area->get_size();
00985         else size = Vec3d(nx,ny,nz);
00986 
00987         //size_t area_sec_size = xlen * ylen * mode_size;
00988         size_t memory_sec_size = size[0] * size[1] * mode_size;
00989         size_t img_row_size = nx * mode_size + pre_row + post_row;
00990         size_t area_row_size = xlen * mode_size;
00991         size_t memory_row_size = size[0] * mode_size;
00992 
00993         if ( area_row_size <= 0 ) {
00994                 cout << "Xlen was too small " << xlen << " mode_size " << mode_size << endl;
00995                 return;
00996         }
00997 
00998         size_t x_pre_gap = fx0 * mode_size;
00999         size_t x_post_gap = (nx - fx0 - xlen) * mode_size;
01000 
01001         size_t y_pre_gap = fy0 * img_row_size;
01002         size_t y_post_gap = (ny - fy0 - ylen) * img_row_size;
01003 
01004         portable_fseek(file, img_row_size * ny * fz0, SEEK_CUR);
01005 
01006         float nxlendata[1];
01007         int floatsize = (int) sizeof(float);
01008         nxlendata[0] = (float)(nx * floatsize);
01009 
01010         for (int k = dz0; k < (dz0+zlen); k++) {
01011                 if (y_pre_gap > 0) {
01012                         portable_fseek(file, y_pre_gap, SEEK_CUR);
01013                 }
01014                 //long k2 = k * area_sec_size;
01015                 long k2 = k*memory_sec_size;
01016 
01017                 for (int j = dy0; j < (dy0+ylen); j++) {
01018                         if (pre_row > 0) {
01019                                 if (imgtype == IMAGE_ICOS && rw_mode != ImageIO::READ_ONLY && !area) {
01020                                         fwrite(nxlendata, floatsize, 1, file);
01021                                 }
01022                                 else {
01023                                         portable_fseek(file, pre_row, SEEK_CUR);
01024                                 }
01025                         }
01026 
01027                         if (x_pre_gap > 0) {
01028                                 portable_fseek(file, x_pre_gap, SEEK_CUR);
01029                         }
01030 
01031                         int jj = j;
01032                         if (need_flip) {
01033                                 jj = (dy0+ylen) - 1 - j;
01034                                 if (dy0 > 0 ) { // region considerations add complications in the flipping scenario (imagic format)
01035                                         jj += dy0;
01036                                 }
01037                         }
01038 
01039                         if (rw_mode == ImageIO::READ_ONLY) {
01040                                 if (fread(&cdata[k2 + jj * memory_row_size+dx0*mode_size],
01041                                                   area_row_size, 1, file) != 1) {
01042                                         cout << jj << " " << k2 << " " << memory_row_size << " " << dx0 << " " << mode_size << " " << area_row_size << " " << cdata << "done" << endl;
01043                                         throw ImageReadException("Unknownfilename", "incomplete data read");
01044                                 }
01045                         }
01046                         else {
01047                                 if (fwrite(&cdata[k2 + jj * memory_row_size+dx0*mode_size],
01048                                                    area_row_size, 1, file) != 1) {
01049                                         throw ImageWriteException("", "incomplete data write");
01050                                 }
01051                         }
01052 
01053                         if (x_post_gap > 0) {
01054                                 portable_fseek(file, x_post_gap, SEEK_CUR);
01055                         }
01056 
01057                         if (post_row > 0) {
01058                                 if (imgtype == IMAGE_ICOS && rw_mode != ImageIO::READ_ONLY && !area) {
01059                                         fwrite(nxlendata, floatsize, 1, file);
01060                                 }
01061                                 else {
01062                                         portable_fseek(file, post_row, SEEK_CUR);
01063                                 }
01064                         }
01065                 }
01066 
01067                 if (y_post_gap > 0) {
01068                         portable_fseek(file, y_post_gap, SEEK_CUR);
01069                 }
01070         }
01071 }
01072 
01073 
01074 void EMUtil::dump_dict(const Dict & dict)
01075 {
01076         vector < string > keys = dict.keys();
01077         vector < EMObject > values = dict.values();
01078 
01079         for (unsigned int i = 0; i < keys.size(); i++) {
01080                 EMObject obj = values[i];
01081                 if( !obj.is_null() ) {
01082                         string val = obj.to_str();
01083 
01084                         if (keys[i] == "datatype") {
01085                                 val = get_datatype_string((EMDataType) (int) obj);
01086                         }
01087 
01088                         fprintf(stdout, "%25s\t%s\n", keys[i].c_str(), val.c_str());
01089                 }
01090         }
01091 }
01092 
01093 
01094 bool EMUtil::is_same_size(const EMData * const em1, const EMData * const em2)
01095 {
01096         if (em1->get_xsize() == em2->get_xsize() &&
01097                 em1->get_ysize() == em2->get_ysize() &&
01098                 em1->get_zsize() == em2->get_zsize()) {
01099                 return true;
01100         }
01101         return false;
01102 }
01103 
01104 bool EMUtil::is_complex_type(EMDataType datatype)
01105 {
01106         if (datatype == EM_SHORT_COMPLEX ||
01107                 datatype == EM_USHORT_COMPLEX ||
01108                 datatype == EM_FLOAT_COMPLEX) {
01109                 return true;
01110         }
01111         return false;
01112 }
01113 
01114 
01115 EMData *EMUtil::vertical_acf(const EMData * image, int maxdy)
01116 {
01117         if (!image) {
01118                 throw NullPointerException("NULL Image");
01119         }
01120 
01121         EMData *ret = new EMData();
01122         int nx = image->get_xsize();
01123         int ny = image->get_ysize();
01124 
01125         if (maxdy <= 1) {
01126                 maxdy = ny / 8;
01127         }
01128 
01129         ret->set_size(nx, maxdy, 1);
01130 
01131         float *data = image->get_data();
01132         float *ret_data = ret->get_data();
01133 
01134         for (int x = 0; x < nx; x++) {
01135                 for (int y = 0; y < maxdy; y++) {
01136                         float dot = 0;
01137                         for (int yy = maxdy; yy < ny - maxdy; yy++) {
01138                                 dot += data[x + (yy + y) * nx] * data[x + (yy - y) * nx];
01139                         }
01140                         ret_data[x + y * nx] = dot;
01141                 }
01142         }
01143 
01144         ret->update();
01145 
01146         return ret;
01147 }
01148 
01149 
01150 
01151 EMData *EMUtil::make_image_median(const vector < EMData * >&image_list)
01152 {
01153         if (image_list.size() == 0) {
01154                 return 0;
01155         }
01156 
01157         EMData *image0 = image_list[0];
01158         int image0_nx = image0->get_xsize();
01159         int image0_ny = image0->get_ysize();
01160         int image0_nz = image0->get_zsize();
01161         size_t size = (size_t)image0_nx * image0_ny * image0_nz;
01162 
01163         EMData *result = new EMData();
01164 
01165         result->set_size(image0_nx, image0_ny, image0_nz);
01166 
01167         float *dest = result->get_data();
01168         int nitems = static_cast < int >(image_list.size());
01169         float *srt = new float[nitems];
01170         float **src = new float *[nitems];
01171 
01172         for (int i = 0; i < nitems; i++) {
01173                 src[i] = image_list[i]->get_data();
01174         }
01175 
01176         for (size_t i = 0; i < size; ++i) {
01177                 for (int j = 0; j < nitems; j++) {
01178                         srt[j] = src[j][i];
01179                 }
01180 
01181                 for (int j = 0; j < nitems; j++) {
01182                         for (int k = j + 1; k < nitems; k++) {
01183                                 if (srt[j] < srt[k]) {
01184                                         float v = srt[j];
01185                                         srt[j] = srt[k];
01186                                         srt[k] = v;
01187                                 }
01188                         }
01189                 }
01190 
01191                 int l = nitems / 2;
01192                 if (nitems < 3) {
01193                         dest[i] = srt[l];
01194                 }
01195                 else {
01196                         dest[i] = (srt[l] + srt[l + 1] + srt[l - 1]) / 3.0f;
01197                 }
01198         }
01199 
01200         if( srt )
01201         {
01202                 delete[]srt;
01203                 srt = 0;
01204         }
01205         if( src )
01206         {
01207                 delete[]src;
01208                 src = 0;
01209         }
01210 
01211         result->update();
01212 
01213         return result;
01214 }
01215 
01216 bool EMUtil::is_same_ctf(const EMData * image1, const EMData * image2)
01217 {
01218         if (!image1) {
01219                 throw NullPointerException("image1 is NULL");
01220         }
01221         if (!image2) {
01222                 throw NullPointerException("image2 is NULL");
01223         }
01224 
01225         Ctf *ctf1 = image1->get_ctf();
01226         Ctf *ctf2 = image2->get_ctf();
01227 
01228         if ((!ctf1 && !ctf2) && (image1->has_ctff() == false && image2->has_ctff() == false)) {
01229                 return true;
01230         }
01231 
01232         if (ctf1 && ctf2) {
01233                 bool result = ctf1->equal(ctf2);
01234                 delete ctf1;
01235                 ctf1 = 0;
01236                 delete ctf2;
01237                 ctf2 = 0;
01238 
01239                 return result;
01240         }
01241         return false;
01242 }
01243 
01244 static int imgscore_cmp(const void *imgscore1, const void *imgscore2)
01245 {
01246         Assert(imgscore1 != 0);
01247         Assert(imgscore2 != 0);
01248 
01249         float c = ((ImageScore *)imgscore1)->score - ((ImageScore *)imgscore2)->score;
01250         if (c<0) {
01251                 return 1;
01252         }
01253         else if (c>0) {
01254                 return -1;
01255         }
01256         return 0;
01257 }
01258 
01259 ImageSort::ImageSort(int nn)
01260 {
01261         Assert(nn > 0);
01262         n = nn;
01263         image_scores = new ImageScore[n];
01264 }
01265 
01266 ImageSort::~ImageSort()
01267 {
01268         if( image_scores )
01269         {
01270                 delete [] image_scores;
01271                 image_scores = 0;
01272         }
01273 }
01274 
01275 void ImageSort::sort()
01276 {
01277         qsort(image_scores, n, sizeof(ImageScore), imgscore_cmp);
01278 
01279 }
01280 
01281 void ImageSort::set(int i, float score)
01282 {
01283         Assert(i >= 0);
01284         image_scores[i] = ImageScore(i, score);
01285 }
01286 
01287 int ImageSort::get_index(int i) const
01288 {
01289         Assert(i >= 0);
01290         return image_scores[i].index;
01291 }
01292 
01293 
01294 float ImageSort::get_score(int i) const
01295 {
01296         Assert(i >= 0);
01297         return image_scores[i].score;
01298 }
01299 
01300 
01301 int ImageSort::size() const
01302 {
01303         return n;
01304 }
01305 
01306 
01307 void EMUtil::process_ascii_region_io(float *data, FILE * file, int rw_mode,
01308                                                                          int , size_t mode_size, int nx, int ny, int nz,
01309                                                                          const Region * area, bool has_index_line,
01310                                                                          int nitems_per_line, const char *outformat)
01311 {
01312         Assert(data != 0);
01313         Assert(file != 0);
01314         Assert(rw_mode == ImageIO::READ_ONLY ||
01315                    rw_mode == ImageIO::READ_WRITE ||
01316                    rw_mode == ImageIO::WRITE_ONLY);
01317 
01318         int xlen = 0, ylen = 0, zlen = 0;
01319         get_region_dims(area, nx, &xlen, ny, &ylen, nz, &zlen);
01320 
01321         int x0 = 0;
01322         int y0 = 0;
01323         int z0 = 0;
01324 
01325         if (area) {
01326                 x0 = (int)area->origin[0];
01327                 y0 = (int)area->origin[1];
01328                 z0 = (int)area->origin[2];
01329         }
01330 
01331         int nlines_per_sec = (nx *ny) / nitems_per_line;
01332         int nitems_last_line = (nx * ny) % nitems_per_line;
01333         if (nitems_last_line != 0) {
01334                 nlines_per_sec++;
01335         }
01336 
01337         if (has_index_line) {
01338                 nlines_per_sec++;
01339         }
01340 
01341         if (z0 > 0) {
01342                 jump_lines(file, z0 * nlines_per_sec);
01343         }
01344 
01345 
01346         int nlines_pre_sec = (y0 * nx + x0) / nitems_per_line;
01347         int gap_nitems = nx - xlen;
01348         int ti = 0;
01349         int rlines = 0;
01350 
01351         for (int k = 0; k < zlen; k++) {
01352                 EMUtil::jump_lines(file, nlines_pre_sec+1);
01353 
01354                 int head_nitems = (y0 * nx + x0) % nitems_per_line;
01355                 int tail_nitems = 0;
01356                 bool is_head_read = false;
01357 
01358                 for (int j = 0; j < ylen; j++) {
01359 
01360                         if (head_nitems > 0 && !is_head_read) {
01361                                 EMUtil::process_numbers_io(file, rw_mode, nitems_per_line, mode_size,
01362                                                                                    nitems_per_line-head_nitems,
01363                                                                                    nitems_per_line-1, data, &ti, outformat);
01364                                 rlines++;
01365                         }
01366 
01367                         EMUtil::process_lines_io(file, rw_mode, nitems_per_line,
01368                                                                          mode_size, (xlen - head_nitems),
01369                                                                          data, &ti, outformat);
01370 
01371                         rlines += ((xlen - head_nitems)/nitems_per_line);
01372 
01373                         tail_nitems = (xlen - head_nitems) % nitems_per_line;
01374 
01375                         if ((gap_nitems + tail_nitems) > 0) {
01376                                 head_nitems = nitems_per_line -
01377                                         (gap_nitems + tail_nitems) % nitems_per_line;
01378                         }
01379                         else {
01380                                 head_nitems = 0;
01381                         }
01382 
01383                         is_head_read = false;
01384 
01385                         if (tail_nitems > 0) {
01386                                 if ((gap_nitems < (nitems_per_line-tail_nitems)) &&
01387                                         (j != (ylen-1))) {
01388                                         EMUtil::exclude_numbers_io(file, rw_mode, nitems_per_line,
01389                                                                                            mode_size, tail_nitems,
01390                                                                                            tail_nitems+gap_nitems-1, data, &ti, outformat);
01391                                         is_head_read = true;
01392                                         rlines++;
01393                                 }
01394                                 else {
01395                                         EMUtil::process_numbers_io(file, rw_mode, nitems_per_line, mode_size,
01396                                                                                            0, tail_nitems-1, data, &ti, outformat);
01397                                         rlines++;
01398                                 }
01399                         }
01400 
01401                         if (gap_nitems > (nitems_per_line-tail_nitems)) {
01402                                 int gap_nlines = (gap_nitems - (nitems_per_line-tail_nitems)) /
01403                                         nitems_per_line;
01404                                 if (gap_nlines > 0 && j != (ylen-1)) {
01405                                         EMUtil::jump_lines(file, gap_nlines);
01406                                 }
01407                         }
01408                 }
01409 
01410                 int ytail_nitems = (ny-ylen-y0) * nx + (nx-xlen-x0) - (nitems_per_line-tail_nitems);
01411                 EMUtil::jump_lines_by_items(file, ytail_nitems, nitems_per_line);
01412         }
01413 }
01414 
01415 
01416 void EMUtil::jump_lines_by_items(FILE * file, int nitems, int nitems_per_line)
01417 {
01418         Assert(file);
01419         Assert(nitems_per_line > 0);
01420 
01421         if (nitems <= 0) {
01422                 return;
01423         }
01424 
01425         int nlines = nitems / nitems_per_line;
01426         if ((nitems % nitems_per_line) != 0) {
01427                 nlines++;
01428         }
01429         if (nlines > 0) {
01430                 jump_lines(file, nlines);
01431         }
01432 }
01433 
01434 
01435 void EMUtil::jump_lines(FILE * file, int nlines)
01436 {
01437         Assert(file);
01438 
01439         if (nlines > 0) {
01440                 char line[MAXPATHLEN];
01441                 for (int l = 0; l < nlines; l++) {
01442                         if (!fgets(line, sizeof(line), file)) {
01443                                 Assert("read xplor file failed");
01444                         }
01445                 }
01446         }
01447 }
01448 
01449 void EMUtil::process_numbers_io(FILE * file, int rw_mode,
01450                                                                 int nitems_per_line, size_t mode_size, int start,
01451                                                                 int end, float *data, int *p_i, const char * outformat)
01452 {
01453         Assert(file);
01454         Assert(start >= 0);
01455         Assert(start <= end);
01456         Assert(end <= nitems_per_line);
01457         Assert(data);
01458         Assert(p_i);
01459         Assert(outformat);
01460 
01461         char line[MAXPATHLEN];
01462 
01463         if (rw_mode == ImageIO::READ_ONLY) {
01464                 if (!fgets(line, sizeof(line), file)) {
01465                         Assert("read xplor file failed");
01466                 }
01467 
01468                 int nitems_in_line = (int) (strlen(line) / mode_size);
01469                 Assert(end <= nitems_in_line);
01470                 vector<float> d(nitems_in_line);
01471                 char * pline = line;
01472 
01473                 for (int i = 0; i < nitems_in_line; i++) {
01474                         sscanf(pline, "%f", &d[i]);
01475                         pline += (int)mode_size;
01476                 }
01477 
01478 
01479                 for (int i = start; i <= end; i++) {
01480                         data[*p_i] = d[i];
01481                         (*p_i)++;
01482                 }
01483         }
01484         else {
01485                 portable_fseek(file, mode_size * start, SEEK_CUR);
01486                 for (int i = start; i <= end; i++) {
01487                         fprintf(file, outformat, data[*p_i]);
01488                         (*p_i)++;
01489                 }
01490 
01491                 portable_fseek(file, mode_size * (nitems_per_line - end-1)+1, SEEK_CUR);
01492         }
01493 }
01494 
01495 
01496 void EMUtil::exclude_numbers_io(FILE * file, int rw_mode,
01497                                                                 int nitems_per_line, size_t mode_size, int start,
01498                                                                 int end, float * data, int *p_i, const char * outformat)
01499 {
01500         Assert(file);
01501         Assert(mode_size > 0);
01502         Assert(start >= 0);
01503         Assert(end <= nitems_per_line);
01504         Assert(data);
01505         Assert(p_i);
01506         Assert(outformat);
01507 
01508         char line[MAXPATHLEN];
01509 
01510         if (rw_mode == ImageIO::READ_ONLY) {
01511 
01512                 if (!fgets(line, sizeof(line), file)) {
01513                         Assert("read xplor file failed");
01514                 }
01515 
01516                 int nitems_in_line =  (int) (strlen(line) / mode_size);
01517                 Assert(end <= nitems_in_line);
01518 
01519                 vector<float> d(nitems_in_line);
01520                 char *pline = line;
01521 
01522                 for (int i = 0; i < nitems_in_line; i++) {
01523                         sscanf(pline, "%f", &d[i]);
01524                         pline = pline + (int)mode_size;
01525                 }
01526 
01527 
01528                 for (int i = 0; i < start; i++) {
01529                         data[*p_i] = d[i];
01530                         (*p_i)++;
01531                 }
01532 
01533                 for (int i = end+1; i < nitems_in_line; i++) {
01534                         data[*p_i] = d[i];
01535                         (*p_i)++;
01536                 }
01537         }
01538         else {
01539                 for (int i = 0; i < start; i++) {
01540                         fprintf(file, outformat, data[*p_i]);
01541                         (*p_i)++;
01542                 }
01543 
01544                 portable_fseek(file, (end-start+1) * mode_size, SEEK_CUR);
01545 
01546                 for (int i = end+1; i < nitems_per_line; i++) {
01547                         fprintf(file, outformat, data[*p_i]);
01548                         (*p_i)++;
01549                 }
01550                 portable_fseek(file, 1, SEEK_CUR);
01551         }
01552 }
01553 
01554 void EMUtil::process_lines_io(FILE * file, int rw_mode,
01555                                                           int nitems_per_line, size_t mode_size,
01556                                                           int nitems, float *data, int *p_i,
01557                                                           const char * outformat)
01558 {
01559         Assert(file);
01560         Assert(data);
01561         Assert(p_i);
01562 
01563         if (nitems > 0) {
01564                 int nlines = nitems / nitems_per_line;
01565                 for (int i = 0; i < nlines; i++) {
01566                         EMUtil::process_numbers_io(file, rw_mode, nitems_per_line, mode_size, 0,
01567                                                                            nitems_per_line-1, data, p_i, outformat);
01568                 }
01569         }
01570 }
01571 
01572 vector<string> EMUtil::get_euler_names(const string & euler_type)
01573 {
01574     vector<string> v;
01575     string b = "euler_";
01576 
01577     if (euler_type == "EMAN") {
01578         v.push_back(b + "alt");
01579         v.push_back(b + "az");
01580         v.push_back(b + "phi");
01581     }
01582     else if (euler_type == "MRC") {
01583         v.push_back(b + "theta");
01584         v.push_back(b + "phi");
01585         v.push_back(b + "omega");
01586     }
01587     else if (euler_type == "IMAGIC") {
01588         v.push_back(b + "alpha");
01589         v.push_back(b + "beta");
01590         v.push_back(b + "gamma");
01591     }
01592     else if (euler_type == "SPIDER") {
01593         v.push_back(b + "phi");
01594         v.push_back(b + "theta");
01595         v.push_back(b + "gamma");
01596     }
01597     else if (euler_type == "SPIN" ||
01598              euler_type == "SGIROT") {
01599         v.push_back(b + "q");
01600         v.push_back(b + "n1");
01601         v.push_back(b + "n2");
01602         v.push_back(b + "n3");
01603     }
01604 
01605     else if (euler_type == "QUATERNION") {
01606         v.push_back(b + "e0");
01607         v.push_back(b + "e1");
01608         v.push_back(b + "e2");
01609         v.push_back(b + "e3");
01610     }
01611 
01612     return v;
01613 }
01614 
01615 
01616 vector<EMObject> EMUtil::get_all_attributes(const string & file_name, const string & attr_name)
01617 {
01618         vector<EMObject> v;
01619 
01620         Assert(file_name != "");
01621         Assert(attr_name != "");
01622 
01623         vector< shared_ptr<EMData> > vpImg = EMData::read_images(file_name, vector<int>(), true);
01624         vector< shared_ptr<EMData> >::const_iterator iter;
01625         for(iter = vpImg.begin(); iter!=vpImg.end(); ++iter) {
01626                 v.push_back((*iter)->get_attr_default(attr_name));
01627         }
01628 
01629         return v;
01630 }
01631 
01632 void EMUtil::getRenderMinMax(float * data, const int nx, const int ny, float& rendermin, float& rendermax, const int nz)
01633 {
01634 #ifdef _WIN32
01635         if (rendermax<=rendermin || _isnan(rendermin) || _isnan(rendermax)) {
01636 #else
01637         if (rendermax<=rendermin || std::isnan(rendermin) || std::isnan(rendermax)) {
01638 #endif
01639                 float m=0.0f,s=0.0f;
01640 
01641                 size_t size = (size_t)nx*ny*nz;
01642                 float min=data[0],max=data[0];
01643 
01644                 for (size_t i=0; i<size; ++i) { m+=data[i]; s+=data[i]*data[i]; min=data[i]<min?data[i]:min; max=data[i]>max?data[i]:max; }
01645                 m/=(float)(size);
01646                 s=sqrt(s/(float)(size)-m*m);
01647 #ifdef _WIN32
01648                 if (s<=0 || _isnan(s)) s=1.0;   // this means all data values are the same
01649 #else
01650                 if (s<=0 || std::isnan(s)) s=1.0;       // this means all data values are the same
01651 #endif  //_WIN32
01652                 rendermin=m-s*5.0f;
01653                 rendermax=m+s*5.0f;
01654                 if (rendermin<=min) rendermin=min;
01655                 if (rendermax>=max) rendermax=max;
01656         }
01657 }
01658 
01659 #ifdef EM_HDF5
01660 EMObject EMUtil::read_hdf_attribute(const string & filename, const string & key, int image_index)
01661 {
01662         ImageType image_type = get_image_type(filename);
01663         if(image_type != IMAGE_HDF) {
01664                 throw ImageFormatException("This function only applies to HDF5 file.");
01665         }
01666 
01667         HdfIO2* imageio = new HdfIO2(filename, ImageIO::READ_ONLY);
01668         imageio->init();
01669 
01670         // Each image is in a group for later expansion. Open the group
01671         hid_t file = imageio->get_fileid();
01672         char ipath[50];
01673         sprintf(ipath,"/MDF/images/%d",image_index);
01674         hid_t igrp=H5Gopen(file,ipath);
01675 
01676         if (igrp<0) {   //group not existed
01677                 throw _NotExistingObjectException(string(ipath));
01678         }
01679 
01680         string s("EMAN.");
01681         s += key;
01682         hid_t attr = H5Aopen_name(igrp, s.c_str());
01683         EMObject emobj = imageio->read_attr(attr);
01684 
01685         H5Aclose(attr);
01686         H5Gclose(igrp);
01687         delete imageio;
01688 
01689         return emobj;
01690 }
01691 
01692 int EMUtil::write_hdf_attribute(const string & filename, const string & key, EMObject value, int image_index)
01693 {
01694         ImageType image_type = get_image_type(filename);
01695         if(image_type != IMAGE_HDF) {
01696                 throw ImageFormatException("This function only applies to HDF5 file.");
01697         }
01698 
01699         HdfIO2* imageio = new HdfIO2(filename, ImageIO::WRITE_ONLY);
01700         imageio->init();
01701 
01702         // Each image is in a group for later expansion. Open the group
01703         hid_t file = imageio->get_fileid();
01704         char ipath[50];
01705         sprintf(ipath,"/MDF/images/%d",image_index);
01706         hid_t igrp=H5Gopen(file,ipath);
01707 
01708         if (igrp<0) {   //group not existed
01709                 throw _NotExistingObjectException(string(ipath));
01710         }
01711 
01712         string s("EMAN.");
01713         s += key;
01714         int ret = imageio->write_attr(igrp, s.c_str(), value);
01715 
01716         H5Gclose(igrp);
01717         delete imageio;
01718 
01719         return ret;
01720 }
01721 
01722 int EMUtil::delete_hdf_attribute(const string & filename, const string & key, int image_index)
01723 {
01724         ImageType image_type = get_image_type(filename);
01725         if(image_type != IMAGE_HDF) {
01726                 throw ImageFormatException("This function only applies to HDF5 file.");
01727         }
01728 
01729         HdfIO2* imageio = new HdfIO2(filename, ImageIO::READ_WRITE);
01730         imageio->init();
01731 
01732         // Each image is in a group for later expansion. Open the group
01733         hid_t file = imageio->get_fileid();
01734         char ipath[50];
01735         sprintf(ipath,"/MDF/images/%d",image_index);
01736         hid_t igrp=H5Gopen(file,ipath);
01737 
01738         if (igrp<0) {   //group not existed
01739                 throw _NotExistingObjectException(string(ipath));
01740         }
01741 
01742         string s("EMAN.");
01743         s += key;
01744         herr_t ret = H5Adelete(igrp, s.c_str());
01745 
01746         H5Gclose(igrp);
01747         delete imageio;
01748 
01749         if(ret >= 0) return 0;
01750         else return -1;
01751 }
01752 #endif  //EM_HDF5

Generated on Tue Jun 11 13:46:14 2013 for EMAN2 by  doxygen 1.3.9.1