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