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 <cstring>
00037 #include "emio.h"
00038 #include "portable_fileio.h"
00039 #include "geometry.h"
00040
00041 using namespace EMAN;
00042
00043 EmIO::EmIO(const string & file, IOMode rw)
00044 : filename(file), rw_mode(rw), em_file(0), initialized(false)
00045 {
00046 mode_size = 0;
00047 mode = EM_EM_UNKNOWN;
00048 is_big_endian = ByteOrder::is_host_big_endian();
00049 is_new_file = false;
00050 memset(&emh, 0, sizeof(EMHeader));
00051 }
00052
00053 EmIO::~EmIO()
00054 {
00055 if (em_file) {
00056 fclose(em_file);
00057 em_file = 0;
00058 }
00059 }
00060
00061 void EmIO::init()
00062 {
00063 ENTERFUNC;
00064
00065 if (initialized) {
00066 return;
00067 }
00068
00069
00070 initialized = true;
00071 em_file = sfopen(filename, rw_mode, &is_new_file);
00072
00073 if (!is_new_file) {
00074 if (fread(&emh, sizeof(EMHeader), 1, em_file) != 1) {
00075 throw ImageReadException(filename, "EM header");
00076 }
00077 if (!is_valid(&emh)) {
00078 throw ImageReadException(filename, "invalid EM image");
00079 }
00080
00081 is_big_endian = ByteOrder::is_data_big_endian(&emh.nz);
00082 become_host_endian(&emh.nx);
00083 become_host_endian(&emh.ny);
00084 become_host_endian(&emh.nz);
00085
00086 mode = (DataType) emh.data_type;
00087
00088 if (mode == EM_EM_DOUBLE) {
00089 throw ImageReadException(filename, "DOUBLE data type not supported for EM image");
00090 }
00091
00092 mode_size = get_mode_size(emh.data_type);
00093 if (is_complex_mode()) {
00094 emh.nx *= 2;
00095 }
00096 }
00097 EXITFUNC;
00098 }
00099
00100 bool EmIO::is_valid(const void *first_block, off_t file_size)
00101 {
00102 ENTERFUNC;
00103
00104 if (!first_block) {
00105 return false;
00106 }
00107
00108 const char *data = static_cast < const char *>(first_block);
00109 char machine = data[0];
00110 char is_new_ver = data[1];
00111 char data_type = data[3];
00112
00113 const int *data1 = static_cast < const int *>(first_block);
00114 int nx = data1[1];
00115 int ny = data1[2];
00116 int nz = data1[3];
00117
00118 bool data_big_endian = ByteOrder::is_data_big_endian(&nz);
00119 if (data_big_endian != ByteOrder::is_host_big_endian()) {
00120 ByteOrder::swap_bytes(&nx);
00121 ByteOrder::swap_bytes(&ny);
00122 ByteOrder::swap_bytes(&nz);
00123 }
00124
00125 const int max_dim = 1 << 20;
00126
00127 if (((int) machine >= EM_OS8 && machine <= EM_PC) &&
00128 (is_new_ver == 0 || is_new_ver == 1) &&
00129 (data_type >= EM_EM_CHAR && data_type <= EM_EM_DOUBLE) &&
00130 (nx > 1 && nx < max_dim) && (ny > 0 && ny < max_dim) && (nz > 0 && nz < max_dim)) {
00131 if (file_size > 0) {
00132 off_t file_size1 = (off_t)nx * (off_t)ny * (off_t)nz * (off_t)get_mode_size(data_type) + (off_t)sizeof(EMHeader);
00133 if (file_size == file_size1) {
00134 return true;
00135 }
00136 }
00137 else {
00138 return true;
00139 }
00140 }
00141
00142 return false;
00143 }
00144
00145 int EmIO::read_header(Dict & dict, int image_index, const Region * area, bool)
00146 {
00147 ENTERFUNC;
00148
00149
00150 if(image_index == -1) {
00151 image_index = 0;
00152 }
00153
00154 if(image_index != 0) {
00155 throw ImageReadException(filename, "no stack allowed for MRC image. For take 2D slice out of 3D image, read the 3D image first, then use get_clip().");
00156 }
00157
00158 init();
00159 check_region(area, IntSize(emh.nx, emh.ny, emh.nz),false,false);
00160
00161 int xlen = 0, ylen = 0, zlen = 0;
00162 EMUtil::get_region_dims(area, emh.nx, &xlen, emh.ny, &ylen, emh.nz, &zlen);
00163
00164 dict["nx"] = xlen;
00165 dict["ny"] = ylen;
00166 dict["nz"] = zlen;
00167 dict["datatype"] = to_em_datatype(emh.data_type);
00168 EXITFUNC;
00169 return 0;
00170 }
00171
00172 int EmIO::write_header(const Dict & dict, int image_index, const Region* area,
00173 EMUtil::EMDataType, bool)
00174 {
00175 ENTERFUNC;
00176
00177 if(image_index == -1) {
00178 image_index = 0;
00179 }
00180 if(image_index != 0) {
00181 throw ImageWriteException(filename, "EM file does not support stack.");
00182 }
00183 check_write_access(rw_mode, image_index, 1);
00184 if (area) {
00185 check_region(area, FloatSize(emh.nx, emh.ny, emh.nz), is_new_file);
00186 EXITFUNC;
00187 return 0;
00188 }
00189
00190 emh.machine = static_cast < char >(get_machine_type());
00191 emh.nx = dict["nx"];
00192 emh.ny = dict["ny"];
00193 emh.nz = dict["nz"];
00194 emh.data_type = EM_EM_FLOAT;
00195
00196 rewind(em_file);
00197 if (fwrite(&emh, sizeof(EMHeader), 1, em_file) != 1) {
00198 throw ImageWriteException(filename, "EM Header");
00199 }
00200
00201 EXITFUNC;
00202 return 0;
00203 }
00204
00205 int EmIO::read_data(float *data, int image_index, const Region * area, bool)
00206 {
00207 ENTERFUNC;
00208
00209
00210 image_index = 0;
00211 check_read_access(image_index, data);
00212 check_region(area, IntSize(emh.nx, emh.ny, emh.nz),false,false);
00213
00214 portable_fseek(em_file, sizeof(EMHeader), SEEK_SET);
00215
00216 unsigned char *cdata = (unsigned char *) data;
00217 EMUtil::process_region_io(cdata, em_file, READ_ONLY, image_index, mode_size,
00218 emh.nx, emh.ny, emh.nz, area);
00219
00220 int xlen = 0, ylen = 0, zlen = 0;
00221 EMUtil::get_region_dims(area, emh.nx, &xlen, emh.ny, &ylen, emh.nz, &zlen);
00222
00223 int total_sz = xlen * ylen * zlen;
00224
00225 if (mode_size == sizeof(short)) {
00226 become_host_endian((short *) cdata, total_sz);
00227 }
00228 else if (mode_size == sizeof(int)) {
00229 become_host_endian((int *) cdata, total_sz);
00230 }
00231 else if (mode_size == sizeof(double)) {
00232 throw ImageReadException(filename, "double type image is not supported");
00233 }
00234
00235 for (int k = total_sz - 1; k >= 0; k--) {
00236 float curr_data = 0;
00237
00238 if (mode == EM_EM_CHAR) {
00239 curr_data = static_cast < float >(cdata[k]);
00240 }
00241 else if (mode == EM_EM_SHORT) {
00242 curr_data = static_cast < float >(((short *) cdata)[k]);
00243 }
00244 else if (mode == EM_EM_INT) {
00245 curr_data = static_cast < float >(((int *) cdata)[k]);
00246 }
00247 else if (mode == EM_EM_FLOAT || mode == EM_EM_COMPLEX) {
00248 curr_data = ((float *) cdata)[k];
00249 }
00250 else if (mode_size == sizeof(double)) {
00251 throw ImageReadException(filename, "double type image is not supported");
00252 }
00253
00254 data[k] = curr_data;
00255 }
00256
00257 EXITFUNC;
00258 return 0;
00259 }
00260
00261 int EmIO::write_data(float *data, int image_index, const Region* area, EMUtil::EMDataType, bool)
00262 {
00263 ENTERFUNC;
00264
00265 image_index = 0;
00266 check_write_access(rw_mode, image_index, 1, data);
00267 portable_fseek(em_file, sizeof(EMHeader), SEEK_SET);
00268
00269 EMUtil::process_region_io(data, em_file, rw_mode,
00270 image_index, sizeof(float),
00271 emh.nx, emh.ny, emh.nz, area);
00272 #if 0
00273 int sec_size = emh.nx * emh.ny;
00274 int row_size = sizeof(float) * emh.nx;
00275
00276 for (int i = 0; i < emh.nz; i++) {
00277 int k = i * sec_size;
00278 for (int j = 0; j < emh.ny; j++) {
00279 fwrite(&data[k + j * emh.nx], row_size, 1, em_file);
00280 }
00281 }
00282 #endif
00283
00284 EXITFUNC;
00285 return 0;
00286 }
00287
00288 void EmIO::flush()
00289 {
00290 fflush(em_file);
00291 }
00292
00293 bool EmIO::is_complex_mode()
00294 {
00295 init();
00296 if (emh.data_type == EM_EM_COMPLEX) {
00297 return true;
00298 }
00299 return false;
00300 }
00301
00302 bool EmIO::is_image_big_endian()
00303 {
00304 init();
00305 return is_big_endian;
00306 }
00307
00308
00309 int EmIO::to_em_datatype(char t)
00310 {
00311 DataType type = static_cast < DataType > (t);
00312 switch (type) {
00313 case EM_EM_CHAR:
00314 return EMUtil::EM_CHAR;
00315 case EM_EM_SHORT:
00316 return EMUtil::EM_SHORT;
00317 case EM_EM_INT:
00318 return EMUtil::EM_INT;
00319 case EM_EM_FLOAT:
00320 return EMUtil::EM_FLOAT;
00321 case EM_EM_DOUBLE:
00322 return EMUtil::EM_DOUBLE;
00323 case EM_EM_COMPLEX:
00324 return EMUtil::EM_FLOAT_COMPLEX;
00325 default:
00326 break;
00327 }
00328 return EMUtil::EM_UNKNOWN;
00329 }
00330
00331 int EmIO::get_machine_type()
00332 {
00333 int m = EM_UNKNOWN_MACHINE;
00334 #ifdef __sgi
00335 m = EM_SGI;
00336 #elif defined __linux__
00337 m = EM_PC;
00338 #elif defined __CYGWIN__
00339 m = EM_PC;
00340 #elif defined WIN32
00341 m = EM_PC;
00342 #elif defined MACOS
00343 m = EM_MAC;
00344 #elif defined macintosh
00345 m = EM_MAC;
00346 #elif defined __darwin__
00347 m = EM_MAC;
00348 #elif defined __APPLE__
00349 m = EM_MAC;
00350 #else
00351 m = EM_UNKNOWN_MACHINE;
00352 #endif
00353 return m;
00354 }
00355
00356 size_t EmIO::get_mode_size(char data_type)
00357 {
00358 int mode = (int) data_type;
00359 switch (mode) {
00360 case EM_EM_CHAR:
00361 return sizeof(char);
00362 case EM_EM_SHORT:
00363 return sizeof(short);
00364 case EM_EM_INT:
00365 case EM_EM_FLOAT:
00366 case EM_EM_COMPLEX:
00367 return sizeof(int);
00368 case EM_EM_DOUBLE:
00369 return sizeof(double);
00370 }
00371 return 0;
00372 }