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 <climits>
00038 #include "vtkio.h"
00039 #include "util.h"
00040 #include "geometry.h"
00041 #include "portable_fileio.h"
00042
00043 using namespace EMAN;
00044
00045 const char *VtkIO::MAGIC = "# vtk DataFile Version";
00046
00047 VtkIO::VtkIO(const string & vtk_filename, IOMode rw)
00048 : filename(vtk_filename), rw_mode(rw), vtk_file(0), initialized(false)
00049 {
00050 is_big_endian = ByteOrder::is_host_big_endian();
00051 is_new_file = false;
00052
00053 datatype = DATATYPE_UNKNOWN;
00054 filetype = VTK_UNKNOWN;
00055 nx = 0;
00056 ny = 0;
00057 nz = 0;
00058 originx = 0;
00059 originy = 0;
00060 originz = 0;
00061 spacingx = 0;
00062 spacingy = 0;
00063 spacingz = 0;
00064 file_offset = 0;
00065 }
00066
00067 VtkIO::~VtkIO()
00068 {
00069 if (vtk_file) {
00070 fclose(vtk_file);
00071 vtk_file = 0;
00072 }
00073 }
00074
00075 static int samestr(const char *s1, const char *s2)
00076 {
00077 return (strncmp(s1, s2, strlen(s2)) == 0);
00078 }
00079
00080
00081 void VtkIO::init()
00082 {
00083 if (initialized) {
00084 return;
00085 }
00086 ENTERFUNC;
00087 initialized = true;
00088
00089 vtk_file = sfopen(filename, rw_mode, &is_new_file);
00090
00091 if (!is_new_file) {
00092 char buf[1024];
00093 int bufsz = sizeof(buf);
00094 if (fgets(buf, bufsz, vtk_file) == 0) {
00095 throw ImageReadException(filename, "first block");
00096 }
00097
00098 if (!is_valid(&buf)) {
00099 throw ImageReadException(filename, "invalid VTK");
00100 }
00101
00102 if (fgets(buf, bufsz, vtk_file) == 0) {
00103 throw ImageReadException(filename, "read VTK file failed");
00104 }
00105
00106 if (fgets(buf, bufsz, vtk_file)) {
00107 if (samestr(buf, "ASCII")) {
00108 filetype = VTK_ASCII;
00109 }
00110 else if (samestr(buf, "BINARY")) {
00111 filetype = VTK_BINARY;
00112 }
00113 }
00114 else {
00115 throw ImageReadException(filename, "read VTK file failed");
00116 }
00117
00118 if (fgets(buf, bufsz, vtk_file)) {
00119 if (samestr(buf, "DATASET")) {
00120 char dataset_name[128];
00121 sscanf(buf, "DATASET %s", dataset_name);
00122 DatasetType ds_type = get_datasettype_from_name(dataset_name);
00123 read_dataset(ds_type);
00124 }
00125 }
00126 else {
00127 throw ImageReadException(filename, "read VTK file failed");
00128 }
00129
00130 while (fgets(buf, bufsz, vtk_file)) {
00131 if (samestr(buf, "SCALARS")) {
00132 char datatypestr[32];
00133 char scalartype[32];
00134 sscanf(buf, "SCALARS %s %s", scalartype, datatypestr);
00135
00136 datatype = get_datatype_from_name(datatypestr);
00137 if (datatype != UNSIGNED_SHORT && datatype != FLOAT) {
00138 string desc = "unknown data type: " + string(datatypestr);
00139 throw ImageReadException(filename, desc);
00140 }
00141 }
00142 else if (samestr(buf, "LOOKUP_TABLE")) {
00143 char tablename[128];
00144 sscanf(buf, "LOOKUP_TABLE %s", tablename);
00145 if (!samestr(tablename, "default")) {
00146 throw ImageReadException(filename, "only default LOOKUP_TABLE supported");
00147 }
00148 else {
00149 break;
00150 }
00151 }
00152 }
00153 #if 0
00154 if (filetype == VTK_BINARY) {
00155 throw ImageReadException(filename, "binary VTK is not supported");
00156 }
00157 #endif
00158 file_offset = portable_ftell(vtk_file);
00159 }
00160 EXITFUNC;
00161 }
00162
00163
00164 bool VtkIO::is_valid(const void *first_block)
00165 {
00166 ENTERFUNC;
00167 bool result = false;
00168 if (first_block) {
00169 result = Util::check_file_by_magic(first_block, MAGIC);
00170 }
00171 EXITFUNC;
00172 return result;
00173 }
00174
00175 int VtkIO::read_header(Dict & dict, int image_index, const Region * area, bool)
00176 {
00177 ENTERFUNC;
00178
00179
00180 if(image_index == -1) {
00181 image_index = 0;
00182 }
00183
00184 if(image_index != 0) {
00185 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().");
00186 }
00187
00188 init();
00189 check_region(area, IntSize(nx, ny, nz));
00190
00191 int xlen = 0, ylen = 0, zlen = 0;
00192 EMUtil::get_region_dims(area, nx, &xlen, ny, &ylen, nz, &zlen);
00193
00194 dict["nx"] = xlen;
00195 dict["ny"] = ylen;
00196 dict["nz"] = zlen;
00197
00198 dict["datatype"] = to_em_datatype(datatype);
00199
00200 dict["apix_x"] = spacingx;
00201 dict["apix_y"] = spacingy;
00202 dict["apix_z"] = spacingz;
00203
00204 dict["origin_x"] = originx;
00205 dict["origin_y"] = originy;
00206 dict["origin_z"] = originz;
00207
00208 EXITFUNC;
00209 return 0;
00210 }
00211
00212 int VtkIO::write_header(const Dict & dict, int image_index, const Region*,
00213 EMUtil::EMDataType, bool)
00214 {
00215 ENTERFUNC;
00216
00217
00218 if(image_index == -1) {
00219 image_index = 0;
00220 }
00221 if(image_index != 0) {
00222 throw ImageWriteException(filename, "VTK file does not support stack.");
00223 }
00224 check_write_access(rw_mode, image_index);
00225
00226 nx = dict["nx"];
00227 ny = dict["ny"];
00228 nz = dict["nz"];
00229
00230 originx = dict["origin_x"];
00231 originy = dict["origin_y"];
00232 originz = dict["origin_z"];
00233
00234 spacingx = dict["apix_x"];
00235 spacingy = dict["apix_y"];
00236 spacingz = dict["apix_z"];
00237
00238 fprintf(vtk_file, "# vtk DataFile Version 2.0\n");
00239 fprintf(vtk_file, "EMAN\n");
00240 fprintf(vtk_file, "BINARY\n");
00241 fprintf(vtk_file, "DATASET STRUCTURED_POINTS\n");
00242 fprintf(vtk_file, "DIMENSIONS %0d %0d %0d\nORIGIN %f %f %f\nSPACING %f %f %f\n",
00243 nx, ny, nz, originx, originy, originz, spacingx, spacingy, spacingz);
00244
00245
00246 fprintf(vtk_file, "POINT_DATA %0lu\nSCALARS density float 1\nLOOKUP_TABLE default\n",
00247 (size_t)nx * ny * nz);
00248 EXITFUNC;
00249 return 0;
00250 }
00251
00252 int VtkIO::read_data(float *data, int image_index, const Region * area, bool)
00253 {
00254 ENTERFUNC;
00255
00256
00257 image_index = 0;
00258 check_read_access(image_index, data);
00259
00260 if (area) {
00261 LOGWARN("read VTK region is not supported yet. Read whole image instead.");
00262 }
00263
00264 portable_fseek(vtk_file, file_offset, SEEK_SET);
00265
00266 int xlen = 0, ylen = 0, zlen = 0;
00267 int x0 = 0, y0 = 0, z0 = 0;
00268 EMUtil::get_region_dims(area, nx, &xlen, ny, &ylen, nz, &zlen);
00269 EMUtil::get_region_origins(area, &x0, &y0, &z0, nz, image_index);
00270
00271 if (filetype == VTK_ASCII) {
00272
00273 int bufsz = nx * get_mode_size(datatype) * CHAR_BIT;
00274 char *buf = new char[bufsz];
00275 int i = 0;
00276
00277 while (fgets(buf, bufsz, vtk_file)) {
00278 size_t bufslen = strlen(buf) - 1;
00279 char numstr[32];
00280 int k = 0;
00281 for (size_t j = 0; j < bufslen; j++) {
00282 if (!isspace(buf[j])) {
00283 numstr[k++] = buf[j];
00284 }
00285 else {
00286 numstr[k] = '\0';
00287 data[i++] = (float)atoi(numstr);
00288 k = 0;
00289 }
00290 }
00291 }
00292 if( buf )
00293 {
00294 delete[]buf;
00295 buf = 0;
00296 }
00297 }
00298 else if (filetype == VTK_BINARY) {
00299 int nxy = nx * ny;
00300 int row_size = nx * get_mode_size(datatype);
00301
00302 for (int i = 0; i < nz; i++) {
00303 int i2 = i * nxy;
00304 for (int j = 0; j < ny; j++) {
00305 fread(&data[i2 + j * nx], row_size, 1, vtk_file);
00306 }
00307 }
00308
00309 if (!ByteOrder::is_host_big_endian()) {
00310 ByteOrder::swap_bytes(data, (size_t)nx * ny * nz);
00311 }
00312 }
00313
00314 EXITFUNC;
00315 return 0;
00316 }
00317
00318 int VtkIO::write_data(float *data, int image_index, const Region* ,
00319 EMUtil::EMDataType, bool)
00320 {
00321 ENTERFUNC;
00322
00323
00324 image_index = 0;
00325 check_write_access(rw_mode, image_index, 1, data);
00326
00327 bool swapped = false;
00328 if (!ByteOrder::is_host_big_endian()) {
00329 ByteOrder::swap_bytes(data, (size_t)nx * ny * nz);
00330 swapped = true;
00331 }
00332
00333 fwrite(data, nx * nz, ny * sizeof(float), vtk_file);
00334
00335 if (swapped) {
00336 ByteOrder::swap_bytes(data, (size_t)nx * ny * nz);
00337 }
00338 EXITFUNC;
00339 return 0;
00340 }
00341
00342 void VtkIO::flush()
00343 {
00344 fflush(vtk_file);
00345 }
00346
00347 bool VtkIO::is_complex_mode()
00348 {
00349 return false;
00350 }
00351
00352 bool VtkIO::is_image_big_endian()
00353 {
00354 return true;
00355 }
00356
00357 int VtkIO::to_em_datatype(int vtk_datatype)
00358 {
00359 DataType d = static_cast < DataType > (vtk_datatype);
00360 switch (d) {
00361 case UNSIGNED_SHORT:
00362 return EMUtil::EM_USHORT;
00363 case FLOAT:
00364 return EMUtil::EM_FLOAT;
00365 default:
00366 break;
00367 }
00368 return EMUtil::EM_UNKNOWN;
00369 }
00370
00371
00372 int VtkIO::get_mode_size(DataType d)
00373 {
00374 switch (d) {
00375 case UNSIGNED_CHAR:
00376 case CHAR:
00377 return sizeof(char);
00378 case UNSIGNED_SHORT:
00379 case SHORT:
00380 return sizeof(short);
00381 case UNSIGNED_INT:
00382 case INT:
00383 return sizeof(int);
00384 case UNSIGNED_LONG:
00385 case LONG:
00386 return sizeof(long);
00387 case FLOAT:
00388 return sizeof(float);
00389 case DOUBLE:
00390 return sizeof(double);
00391 default:
00392 LOGERR("don't support this data type '%d'", d);
00393 break;
00394 }
00395 return 0;
00396 }
00397
00398 VtkIO::DataType VtkIO::get_datatype_from_name(const string& datatype_name)
00399 {
00400 static bool initialized = false;
00401 static map < string, VtkIO::DataType > datatypes;
00402
00403 if (!initialized) {
00404 datatypes["bit"] = BIT;
00405
00406 datatypes["unsigned_char"] = UNSIGNED_CHAR;
00407 datatypes["char"] = CHAR;
00408
00409 datatypes["unsigned_short"] = UNSIGNED_SHORT;
00410 datatypes["short"] = SHORT;
00411
00412 datatypes["unsigned_int"] = UNSIGNED_INT;
00413 datatypes["int"] = INT;
00414
00415 datatypes["unsigned_long"] = UNSIGNED_LONG;
00416 datatypes["long"] = LONG;
00417
00418 datatypes["float"] = FLOAT;
00419 datatypes["double"] = DOUBLE;
00420 initialized = true;
00421 }
00422
00423 DataType result = DATATYPE_UNKNOWN;
00424
00425 if (datatypes.find(datatype_name) != datatypes.end()) {
00426 result = datatypes[datatype_name];
00427 }
00428 return result;
00429 }
00430
00431 VtkIO::DatasetType VtkIO::get_datasettype_from_name(const string& dataset_name)
00432 {
00433
00434 static bool initialized = false;
00435 static map < string, DatasetType > types;
00436
00437 if (!initialized) {
00438 types["STRUCTURED_POINTS"] = STRUCTURED_POINTS;
00439 types["STRUCTURED_GRID"] = STRUCTURED_GRID;
00440 types["RECTILINEAR_GRID"] = RECTILINEAR_GRID;
00441 types["UNSTRUCTURED_GRID"] = UNSTRUCTURED_GRID;
00442 types["POLYDATA"] = POLYDATA;
00443 }
00444
00445 DatasetType result = DATASET_UNKNOWN;
00446 if (types.find(dataset_name) != types.end()) {
00447 result = types[dataset_name];
00448 }
00449 return result;
00450 }
00451
00452 void VtkIO::read_dataset(DatasetType dstype)
00453 {
00454 char buf[1024];
00455 int bufsz = sizeof(buf);
00456
00457 if (dstype == STRUCTURED_POINTS) {
00458 int nlines = 3;
00459 int i = 0;
00460 while (i < nlines && fgets(buf, bufsz, vtk_file)) {
00461 if (samestr(buf, "DIMENSIONS")) {
00462 sscanf(buf, "DIMENSIONS %d %d %d", &nx, &ny, &nz);
00463 }
00464 else if (samestr(buf, "ORIGIN")) {
00465 sscanf(buf, "ORIGIN %f %f %f", &originx, &originy, &originz);
00466 }
00467 else if (samestr(buf, "SPACING") || samestr(buf, "ASPECT_RATIO")) {
00468 if (samestr(buf, "SPACING")) {
00469 sscanf(buf, "SPACING %f %f %f", &spacingx, &spacingy, &spacingz);
00470 }
00471 else {
00472 sscanf(buf, "ASPECT_RATIO %f %f %f", &spacingx, &spacingy, &spacingz);
00473 }
00474
00475 if (spacingx != spacingy || spacingx != spacingz || spacingy != spacingz) {
00476 throw ImageReadException(filename,
00477 "not support non-uniform spacing VTK so far\n");
00478 }
00479 }
00480 i++;
00481 }
00482
00483 if (i != nlines) {
00484 throw ImageReadException(filename, "read VTK file failed");
00485 }
00486 }
00487 else {
00488 throw ImageReadException(filename, "only STRUCTURED_POINTS is supported so far");
00489 }
00490 }
00491