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

Generated on Mon Mar 7 17:58:47 2011 for EMAN2 by  doxygen 1.4.7