vtkio.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 <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         //single image format, index can only be zero
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         //single image format, index can only be zero
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         //single image format, index can only be zero
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         //single image format, index can only be zero
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 

Generated on Tue Jun 11 12:40:27 2013 for EMAN2 by  doxygen 1.4.7