EMAN2
pifio.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 "pifio.h"
00037 #include "portable_fileio.h"
00038 #include "geometry.h"
00039 #include "imageio.h"
00040 #include <cstring>
00041 
00042 #ifdef WIN32
00043 #include <time.h>
00044 #endif
00045 
00046 using namespace EMAN;
00047 
00048 PifIO::PifIO(const string & pif_filename, IOMode rw)
00049 :       filename(pif_filename), rw_mode(rw)
00050 {
00051         pif_file = 0;
00052         mode_size = 0;
00053         is_big_endian = ByteOrder::is_host_big_endian();
00054         initialized = false;
00055         real_scale_factor = 1;
00056         is_new_file = false;
00057         memset(&pfh, 0, sizeof(PifFileHeader));
00058 }
00059 
00060 PifIO::~PifIO()
00061 {
00062         if (pif_file) {
00063                 fclose(pif_file);
00064                 pif_file = 0;
00065         }
00066 }
00067 
00068 
00069 int PifIO::get_mode_size(PifDataMode mode)
00070 {
00071         int size = 0;
00072 
00073         switch (mode) {
00074         case PIF_CHAR:
00075         case PIF_BOXED_DATA:
00076                 size = sizeof(char);
00077                 break;
00078         case PIF_SHORT:
00079         case PIF_SHORT_FLOAT:
00080         case PIF_SHORT_COMPLEX:
00081         case PIF_SHORT_FLOAT_COMPLEX:
00082         case PIF_MAP_FLOAT_SHORT:
00083                 size = sizeof(short);
00084                 break;
00085         case PIF_FLOAT:
00086         case PIF_FLOAT_COMPLEX:
00087                 size = sizeof(float);
00088                 break;
00089         case PIF_FLOAT_INT:
00090         case PIF_FLOAT_INT_COMPLEX:
00091         case PIF_MAP_FLOAT_INT:
00092         case PIF_MAP_FLOAT_INT_2:
00093         case PIF_BOXED_FLOAT_INT:
00094                 size = sizeof(int);
00095                 break;
00096         default:
00097                 break;
00098         }
00099         return size;
00100 }
00101 
00102 bool PifIO::is_float_int(int m)
00103 {
00104         PifDataMode mode = static_cast < PifDataMode > (m);
00105         switch (mode) {
00106         case PIF_SHORT_FLOAT:
00107         case PIF_SHORT_FLOAT_COMPLEX:
00108         //case PIF_FLOAT:
00109         case PIF_FLOAT_INT:
00110         //case PIF_FLOAT_COMPLEX:
00111         case PIF_FLOAT_INT_COMPLEX:
00112         case PIF_MAP_FLOAT_SHORT:
00113         case PIF_MAP_FLOAT_INT:
00114         case PIF_MAP_FLOAT_INT_2:
00115         case PIF_BOXED_FLOAT_INT:
00116                 return true;
00117         default:
00118                 break;
00119         }
00120         return false;
00121 }
00122 
00123 void PifIO::init()
00124 {
00125         ENTERFUNC;
00126         if (initialized) {
00127                 return;
00128         }
00129 
00130         initialized = true;
00131         pif_file = sfopen(filename, rw_mode, &is_new_file);
00132 
00133         if (!is_new_file) {
00134                 if (fread(&pfh, sizeof(PifFileHeader), 1, pif_file) != 1) {
00135                         throw ImageReadException(filename, "PIF file header");
00136                 }
00137 
00138                 if (!is_valid(&pfh)) {
00139                         throw ImageReadException(filename, "invalid PIF file");
00140                 }
00141 
00142                 is_big_endian = ByteOrder::is_data_big_endian(&pfh.nz);
00143                 become_host_endian(&pfh.htype);
00144 
00145                 if (pfh.htype != 1) {
00146                         string desc = "only support PIF with all projects having the same dimensions";
00147                         throw ImageReadException(filename, desc);
00148                 }
00149 
00150                 become_host_endian(&pfh.mode);
00151                 become_host_endian(&pfh.nx);
00152                 become_host_endian(&pfh.ny);
00153                 become_host_endian(&pfh.nz);
00154                 become_host_endian(&pfh.nimg);
00155 
00156                 if (is_float_int(pfh.mode)) {
00157                         real_scale_factor = (float) atof(pfh.scalefactor);
00158                 }
00159 
00160                 mode_size = get_mode_size(static_cast < PifDataMode > (pfh.mode));
00161 
00162                 if (is_complex_mode()) {
00163                         pfh.nx *= 2;
00164                 }
00165         }
00166         EXITFUNC;
00167 }
00168 
00169 bool PifIO::is_valid(const void *first_block)
00170 {
00171         ENTERFUNC;
00172         bool result = false;
00173 
00174         if (first_block) {
00175                 const int *data = static_cast < const int *>(first_block);
00176                 int m1 = data[0];
00177                 int m2 = data[1];
00178                 int endian = data[7];
00179                 bool data_big_endian = false;
00180                 if (endian) {
00181                         data_big_endian = true;
00182                 }
00183 
00184                 if (data_big_endian != ByteOrder::is_host_big_endian()) {
00185                         ByteOrder::swap_bytes(&m1);
00186                         ByteOrder::swap_bytes(&m2);
00187                 }
00188 
00189                 if (m1 == PIF_MAGIC_NUM && m2 == PIF_MAGIC_NUM) {
00190                          result = true;
00191                 }
00192         }
00193 
00194         EXITFUNC;
00195         return result;
00196 }
00197 
00198 void PifIO::fseek_to(int image_index)
00199 {
00200         int pih_sz = sizeof(PifImageHeader);
00201         int image_size = 0;
00202 
00203 #if 0
00204         // this works for some images that PURDUE people gave to me.
00205         // But those images don't follow the PIF specification. So
00206         // I believe they are in wrong format.
00207         if (pfh.nimg == 1) {
00208                 image_size = pfh.nx * pfh.ny * pfh.nz;
00209         }
00210         else {
00211                 image_size = pfh.nx * pfh.ny;
00212         }
00213 #endif
00214         image_size = pfh.nx * pfh.ny * pfh.nz;
00215 
00216         size_t file_offset = sizeof(PifFileHeader) +
00217                 (pih_sz + image_size * mode_size) * image_index;
00218 
00219         portable_fseek(pif_file, file_offset, SEEK_SET);
00220 }
00221 
00222 
00223 int PifIO::read_header(Dict & dict, int image_index, const Region * area, bool)
00224 {
00225         ENTERFUNC;
00226 
00227         check_read_access(image_index);
00228         fseek_to(image_index);
00229 
00230         int pih_sz = sizeof(PifImageHeader);
00231         PifImageHeader pih;
00232 
00233         if (fread(&pih, pih_sz, 1, pif_file) != 1) {
00234                 throw ImageReadException(filename, "PIF Image header");
00235         }
00236         else {
00237                 check_region(area, FloatSize(pih.nx, pih.ny, pih.nz), is_new_file);
00238                 int xlen = 0, ylen = 0, zlen = 0;
00239                 EMUtil::get_region_dims(area, pih.nx, &xlen, pih.ny, &ylen, pih.nz, &zlen);
00240 
00241                 dict["nx"] = xlen;
00242                 dict["ny"] = ylen;
00243                 dict["nz"] = zlen;
00244 
00245                 dict["datatype"] = to_em_datatype(pih.mode);
00246 
00247                 dict["apix_x"] = static_cast < float >(pih.xlen);
00248                 dict["apix_y"] = static_cast < float >(pih.ylen);
00249                 dict["apix_z"] = static_cast < float >(pih.zlen);
00250 
00251                 dict["minimum"] = static_cast < float >(pih.min);
00252                 dict["maximum"] = static_cast < float >(pih.max);
00253                 dict["mean"] = static_cast < float >(pih.mean);
00254                 dict["sigma"] = static_cast < float >(pih.sigma);
00255 
00256                 dict["origin_x"] = static_cast < float >(pih.xorigin);
00257                 dict["origin_y"] = static_cast < float >(pih.yorigin);
00258         }
00259 
00260         EXITFUNC;
00261 
00262         return 0;
00263 }
00264 
00265 int PifIO::write_header(const Dict & dict, int image_index, const Region* area,
00266                                                 EMUtil::EMDataType, bool)
00267 {
00268         ENTERFUNC;
00269 
00270         if(image_index<0) {
00271                 image_index = get_nimg();
00272         }
00273 
00274         check_write_access(rw_mode, image_index);
00275 
00276         if (area) {
00277                 check_region(area, FloatSize(pfh.nx, pfh.ny, pfh.nz), is_new_file);
00278                 EXITFUNC;
00279                 return 0;
00280         }
00281 
00282         time_t t0 = time(0);
00283         struct tm *t = localtime(&t0);
00284 
00285         if (!is_new_file) {
00286                 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00287                         throw ImageWriteException(filename, "writing to opposite byteorder file");
00288                 }
00289 
00290                 int nx1 = dict["nx"];
00291                 int ny1 = dict["ny"];
00292                 int nz1 = dict["nz"];
00293 
00294                 int mode1 = to_pif_datatype(dict["datatype"]);
00295 
00296                 if (nx1 != pfh.nx || ny1 != pfh.ny || nz1 != pfh.nz) {
00297                         LOGERR("cannot write to different size file %s (%dx%dx%d != %dx%dx%d)",
00298                                    filename.c_str(), nx1, ny1, nz1, pfh.nx, pfh.ny, pfh.nz);
00299                         throw ImageWriteException(filename, "different image size");
00300                 }
00301 
00302                 if (mode1 != pfh.mode) {
00303                         throw ImageWriteException(filename, "different data type");
00304                 }
00305 
00306                 fseek_to(image_index);
00307         }
00308         else {
00309                 pfh.magic[0] = PIF_MAGIC_NUM;
00310                 pfh.magic[1] = PIF_MAGIC_NUM;
00311                 sprintf(pfh.scalefactor, "1.0");
00312 
00313                 pfh.mode = PIF_FLOAT_INT;
00314                 sprintf(pfh.program, "EMAN %d/%02d/%02d", t->tm_mon, t->tm_mday, t->tm_year);
00315 
00316                 pfh.htype = 1;
00317                 pfh.nx = dict["nx"];
00318                 pfh.ny = dict["ny"];
00319                 pfh.nz = dict["nz"];
00320                 pfh.nimg += 1;
00321                 pfh.endian = (int) ByteOrder::is_host_big_endian();
00322                 fwrite(&pfh, sizeof(PifFileHeader), 1, pif_file);
00323         }
00324 
00325         PifImageHeader pih;
00326         memset(&pih, 0, sizeof(PifImageHeader));
00327         pih.nx = dict["nx"];
00328         pih.ny = dict["ny"];
00329         pih.nz = dict["nz"];
00330 
00331         pih.mode = PIF_FLOAT;
00332         pih.xlen = dict["apix_x"];
00333         pih.ylen = dict["apix_y"];
00334         pih.zlen = dict["apix_z"];
00335         pih.alpha = 90;
00336         pih.beta = 90;
00337         pih.gamma = 90;
00338         pih.mapc = 1;
00339         pih.mapr = 2;
00340         pih.maps = 3;
00341         pih.min = dict["minimum"];
00342         pih.max = dict["maximum"];
00343         pih.mean = dict["mean"];
00344         pih.sigma = dict["sigma"];
00345 
00346         pih.xorigin = dict["origin_x"];
00347         pih.yorigin = dict["origin_y"];
00348 
00349         sprintf(pih.time, "%d/%02d/%02d %d:%02d",
00350                         t->tm_mon, t->tm_mday, t->tm_year, t->tm_hour, t->tm_min);
00351         fwrite(&pih, sizeof(PifImageHeader), 1, pif_file);
00352 
00353         EXITFUNC;
00354         return 0;
00355 }
00356 
00357 int PifIO::read_data(float *data, int image_index, const Region *area, bool)
00358 {
00359         ENTERFUNC;
00360 
00361         check_read_access(image_index, data);
00362         fseek_to(image_index);
00363 
00364         int pih_sz = sizeof(PifImageHeader);
00365         PifImageHeader pih;
00366 
00367         if (fread(&pih, pih_sz, 1, pif_file) != 1) {
00368                 throw ImageReadException(filename, "PIF Image header");
00369         }
00370 
00371         if (area) {
00372                 check_region(area, FloatSize(pih.nx, pih.ny, pih.nz), is_new_file);
00373         }
00374 
00375         PifDataMode data_mode = static_cast < PifDataMode > (pih.mode);
00376         int num_layers = pih.nz;
00377 #if 0
00378         if (pfh.nz == pfh.nimg) {
00379                 num_layers = 1;
00380         }
00381 #endif
00382         // new way to read PIF data. The new way includes region reading.
00383         // If it is tested to be OK, remove the code in #if 0 ... #endif
00384         unsigned char * cdata = (unsigned char*)data;
00385         short *sdata = (short*) data;
00386 
00387         EMUtil::process_region_io(cdata, pif_file, READ_ONLY,
00388                                                           0, mode_size, pih.nx, pih.ny,
00389                                                           num_layers, area);
00390 
00391         int xlen = 0, ylen = 0, zlen = 0;
00392         EMUtil::get_region_dims(area, pih.nx, &xlen, pih.ny, &ylen, pih.nz, &zlen);
00393         size_t size = xlen * ylen * zlen;
00394 
00395         if(data_mode == PIF_FLOAT || data_mode == PIF_FLOAT_COMPLEX)
00396         {
00397                 become_host_endian< float >(data, size);
00398         }
00399         else {
00400                 if (mode_size == sizeof(short)) {
00401                         become_host_endian((short *) sdata, size);
00402                 }
00403                 else if (mode_size == sizeof(int)) {
00404                         become_host_endian((int *) data, size);
00405                 }
00406 
00407                 if (mode_size == sizeof(char)) {
00408                         for (size_t i = 0; i < size; i++) {
00409                                 size_t j = size - 1 - i;
00410                                 data[j] = (float)(cdata[j]) * real_scale_factor;
00411                         }
00412                 }
00413                 else if (mode_size == sizeof(short)) {
00414                         for (size_t i = 0; i < size; i++) {
00415                                 size_t j = size - 1 - i;
00416                                 data[j] = (float)(sdata[j]) * real_scale_factor;
00417                         }
00418                 }
00419                 else if (mode_size == sizeof(int)) {
00420                         for (size_t i = 0; i < size; i++) {
00421                                 size_t j = size - 1 - i;
00422                                 data[j] = (float) ((int *)data)[j] * real_scale_factor;
00423                         }
00424                 }
00425         }
00426 
00427         // end of new way for region reading
00428 
00429 #if 0
00430         int buf_size = pih.nx * mode_size;
00431         unsigned char *buf = new unsigned char[buf_size];
00432 
00433         for (int l = 0; l < num_layers; l++) {
00434                 int offset1 = l * pfh.nx * pfh.ny;
00435                 for (int j = 0; j < pfh.ny; j++) {
00436                         if (fread(buf, mode_size, pfh.nx, pif_file) != (unsigned int) pfh.nx) {
00437                                 if( buf )
00438                                 {
00439                                         delete[]buf;
00440                                         buf = 0;
00441                                 }
00442                                 throw ImageReadException(filename, "incomplete PIF read");
00443                         }
00444 
00445                         if (mode_size == sizeof(short)) {
00446                                 become_host_endian((short *) buf, pfh.nx);
00447                         }
00448                         else if (mode_size == sizeof(int)) {
00449                                 become_host_endian((int *) buf, pfh.nx);
00450                         }
00451 
00452                         int offset2 = offset1 + j * pfh.nx;
00453 
00454                         for (int k = 0; k < pfh.nx; k++) {
00455                                 float curr_data = 0;
00456 
00457                                 if (mode_size == sizeof(char)) {
00458                                         curr_data = (float) buf[k];
00459                                 }
00460                                 else if (mode_size == sizeof(short)) {
00461                                         curr_data = (float) ((short *) buf)[k];
00462                                 }
00463                                 else if (mode_size == sizeof(int)) {
00464                                         curr_data = (float) ((int *) buf)[k];
00465                                 }
00466                                 data[offset2 + k] = curr_data * real_scale_factor;
00467                         }
00468                 }
00469         }
00470         if( buf )
00471         {
00472                 delete[]buf;
00473                 buf = 0;
00474         }
00475 #endif
00476         EXITFUNC;
00477         return 0;
00478 }
00479 
00480 
00481 int PifIO::write_data(float *data, int image_index, const Region* area,
00482                                           EMUtil::EMDataType, bool)
00483 {
00484         ENTERFUNC;
00485 
00486         check_write_access(rw_mode, image_index, 0, data);
00487         fseek_to(image_index);
00488 
00489         int nx = pfh.nx;
00490         int ny = pfh.ny;
00491         int nz = pfh.nz;
00492 
00493         check_region(area, FloatSize(nx, ny, nz), is_new_file);
00494 
00495         // to write a region; if it works, please remove the code
00496         // in #if 0 ... #endif
00497         EMUtil::process_region_io(data, pif_file, WRITE_ONLY, 0,
00498                                                           mode_size, nx, ny, nz, area);
00499 
00500 #if 0
00501         size_t idx;
00502         int *buf = new int[nx];
00503         for (int i = 0; i < nz; i++) {
00504                 for (int j = 0; j < ny; j++) {
00505                         for (int k = 0; k < pfh.nx; k++) {
00506                                 idx = i * nx * ny + j * nx + k;
00507                                 buf[k] = (int) data[idx];
00508                         }
00509                         fwrite(buf, sizeof(int) * nx, 1, pif_file);
00510                 }
00511         }
00512         if( buf )
00513         {
00514                 delete[]buf;
00515                 buf = 0;
00516         }
00517 #endif
00518 
00519         EXITFUNC;
00520         return 0;
00521 }
00522 
00523 void PifIO::flush()
00524 {
00525         fflush(pif_file);
00526 }
00527 
00528 bool PifIO::is_complex_mode()
00529 {
00530         init();
00531         if (pfh.mode == PIF_SHORT_COMPLEX ||
00532                 pfh.mode == PIF_FLOAT_INT_COMPLEX ||
00533                 pfh.mode == PIF_FLOAT_COMPLEX || pfh.mode == PIF_SHORT_FLOAT_COMPLEX) {
00534                 return true;
00535         }
00536         return false;
00537 }
00538 
00539 bool PifIO::is_image_big_endian()
00540 {
00541         init();
00542         return is_big_endian;
00543 }
00544 
00545 int PifIO::get_nimg()
00546 {
00547         init();
00548         return pfh.nimg;
00549 }
00550 
00551 
00552 int PifIO::to_em_datatype(int p)
00553 {
00554         PifDataMode mode = static_cast < PifDataMode > (p);
00555         EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
00556 
00557         switch (mode) {
00558         case PIF_CHAR:
00559         case PIF_BOXED_DATA:
00560                 e = EMUtil::EM_CHAR;
00561                 break;
00562 
00563         case PIF_SHORT:
00564         case PIF_SHORT_FLOAT:
00565         case PIF_MAP_FLOAT_SHORT:
00566                 e = EMUtil::EM_SHORT;
00567                 break;
00568 
00569         case PIF_SHORT_COMPLEX:
00570         case PIF_SHORT_FLOAT_COMPLEX:
00571                 e = EMUtil::EM_SHORT_COMPLEX;
00572                 break;
00573 
00574         case PIF_FLOAT:
00575         case PIF_FLOAT_INT:
00576         case PIF_MAP_FLOAT_INT:
00577         case PIF_MAP_FLOAT_INT_2:
00578         case PIF_BOXED_FLOAT_INT:
00579                 e = EMUtil::EM_FLOAT;
00580                 break;
00581         case PIF_FLOAT_COMPLEX:
00582         case PIF_FLOAT_INT_COMPLEX:
00583                 e = EMUtil::EM_FLOAT_COMPLEX;
00584                 break;
00585         case PIF_INVALID:
00586                 e = EMUtil::EM_UNKNOWN;
00587                 break;
00588         }
00589         return e;
00590 }
00591 
00592 int PifIO::to_pif_datatype(int e)
00593 {
00594         PifDataMode m = PIF_INVALID;
00595 
00596         switch (e) {
00597         case EMUtil::EM_CHAR:
00598                 m = PIF_BOXED_DATA;
00599                 break;
00600         case EMUtil::EM_SHORT:
00601                 m = PIF_SHORT;
00602                 break;
00603         case EMUtil::EM_SHORT_COMPLEX:
00604                 m = PIF_SHORT_COMPLEX;
00605                 break;
00606         case EMUtil::EM_FLOAT:
00607                 m = PIF_FLOAT_INT;
00608                 break;
00609         case EMUtil::EM_FLOAT_COMPLEX:
00610                 m = PIF_FLOAT_COMPLEX;
00611                 break;
00612         default:
00613                 LOGERR("unknown PIF mode: %d", e);
00614         }
00615 
00616         return m;
00617 }