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