1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkVolume16Reader.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 #include "vtkVolume16Reader.h"
16
17 #include "vtkImageData.h"
18 #include "vtkInformation.h"
19 #include "vtkInformationVector.h"
20 #include "vtkObjectFactory.h"
21 #include "vtkPointData.h"
22 #include "vtkStreamingDemandDrivenPipeline.h"
23 #include "vtkTransform.h"
24 #include "vtkUnsignedShortArray.h"
25
26 vtkStandardNewMacro(vtkVolume16Reader);
27
28 vtkCxxSetObjectMacro(vtkVolume16Reader,Transform,vtkTransform);
29
30 //----------------------------------------------------------------------------
31 // Construct object with NULL file prefix; file pattern "%s.%d"; image range
32 // set to (1,1); data origin (0,0,0); data spacing (1,1,1); no data mask;
33 // header size 0; and byte swapping turned off.
vtkVolume16Reader()34 vtkVolume16Reader::vtkVolume16Reader()
35 {
36 this->DataMask = 0x0000;
37 this->HeaderSize = 0;
38 this->SwapBytes = 0;
39 this->DataDimensions[0] = this->DataDimensions[1] = 0;
40 this->Transform = NULL;
41 }
42
43 //----------------------------------------------------------------------------
~vtkVolume16Reader()44 vtkVolume16Reader::~vtkVolume16Reader()
45 {
46 this->SetTransform(NULL);
47 }
48
49 //----------------------------------------------------------------------------
SetDataByteOrderToBigEndian()50 void vtkVolume16Reader::SetDataByteOrderToBigEndian()
51 {
52 #ifndef VTK_WORDS_BIGENDIAN
53 this->SwapBytesOn();
54 #else
55 this->SwapBytesOff();
56 #endif
57 }
58
59 //----------------------------------------------------------------------------
SetDataByteOrderToLittleEndian()60 void vtkVolume16Reader::SetDataByteOrderToLittleEndian()
61 {
62 #ifdef VTK_WORDS_BIGENDIAN
63 this->SwapBytesOn();
64 #else
65 this->SwapBytesOff();
66 #endif
67 }
68
69 //----------------------------------------------------------------------------
SetDataByteOrder(int byteOrder)70 void vtkVolume16Reader::SetDataByteOrder(int byteOrder)
71 {
72 if ( byteOrder == VTK_FILE_BYTE_ORDER_BIG_ENDIAN )
73 {
74 this->SetDataByteOrderToBigEndian();
75 }
76 else
77 {
78 this->SetDataByteOrderToLittleEndian();
79 }
80 }
81
82 //----------------------------------------------------------------------------
GetDataByteOrder()83 int vtkVolume16Reader::GetDataByteOrder()
84 {
85 #ifdef VTK_WORDS_BIGENDIAN
86 if ( this->SwapBytes )
87 {
88 return VTK_FILE_BYTE_ORDER_LITTLE_ENDIAN;
89 }
90 else
91 {
92 return VTK_FILE_BYTE_ORDER_BIG_ENDIAN;
93 }
94 #else
95 if ( this->SwapBytes )
96 {
97 return VTK_FILE_BYTE_ORDER_BIG_ENDIAN;
98 }
99 else
100 {
101 return VTK_FILE_BYTE_ORDER_LITTLE_ENDIAN;
102 }
103 #endif
104 }
105
106 //----------------------------------------------------------------------------
GetDataByteOrderAsString()107 const char *vtkVolume16Reader::GetDataByteOrderAsString()
108 {
109 #ifdef VTK_WORDS_BIGENDIAN
110 if ( this->SwapBytes )
111 {
112 return "LittleEndian";
113 }
114 else
115 {
116 return "BigEndian";
117 }
118 #else
119 if ( this->SwapBytes )
120 {
121 return "BigEndian";
122 }
123 else
124 {
125 return "LittleEndian";
126 }
127 #endif
128 }
129
130
131 //----------------------------------------------------------------------------
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)132 int vtkVolume16Reader::RequestInformation(
133 vtkInformation *vtkNotUsed(request),
134 vtkInformationVector **vtkNotUsed(inputVector),
135 vtkInformationVector *outputVector)
136 {
137 int dim[3];
138
139 this->ComputeTransformedDimensions(dim);
140
141 vtkInformation *outInfo = outputVector->GetInformationObject(0);
142 outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
143 0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1);
144
145 vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_UNSIGNED_SHORT, 1);
146 outInfo->Set(vtkDataObject::SPACING(), this->DataSpacing, 3);
147 outInfo->Set(vtkDataObject::ORIGIN(), this->DataOrigin, 3);
148 // spacing and origin ?
149
150 return 1;
151 }
152
153 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)154 int vtkVolume16Reader::RequestData(
155 vtkInformation *vtkNotUsed(request),
156 vtkInformationVector **vtkNotUsed(inputVector),
157 vtkInformationVector *outputVector)
158 {
159 int first, last;
160 int *dim;
161 int dimensions[3];
162 double Spacing[3];
163 double origin[3];
164
165 vtkInformation* outInfo = outputVector->GetInformationObject(0);
166 vtkDataObject* output_do = outInfo->Get(vtkDataObject::DATA_OBJECT());
167 vtkImageData *output = this->AllocateOutputData(output_do,
168 outInfo);
169 vtkUnsignedShortArray *newScalars =
170 vtkUnsignedShortArray::SafeDownCast(output->GetPointData()->GetScalars());
171
172 // Validate instance variables
173 if (this->FilePrefix == NULL)
174 {
175 vtkErrorMacro(<< "FilePrefix is NULL");
176 return 1;
177 }
178
179 if (this->HeaderSize < 0)
180 {
181 vtkErrorMacro(<< "HeaderSize " << this->HeaderSize << " must be >= 0");
182 return 1;
183 }
184
185 dim = this->DataDimensions;
186
187 if (dim[0] <= 0 || dim[1] <= 0)
188 {
189 vtkErrorMacro(<< "x, y dimensions " << dim[0] << ", " << dim[1]
190 << "must be greater than 0.");
191 return 1;
192 }
193
194 if ( (this->ImageRange[1]-this->ImageRange[0]) <= 0 )
195 {
196 this->ReadImage(this->ImageRange[0], newScalars);
197 }
198 else
199 {
200 first = this->ImageRange[0];
201 last = this->ImageRange[1];
202 this->ReadVolume(first, last, newScalars);
203 }
204
205 // calculate dimensions of output from data dimensions and transform
206 this->ComputeTransformedDimensions (dimensions);
207 output->SetDimensions(dimensions);
208
209 // calculate spacing of output from data spacing and transform
210 this->ComputeTransformedSpacing(Spacing);
211
212 // calculate origin of output from data origin and transform
213 this->ComputeTransformedOrigin(origin);
214
215 // adjust spacing and origin if spacing is negative
216 this->AdjustSpacingAndOrigin(dimensions, Spacing, origin);
217
218 output->SetSpacing(Spacing);
219 output->SetOrigin(origin);
220
221 return 1;
222 }
223
224 //----------------------------------------------------------------------------
GetImage(int ImageNumber)225 vtkImageData *vtkVolume16Reader::GetImage(int ImageNumber)
226 {
227 vtkUnsignedShortArray *newScalars;
228 int *dim;
229 int dimensions[3];
230 vtkImageData *result;
231
232 // Validate instance variables
233 if (this->FilePrefix == NULL)
234 {
235 vtkErrorMacro(<< "FilePrefix is NULL");
236 return NULL;
237 }
238
239 if (this->HeaderSize < 0)
240 {
241 vtkErrorMacro(<< "HeaderSize " << this->HeaderSize << " must be >= 0");
242 return NULL;
243 }
244
245 dim = this->DataDimensions;
246
247 if (dim[0] <= 0 || dim[1] <= 0)
248 {
249 vtkErrorMacro(<< "x, y dimensions " << dim[0] << ", " << dim[1]
250 << "must be greater than 0.");
251 return NULL;
252 }
253
254 result = vtkImageData::New();
255 newScalars = vtkUnsignedShortArray::New();
256 this->ReadImage(ImageNumber, newScalars);
257 dimensions[0] = dim[0]; dimensions[1] = dim[1];
258 dimensions[2] = 1;
259 result->SetDimensions(dimensions);
260 result->SetSpacing(this->DataSpacing);
261 result->SetOrigin(this->DataOrigin);
262 if ( newScalars )
263 {
264 result->GetPointData()->SetScalars(newScalars);
265 newScalars->Delete();
266 }
267 return result;
268 }
269
270 //----------------------------------------------------------------------------
271 // Read a slice of volume data.
ReadImage(int sliceNumber,vtkUnsignedShortArray * scalars)272 void vtkVolume16Reader::ReadImage(int sliceNumber,
273 vtkUnsignedShortArray *scalars)
274 {
275 unsigned short *pixels;
276 FILE *fp;
277 int numPts;
278 char filename[1024];
279
280 // build the file name. if there is no prefix, just use the slice number
281 if (this->FilePrefix)
282 {
283 sprintf (filename, this->FilePattern, this->FilePrefix, sliceNumber);
284 }
285 else
286 {
287 sprintf (filename, this->FilePattern, sliceNumber);
288 }
289 if ( !(fp = fopen(filename,"rb")) )
290 {
291 vtkErrorMacro(<<"Can't open file: " << filename);
292 return;
293 }
294
295 numPts = this->DataDimensions[0] * this->DataDimensions[1];
296
297 // get a pointer to the data
298 pixels = scalars->WritePointer(0, numPts);
299
300 // read the image data
301 this->Read16BitImage (fp, pixels, this->DataDimensions[0],
302 this->DataDimensions[1], this->HeaderSize,
303 this->SwapBytes);
304
305 // close the file
306 fclose (fp);
307 }
308
309 //----------------------------------------------------------------------------
310 // Read a volume of data.
ReadVolume(int first,int last,vtkUnsignedShortArray * scalars)311 void vtkVolume16Reader::ReadVolume(int first, int last,
312 vtkUnsignedShortArray *scalars)
313 {
314 unsigned short *pixels;
315 unsigned short *slice;
316 FILE *fp;
317 int numPts;
318 int fileNumber;
319 int status;
320 int numberSlices = last - first + 1;
321 char filename[1024];
322 int dimensions[3];
323 int bounds[6];
324
325 // calculate the number of points per image
326 numPts = this->DataDimensions[0] * this->DataDimensions[1];
327
328 // compute transformed dimensions
329 this->ComputeTransformedDimensions (dimensions);
330
331 // compute transformed bounds
332 this->ComputeTransformedBounds (bounds);
333
334 // get memory for slice
335 slice = new unsigned short[numPts];
336
337 // get a pointer to the scalar data
338 pixels = scalars->WritePointer(0, numPts*numberSlices);
339
340 vtkDebugMacro (<< "Creating scalars with " << numPts * numberSlices
341 << " points.");
342
343 // build each file name and read the data from the file
344 for (fileNumber = first; fileNumber <= last; fileNumber++)
345 {
346 // build the file name. if there is no prefix, just use the slice number
347 if (this->FilePattern)
348 {
349 sprintf (filename, this->FilePattern, this->FilePrefix, fileNumber);
350 }
351 else
352 {
353 sprintf (filename, this->FilePattern, fileNumber);
354 }
355 if ( !(fp = fopen(filename,"rb")) )
356 {
357 vtkErrorMacro(<<"Can't find file: " << filename);
358 return;
359 }
360
361 vtkDebugMacro ( << "Reading " << filename );
362
363 // read the image data
364 status = this->Read16BitImage (fp, slice, this->DataDimensions[0],
365 this->DataDimensions[1], this->HeaderSize, this->SwapBytes);
366
367 fclose (fp);
368
369 if (status == 0)
370 {
371 break;
372 }
373
374 // transform slice
375 this->TransformSlice (slice, pixels, fileNumber - first, dimensions, bounds);
376 }
377
378 delete []slice;
379 }
380
381 //----------------------------------------------------------------------------
Read16BitImage(FILE * fp,unsigned short * pixels,int xsize,int ysize,int skip,int swapBytes)382 int vtkVolume16Reader:: Read16BitImage (FILE *fp, unsigned short *pixels, int xsize,
383 int ysize, int skip, int swapBytes)
384 {
385 unsigned short *shortPtr;
386 int numShorts = xsize * ysize;
387
388 if (skip)
389 {
390 fseek (fp, skip, 0);
391 }
392
393 shortPtr = pixels;
394 shortPtr += xsize*(ysize - 1);
395 for (int j=0; j<ysize; j++, shortPtr -= xsize)
396 {
397 if ( ! fread(shortPtr,sizeof (unsigned short),xsize,fp) )
398 {
399 vtkErrorMacro(<<"Error reading raw pgm data!");
400 return 0;
401 }
402 }
403
404 if (swapBytes)
405 {
406 unsigned char *bytes = (unsigned char *) pixels;
407 unsigned char tmp;
408 int i;
409 for (i = 0; i < numShorts; i++, bytes += 2)
410 {
411 tmp = *bytes;
412 *bytes = *(bytes + 1);
413 *(bytes + 1) = tmp;
414 }
415 }
416
417 if (this->DataMask != 0x0000 )
418 {
419 unsigned short *dataPtr = pixels;
420 int i;
421 for (i = 0; i < numShorts; i++, dataPtr++)
422 {
423 *dataPtr &= this->DataMask;
424 }
425 }
426
427 return 1;
428 }
429
430 //----------------------------------------------------------------------------
ComputeTransformedSpacing(double spacing[3])431 void vtkVolume16Reader::ComputeTransformedSpacing (double spacing[3])
432 {
433 if (!this->Transform)
434 {
435 memcpy (spacing, this->DataSpacing, 3 * sizeof (double));
436 }
437 else
438 {
439 double transformedSpacing[4];
440 memcpy (transformedSpacing, this->DataSpacing, 3 * sizeof (double));
441 transformedSpacing[3] = 1.0;
442 this->Transform->MultiplyPoint (transformedSpacing, transformedSpacing);
443
444 for (int i = 0; i < 3; i++)
445 {
446 spacing[i] = transformedSpacing[i];
447 }
448 vtkDebugMacro("Transformed Spacing " <<
449 spacing[0] << ", " << spacing[1] << ", " << spacing[2]);
450 }
451 }
452
453 //----------------------------------------------------------------------------
ComputeTransformedOrigin(double origin[3])454 void vtkVolume16Reader::ComputeTransformedOrigin (double origin[3])
455 {
456 if (!this->Transform)
457 {
458 memcpy (origin, this->DataOrigin, 3 * sizeof (double));
459 }
460 else
461 {
462 double transformedOrigin[4];
463 memcpy (transformedOrigin, this->DataOrigin, 3 * sizeof (double));
464 transformedOrigin[3] = 1.0;
465 this->Transform->MultiplyPoint (transformedOrigin, transformedOrigin);
466
467 for (int i = 0; i < 3; i++)
468 {
469 origin[i] = transformedOrigin[i];
470 }
471 vtkDebugMacro("Transformed Origin " <<
472 origin[0] << ", " << origin[1] << ", " << origin[2]);
473 }
474 }
475
476 //----------------------------------------------------------------------------
ComputeTransformedDimensions(int dimensions[3])477 void vtkVolume16Reader::ComputeTransformedDimensions (int dimensions[3])
478 {
479 double transformedDimensions[4];
480 if (!this->Transform)
481 {
482 dimensions[0] = this->DataDimensions[0];
483 dimensions[1] = this->DataDimensions[1];
484 dimensions[2] = this->ImageRange[1] - this->ImageRange[0] + 1;
485 }
486 else
487 {
488 transformedDimensions[0] = this->DataDimensions[0];
489 transformedDimensions[1] = this->DataDimensions[1];
490 transformedDimensions[2] = this->ImageRange[1] - this->ImageRange[0] + 1;
491 transformedDimensions[3] = 1.0;
492 this->Transform->MultiplyPoint (transformedDimensions, transformedDimensions);
493 dimensions[0] = (int) transformedDimensions[0];
494 dimensions[1] = (int) transformedDimensions[1];
495 dimensions[2] = (int) transformedDimensions[2];
496 if (dimensions[0] < 0)
497 {
498 dimensions[0] = -dimensions[0];
499 }
500 if (dimensions[1] < 0)
501 {
502 dimensions[1] = -dimensions[1];
503 }
504 if (dimensions[2] < 0)
505 {
506 dimensions[2] = -dimensions[2];
507 }
508 vtkDebugMacro(<< "Transformed dimensions are:" << dimensions[0] << ", "
509 << dimensions[1] << ", "
510 << dimensions[2]);
511 }
512 }
513
514 //----------------------------------------------------------------------------
ComputeTransformedBounds(int bounds[6])515 void vtkVolume16Reader::ComputeTransformedBounds (int bounds[6])
516 {
517 double transformedBounds[4];
518
519 if (!this->Transform)
520 {
521 bounds[0] = 0;
522 bounds[1] = this->DataDimensions[0] - 1;
523 bounds[2] = 0;
524 bounds[3] = this->DataDimensions[1] - 1;
525 bounds[4] = 0;
526 bounds[5] = this->ImageRange[1] - this->ImageRange[0];
527 }
528 else
529 {
530 transformedBounds[0] = 0;
531 transformedBounds[1] = 0;
532 transformedBounds[2] = 0;
533 transformedBounds[3] = 1.0;
534 this->Transform->MultiplyPoint (transformedBounds, transformedBounds);
535 bounds[0] = (int) transformedBounds[0];
536 bounds[2] = (int) transformedBounds[1];
537 bounds[4] = (int) transformedBounds[2];
538 transformedBounds[0] = this->DataDimensions[0] - 1;
539 transformedBounds[1] = this->DataDimensions[1] - 1;
540 transformedBounds[2] = this->ImageRange[1] - this->ImageRange[0];
541 transformedBounds[3] = 1.0;
542 this->Transform->MultiplyPoint (transformedBounds, transformedBounds);
543 bounds[1] = (int) transformedBounds[0];
544 bounds[3] = (int) transformedBounds[1];
545 bounds[5] = (int) transformedBounds[2];
546 // put bounds in correct order
547 int tmp;
548 for (int i = 0; i < 6; i += 2)
549 {
550 if (bounds[i + 1] < bounds[i])
551 {
552 tmp = bounds[i];
553 bounds[i] = bounds[i + 1];
554 bounds[i + 1] = tmp;
555 }
556 }
557 vtkDebugMacro(<< "Transformed bounds are: "
558 << bounds[0] << ", " << bounds[1] << ", "
559 << bounds[2] << ", " << bounds[3] << ", "
560 << bounds[4] << ", " << bounds[5]);
561 }
562 }
563
564 //----------------------------------------------------------------------------
AdjustSpacingAndOrigin(int dimensions[3],double spacing[3],double origin[3])565 void vtkVolume16Reader::AdjustSpacingAndOrigin (int dimensions[3], double spacing[3],
566 double origin[3])
567 {
568 for (int i = 0; i < 3; i++)
569 {
570 if (spacing[i] < 0)
571 {
572 origin[i] = origin[i] + spacing[i] * dimensions[i];
573 spacing[i] = -spacing[i];
574 }
575 }
576 vtkDebugMacro("Adjusted Spacing " << spacing[0] << ", " << spacing[1] << ", " << spacing[2]);
577 vtkDebugMacro("Adjusted origin " << origin[0] << ", " << origin[1] << ", " << origin[2]);
578 }
579
580 //----------------------------------------------------------------------------
TransformSlice(unsigned short * slice,unsigned short * pixels,int k,int dimensions[3],int bounds[3])581 void vtkVolume16Reader::TransformSlice (unsigned short *slice, unsigned short *pixels, int k, int dimensions[3], int bounds[3])
582 {
583 int iSize = this->DataDimensions[0];
584 int jSize = this->DataDimensions[1];
585
586 if (!this->Transform)
587 {
588 memcpy (pixels + iSize * jSize * k, slice, iSize * jSize * sizeof (unsigned short));
589 }
590 else
591 {
592 double transformedIjk[4], ijk[4];
593 int i, j;
594 int xyz[3];
595 int index;
596 int xSize = dimensions[0];
597 int xySize = dimensions[0] * dimensions[1];
598
599 // now move slice into pixels
600
601 ijk[2] = k;
602 ijk[3] = 1.0;
603 for (j = 0; j < jSize; j++)
604 {
605 ijk[1] = j;
606 for (i = 0; i < iSize; i++, slice++)
607 {
608 ijk[0] = i;
609 this->Transform->MultiplyPoint (ijk, transformedIjk);
610 xyz[0] = (int) ((double)transformedIjk[0] - bounds[0]);
611 xyz[1] = (int) ((double)transformedIjk[1] - bounds[2]);
612 xyz[2] = (int) ((double)transformedIjk[2] - bounds[4]);
613 index = xyz[0] +
614 xyz[1] * xSize +
615 xyz[2] * xySize;
616 *(pixels + index) = *slice;
617 }
618 }
619 }
620 }
621
622 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)623 void vtkVolume16Reader::PrintSelf(ostream& os, vtkIndent indent)
624 {
625 this->Superclass::PrintSelf(os,indent);
626
627 os << indent << "HeaderSize: " << this->HeaderSize << "\n";
628 os << indent << "SwapBytes: " << this->SwapBytes << "\n";
629 os << indent << "Data Dimensions: (" << this->DataDimensions[0] << ", "
630 << this->DataDimensions[1] << ")\n";
631 os << indent << "Data Mask: " << this->DataMask << "\n";
632
633 if ( this->Transform )
634 {
635 os << indent << "Transform:\n";
636 this->Transform->PrintSelf(os,indent.GetNextIndent());
637 }
638 else
639 {
640 os << indent << "Transform: (None)\n";
641 }
642 }
643