1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkMRCReader.cxx
5 
6   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7   All rights reserved.
8   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9 
10      This software is distributed WITHOUT ANY WARRANTY; without even
11      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12      PURPOSE.  See the above copyright notice for more information.
13 
14 =========================================================================*/
15 
16 #include "vtkMRCReader.h"
17 #include "vtkObjectFactory.h"
18 
19 #include "vtkByteSwap.h"
20 #include "vtkDataArray.h"
21 #include "vtkDataObject.h"
22 #include "vtkDataSetAttributes.h"
23 #include "vtkFloatArray.h"
24 #include "vtkImageData.h"
25 #include "vtkInformation.h"
26 #include "vtkInformationVector.h"
27 #include "vtkNew.h"
28 #include "vtkPointData.h"
29 #include "vtkSmartPointer.h"
30 #include "vtkStreamingDemandDrivenPipeline.h"
31 #include "vtkTypeInt16Array.h"
32 #include "vtkTypeInt8Array.h"
33 #include "vtkTypeUInt16Array.h"
34 #include "vtksys/FStream.hxx"
35 
36 #include <cassert>
37 
38 //#define VTK_DEBUG_MRC_HEADER
39 
40 namespace
41 {
42 //
43 // This struct is written based on the description found here:
44 // http://bio3d.colorado.edu/imod/doc/mrc_format.txt
45 //
46 struct mrc_file_header
47 {
48   vtkTypeInt32 nx;
49   vtkTypeInt32 ny;
50   vtkTypeInt32 nz;
51   vtkTypeInt32 mode;
52   vtkTypeInt32 nxstart;
53   vtkTypeInt32 nystart;
54   vtkTypeInt32 nzstart;
55   vtkTypeInt32 mx;
56   vtkTypeInt32 my;
57   vtkTypeInt32 mz;
58   vtkTypeFloat32 xlen;
59   vtkTypeFloat32 ylen;
60   vtkTypeFloat32 zlen;
61   vtkTypeFloat32 alpha;
62   vtkTypeFloat32 beta;
63   vtkTypeFloat32 gamma;
64   vtkTypeInt32 mapc;
65   vtkTypeInt32 mapr;
66   vtkTypeInt32 maps;
67   vtkTypeFloat32 amin;
68   vtkTypeFloat32 amax;
69   vtkTypeFloat32 amean;
70   vtkTypeInt32 ispg;
71   vtkTypeInt32 next;
72   vtkTypeInt16 creatid;
73   vtkTypeInt16 extra1[15];
74   vtkTypeInt16 nint;
75   vtkTypeInt16 nreal;
76   vtkTypeInt32 extra2[5];
77   vtkTypeInt32 imodStamp;
78   vtkTypeInt32 imodFlags;
79   vtkTypeInt16 idtype;
80   vtkTypeInt16 lens;
81   vtkTypeInt16 nd1;
82   vtkTypeInt16 nd2;
83   vtkTypeInt16 vd1;
84   vtkTypeInt16 vd2;
85   vtkTypeFloat32 tiltangles[6];
86   vtkTypeFloat32 xorg;
87   vtkTypeFloat32 yorg;
88   vtkTypeFloat32 zorg;
89   char cmap[4];
90   char stamp[4];
91   vtkTypeFloat32 rms;
92   vtkTypeInt32 nlabl;
93   char labl[10][80];
94 };
95 
96 #ifdef VTK_DEBUG_MRC_HEADER
operator <<(std::ostream & os,const mrc_file_header & hdr)97 std::ostream& operator<<(std::ostream& os, const mrc_file_header& hdr)
98 {
99   os << "extents:" << hdr.nx << " " << hdr.ny << " " << hdr.nz << std::endl;
100   os << "mode: " << hdr.mode << std::endl;
101   os << "start: " << hdr.nxstart << " " << hdr.nystart << " " << hdr.nzstart << std::endl;
102   os << "intervals: " << hdr.mx << " " << hdr.my << " " << hdr.mz << std::endl;
103   os << "len: " << hdr.xlen << " " << hdr.ylen << " " << hdr.zlen << std::endl;
104   os << "abg: " << hdr.alpha << " " << hdr.beta << " " << hdr.gamma << std::endl;
105   os << "map: " << hdr.mapc << " " << hdr.mapr << " " << hdr.maps << std::endl;
106   os << "min: " << hdr.amin << " max: " << hdr.amax << " mean: " << hdr.amean << std::endl;
107   os << "ispg: " << hdr.ispg << " next: " << hdr.next << std::endl;
108   // skipping extra for now
109   os << "nint: " << hdr.nint << " nreal: " << hdr.nreal << std::endl;
110   os << "imodStamp: " << hdr.imodStamp << " imodFlags: " << hdr.imodFlags << std::endl;
111   os << "idtype: " << hdr.idtype << " lens: " << hdr.lens << std::endl;
112   os << "nd1: " << hdr.nd1 << " nd2: " << hdr.nd2 << std::endl;
113   os << "vd1: " << hdr.vd1 << " vd2: " << hdr.vd2 << std::endl;
114   os << "tilt angles: " << hdr.tiltangles[0] << " " << hdr.tiltangles[1] << " " << hdr.tiltangles[2]
115      << " " << hdr.tiltangles[3] << " " << hdr.tiltangles[4] << " " << hdr.tiltangles[5]
116      << std::endl;
117   os << "org: " << hdr.xorg << " " << hdr.yorg << " " << hdr.zorg << std::endl;
118   os << "cmap: '" << hdr.cmap[0] << hdr.cmap[1] << hdr.cmap[2] << hdr.cmap[3] << "'";
119   os << " stamp: '" << hdr.stamp[0] << hdr.stamp[1] << hdr.stamp[2] << hdr.stamp[3] << "'"
120      << std::endl;
121   os << "rms: " << hdr.rms << " nlabl: " << hdr.nlabl << std::endl;
122   for (int i = 0; i < hdr.nlabl; ++i)
123   {
124     os.write(hdr.labl[i], 80);
125     os << std::endl;
126   }
127 }
128 #endif
129 }
130 
131 class vtkMRCReader::vtkInternal
132 {
133 public:
134   vtksys::ifstream* stream;
135   mrc_file_header header;
136 
vtkInternal()137   vtkInternal()
138     : stream(nullptr)
139     , header()
140   {
141   }
142 
~vtkInternal()143   ~vtkInternal() { delete stream; }
144 
openFile(const char * file)145   void openFile(const char* file)
146   {
147     delete stream;
148     stream = new vtksys::ifstream(file, std::ios::binary);
149   }
150 };
151 
152 vtkStandardNewMacro(vtkMRCReader);
153 
vtkMRCReader()154 vtkMRCReader::vtkMRCReader()
155 {
156   this->FileName = nullptr;
157   this->Internals = new vtkInternal;
158   this->SetNumberOfInputPorts(0);
159 }
160 
~vtkMRCReader()161 vtkMRCReader::~vtkMRCReader()
162 {
163   this->SetFileName(nullptr);
164   delete this->Internals;
165 }
166 
PrintSelf(ostream & os,vtkIndent indent)167 void vtkMRCReader::PrintSelf(ostream& os, vtkIndent indent)
168 {
169   Superclass::PrintSelf(os, indent);
170   os << indent << "FileName: " << (this->FileName ? this->FileName : "nullptr") << ", "
171      << std::endl;
172 }
173 
174 namespace
175 {
getFileDataType(int mode)176 int getFileDataType(int mode)
177 {
178   switch (mode)
179   {
180     case 0:
181     case 16:
182       return VTK_TYPE_UINT8;
183     case 2:
184     case 4:
185       return VTK_FLOAT;
186     case 1:
187     case 3:
188       return VTK_TYPE_INT16;
189     case 6:
190       return VTK_TYPE_UINT16;
191     default:
192       return -1;
193   }
194 }
195 
getFileDataNumComponents(int mode)196 int getFileDataNumComponents(int mode)
197 {
198   switch (mode)
199   {
200     case 0:
201     case 1:
202     case 2:
203       return 1;
204     case 3:
205     case 4:
206     case 6:
207       return 2;
208     case 16:
209       return 3;
210     default:
211       return -1;
212   }
213 }
214 }
215 
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)216 int vtkMRCReader::RequestInformation(vtkInformation* vtkNotUsed(request),
217   vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
218 {
219   // If this fails then the packing is wrong and the file's header will not
220   // be read in correctly
221   assert(sizeof(mrc_file_header) == 1024);
222   if (this->FileName)
223   {
224     this->Internals->openFile(this->FileName);
225     if (!this->Internals->stream)
226     {
227       vtkErrorMacro("Error opening input file");
228       return 0;
229     }
230     this->Internals->stream->read((char*)&this->Internals->header, sizeof(mrc_file_header));
231     if (this->Internals->header.stamp[0] != ((char)17))
232     // This is what the big-endian MRC files are supposed to look like.  I don't have one to
233     // test with though.  However, if it does not look like that, assume it is little endian.
234     // There are some non-conformant programs that don't correctly fill in this field, and
235     // assuming little endian is safer.
236     {
237       vtkByteSwap::Swap4LERange(&this->Internals->header, 24);
238       vtkByteSwap::Swap2LERange(&this->Internals->header.creatid, 1);
239       vtkByteSwap::Swap2LERange(&this->Internals->header.nint, 2);
240       vtkByteSwap::Swap4LERange(&this->Internals->header.imodStamp, 2);
241       vtkByteSwap::Swap2LERange(&this->Internals->header.idtype, 6);
242       vtkByteSwap::Swap4LERange(&this->Internals->header.tiltangles, 9);
243       vtkByteSwap::Swap4LERange(&this->Internals->header.rms, 2);
244     }
245     else
246     {
247       vtkByteSwap::Swap4BERange(&this->Internals->header, 24);
248       vtkByteSwap::Swap2BERange(&this->Internals->header.creatid, 1);
249       vtkByteSwap::Swap2BERange(&this->Internals->header.nint, 2);
250       vtkByteSwap::Swap4BERange(&this->Internals->header.imodStamp, 2);
251       vtkByteSwap::Swap2BERange(&this->Internals->header.idtype, 6);
252       vtkByteSwap::Swap4BERange(&this->Internals->header.tiltangles, 9);
253       vtkByteSwap::Swap4BERange(&this->Internals->header.rms, 2);
254     }
255 #ifdef VTK_DEBUG_MRC_HEADER
256     std::cout << this->Internals->header;
257 #endif
258     int extent[6];
259     extent[0] = this->Internals->header.nxstart;
260     extent[1] = extent[0] + this->Internals->header.nx - 1;
261     extent[2] = this->Internals->header.nystart;
262     extent[3] = extent[2] + this->Internals->header.ny - 1;
263     extent[4] = this->Internals->header.nzstart;
264     extent[5] = extent[4] + this->Internals->header.nz - 1;
265     double dataSpacing[3];
266     dataSpacing[0] = this->Internals->header.xlen / this->Internals->header.mx;
267     dataSpacing[1] = this->Internals->header.ylen / this->Internals->header.my;
268     dataSpacing[2] = this->Internals->header.zlen / this->Internals->header.mz;
269     double dataOrigin[3];
270     dataOrigin[0] = this->Internals->header.xorg;
271     dataOrigin[1] = this->Internals->header.yorg;
272     dataOrigin[2] = this->Internals->header.zorg;
273     vtkInformation* outInfo = outputVector->GetInformationObject(0);
274     outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent, 6);
275     outInfo->Set(vtkDataObject::SPACING(), dataSpacing, 3);
276     outInfo->Set(vtkDataObject::ORIGIN(), dataOrigin, 3);
277 
278     vtkDataObject::SetPointDataActiveScalarInfo(outInfo,
279       getFileDataType(this->Internals->header.mode),
280       getFileDataNumComponents(this->Internals->header.mode));
281 
282     outInfo->Set(vtkAlgorithm::CAN_PRODUCE_SUB_EXTENT(), 1);
283     return 1;
284   }
285   else
286   {
287     vtkErrorMacro("No input file set");
288     return 0;
289   }
290 }
291 
292 namespace
293 {
294 
295 typedef void (*ByteSwapFunction)(void*, size_t);
296 
getByteSwapFunction(int vtkType,bool isLittleEndian)297 ByteSwapFunction getByteSwapFunction(int vtkType, bool isLittleEndian)
298 {
299   int size = 0;
300   switch (vtkType)
301   {
302     vtkTemplateMacro(size = sizeof(VTK_TT));
303   }
304   if (size == 2)
305   {
306     return isLittleEndian ? &vtkByteSwap::Swap2LERange : &vtkByteSwap::Swap2BERange;
307   }
308   else if (size == 4)
309   {
310     return isLittleEndian ? &vtkByteSwap::Swap4LERange : &vtkByteSwap::Swap4BERange;
311   }
312   else if (size == 8)
313   {
314     return isLittleEndian ? &vtkByteSwap::Swap8LERange : &vtkByteSwap::Swap8BERange;
315   }
316   return nullptr;
317 }
318 
319 template <typename T>
readData(int numComponents,int * outExt,vtkIdType * outInc,vtkIdType * inOffsets,T * const outPtr,vtksys::ifstream & stream,vtkIdType dataStartPos,ByteSwapFunction byteSwapFunction)320 void readData(int numComponents, int* outExt, vtkIdType* outInc, vtkIdType* inOffsets,
321   T* const outPtr, vtksys::ifstream& stream, vtkIdType dataStartPos,
322   ByteSwapFunction byteSwapFunction)
323 {
324   vtkIdType lineSize = (outExt[1] - outExt[0] + 1) * numComponents;
325   T* ptr = outPtr;
326   for (vtkIdType z = outExt[4]; z <= outExt[5]; ++z)
327   {
328     for (vtkIdType y = outExt[2]; y <= outExt[3]; ++y)
329     {
330       assert(z >= 0 && y >= 0 && outExt[0] >= 0);
331       vtkIdType offset = z * inOffsets[2] + y * inOffsets[1] + outExt[0] * inOffsets[0];
332       offset = dataStartPos + offset * sizeof(T);
333 
334       stream.seekg(offset, vtksys::ifstream::beg);
335       // read the line
336       stream.read((char*)ptr, lineSize * sizeof(T));
337       if (byteSwapFunction)
338       {
339         byteSwapFunction((void*)outPtr, lineSize);
340       }
341       // update the data pointer
342       ptr += (lineSize + outInc[1]);
343     }
344     ptr += outInc[2];
345   }
346 }
347 
348 }
349 
ExecuteDataWithInformation(vtkDataObject * vtkNotUsed (output),vtkInformation * outInfo)350 void vtkMRCReader::ExecuteDataWithInformation(
351   vtkDataObject* vtkNotUsed(output), vtkInformation* outInfo)
352 {
353   vtkIdType outInc[3];
354   vtkIdType inOffsets[3];
355   int* outExt;
356   int modifiedOutExt[6];
357   int* execExt = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT());
358   vtkImageData* data = vtkImageData::GetData(outInfo);
359   this->AllocateOutputData(data, outInfo, execExt);
360 
361   if (data->GetNumberOfPoints() <= 0)
362   {
363     return;
364   }
365   outExt = data->GetExtent();
366   // this should result in the bottom corner of the image having extent
367   // 0,0,0 which makes the 'where in the file is this extent' math easier
368   modifiedOutExt[0] = outExt[0] - this->Internals->header.nxstart;
369   modifiedOutExt[1] = outExt[1] - this->Internals->header.nxstart;
370   modifiedOutExt[2] = outExt[2] - this->Internals->header.nystart;
371   modifiedOutExt[3] = outExt[3] - this->Internals->header.nystart;
372   modifiedOutExt[4] = outExt[4] - this->Internals->header.nzstart;
373   modifiedOutExt[5] = outExt[5] - this->Internals->header.nzstart;
374   data->GetContinuousIncrements(outExt, outInc[0], outInc[1], outInc[2]);
375   void* outPtr = data->GetScalarPointer(outExt[0], outExt[2], outExt[4]);
376 
377   if (!this->Internals->stream)
378   {
379     return;
380   }
381   // data start position is 1024 (the header size) plus the extended header size
382   vtkIdType dataStartPos = 1024 + this->Internals->header.next;
383   this->Internals->stream->seekg(dataStartPos, vtksys::ifstream::beg);
384 
385   int vtkType = getFileDataType(this->Internals->header.mode);
386   int numComponents = getFileDataNumComponents(this->Internals->header.mode);
387   inOffsets[0] = numComponents;
388   inOffsets[1] = this->Internals->header.nx * numComponents;
389   inOffsets[2] = this->Internals->header.ny * this->Internals->header.nx * numComponents;
390 
391   // This is what the big-endian MRC files are supposed to look like.  I don't have one to
392   // test with though.  However, if it does not look like that, assume it is little endian.
393   // There are some non-conformant programs that don't correctly fill in this field, and
394   // assuming little endian is safer.
395   bool fileIsLittleEndian = (this->Internals->header.stamp[0] != ((char)17));
396 
397   ByteSwapFunction byteSwapFunction = getByteSwapFunction(vtkType, fileIsLittleEndian);
398   switch (vtkType)
399   {
400     vtkTemplateMacro(readData<VTK_TT>(numComponents, modifiedOutExt, outInc, inOffsets,
401       static_cast<VTK_TT*>(outPtr), *this->Internals->stream, dataStartPos, byteSwapFunction));
402     default:
403       vtkErrorMacro("Unknown data type");
404   }
405 }
406