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 "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
00109 case PIF_FLOAT_INT:
00110
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
00205
00206
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
00383
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
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
00496
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 }