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