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

Generated on Fri Aug 10 16:35:12 2012 for EMAN2 by  doxygen 1.3.9.1