FixedPointVolumeRayCastMapperCT.cxx 15.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    FixedPointVolumeRayCastMapperCT.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
// VTK includes
#include "vtkBoxWidget.h"
#include "vtkCamera.h"
#include "vtkCommand.h"
#include "vtkColorTransferFunction.h"
#include "vtkDICOMImageReader.h"
#include "vtkImageData.h"
#include "vtkImageResample.h"
#include "vtkMetaImageReader.h"
#include "vtkPiecewiseFunction.h"
#include "vtkPlanes.h"
#include "vtkProperty.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkVolume.h"
#include "vtkVolumeProperty.h"
#include "vtkXMLImageDataReader.h"
#include "vtkFixedPointVolumeRayCastMapper.h"

#define VTI_FILETYPE 1
#define MHA_FILETYPE 2

void PrintUsage()
{
  cout << "Usage: " << endl;
  cout << endl;
  cout << "  FixedPointVolumeRayCastMapperCT <options>" << endl;
  cout << endl;
  cout << "where options may include: " << endl;
  cout << endl;
  cout << "  -DICOM <directory>" << endl;
  cout << "  -VTI <filename>" << endl;
  cout << "  -MHA <filename>" << endl;
  cout << "  -DependentComponents" << endl;
  cout << "  -Clip" << endl;
  cout << "  -MIP <window> <level>" << endl;
  cout << "  -CompositeRamp <window> <level>" << endl;
  cout << "  -CompositeShadeRamp <window> <level>" << endl;
  cout << "  -CT_Skin" << endl;
  cout << "  -CT_Bone" << endl;
  cout << "  -CT_Muscle" << endl;
  cout << "  -FrameRate <rate>" << endl;
  cout << "  -DataReduction <factor>" << endl;
  cout << endl;
  cout << "You must use either the -DICOM option to specify the directory where" << endl;
  cout << "the data is located or the -VTI or -MHA option to specify the path of a .vti file." << endl;
  cout << endl;
  cout << "By default, the program assumes that the file has independent components," << endl;
  cout << "use -DependentComponents to specify that the file has dependent components." << endl;
  cout << endl;
  cout << "Use the -Clip option to display a cube widget for clipping the volume." << endl;
  cout << "Use the -FrameRate option with a desired frame rate (in frames per second)" << endl;
  cout << "which will control the interactive rendering rate." << endl;
  cout << "Use the -DataReduction option with a reduction factor (greater than zero and" << endl;
  cout << "less than one) to reduce the data before rendering." << endl;
  cout << "Use one of the remaining options to specify the blend function" << endl;
  cout << "and transfer functions. The -MIP option utilizes a maximum intensity" << endl;
  cout << "projection method, while the others utilize compositing. The" << endl;
  cout << "-CompositeRamp option is unshaded compositing, while the other" << endl;
  cout << "compositing options employ shading." << endl;
  cout << endl;
  cout << "Note: MIP, CompositeRamp, CompositeShadeRamp, CT_Skin, CT_Bone," << endl;
  cout << "and CT_Muscle are appropriate for DICOM data. MIP, CompositeRamp," << endl;
  cout << "and RGB_Composite are appropriate for RGB data." << endl;
  cout << endl;
  cout << "Example: FixedPointVolumeRayCastMapperCT -DICOM CTNeck -MIP 4096 1024" << endl;
  cout << endl;
}

int main(int argc, char *argv[])
{
  // Parse the parameters

  int count = 1;
  char *dirname = NULL;
  double opacityWindow = 4096;
  double opacityLevel = 2048;
  int blendType = 0;
  int clip = 0;
  double reductionFactor = 1.0;
  double frameRate = 10.0;
  char *fileName=0;
  int fileType=0;

  bool independentComponents=true;

  while ( count < argc )
103
  {
104
    if ( !strcmp( argv[count], "?" ) )
105
    {
106
107
      PrintUsage();
      exit(EXIT_SUCCESS);
108
    }
109
    else if ( !strcmp( argv[count], "-DICOM" ) )
110
    {
111
112
113
      dirname = new char[strlen(argv[count+1])+1];
      sprintf( dirname, "%s", argv[count+1] );
      count += 2;
114
    }
115
    else if ( !strcmp( argv[count], "-VTI" ) )
116
    {
117
118
119
120
      fileName = new char[strlen(argv[count+1])+1];
      fileType = VTI_FILETYPE;
      sprintf( fileName, "%s", argv[count+1] );
      count += 2;
121
    }
122
    else if ( !strcmp( argv[count], "-MHA" ) )
123
    {
124
125
126
127
      fileName = new char[strlen(argv[count+1])+1];
      fileType = MHA_FILETYPE;
      sprintf( fileName, "%s", argv[count+1] );
      count += 2;
128
    }
129
    else if ( !strcmp( argv[count], "-Clip") )
130
    {
131
132
      clip = 1;
      count++;
133
    }
134
    else if ( !strcmp( argv[count], "-MIP" ) )
135
    {
136
137
138
139
      opacityWindow = atof( argv[count+1] );
      opacityLevel  = atof( argv[count+2] );
      blendType = 0;
      count += 3;
140
    }
141
    else if ( !strcmp( argv[count], "-CompositeRamp" ) )
142
    {
143
144
145
146
      opacityWindow = atof( argv[count+1] );
      opacityLevel  = atof( argv[count+2] );
      blendType = 1;
      count += 3;
147
    }
148
    else if ( !strcmp( argv[count], "-CompositeShadeRamp" ) )
149
    {
150
151
152
153
      opacityWindow = atof( argv[count+1] );
      opacityLevel  = atof( argv[count+2] );
      blendType = 2;
      count += 3;
154
    }
155
    else if ( !strcmp( argv[count], "-CT_Skin" ) )
156
    {
157
158
      blendType = 3;
      count += 1;
159
    }
160
    else if ( !strcmp( argv[count], "-CT_Bone" ) )
161
    {
162
163
      blendType = 4;
      count += 1;
164
    }
165
    else if ( !strcmp( argv[count], "-CT_Muscle" ) )
166
    {
167
168
      blendType = 5;
      count += 1;
169
    }
170
    else if ( !strcmp( argv[count], "-RGB_Composite" ) )
171
    {
172
173
      blendType = 6;
      count += 1;
174
    }
175
    else if ( !strcmp( argv[count], "-FrameRate") )
176
    {
177
178
      frameRate = atof( argv[count+1] );
      if ( frameRate < 0.01 || frameRate > 60.0 )
179
      {
180
181
182
183
        cout << "Invalid frame rate - use a number between 0.01 and 60.0" << endl;
        cout << "Using default frame rate of 10 frames per second." << endl;
        frameRate = 10.0;
      }
184
185
      count += 2;
    }
186
    else if ( !strcmp( argv[count], "-ReductionFactor") )
187
    {
188
189
      reductionFactor = atof( argv[count+1] );
      if ( reductionFactor <= 0.0 || reductionFactor >= 1.0 )
190
      {
191
192
193
194
        cout << "Invalid reduction factor - use a number between 0 and 1 (exclusive)" << endl;
        cout << "Using the default of no reduction." << endl;
        reductionFactor = 1.0;
      }
195
196
      count += 2;
    }
197
     else if ( !strcmp( argv[count], "-DependentComponents") )
198
     {
199
200
      independentComponents=false;
      count += 1;
201
     }
202
    else
203
    {
204
205
206
207
208
      cout << "Unrecognized option: " << argv[count] << endl;
      cout << endl;
      PrintUsage();
      exit(EXIT_FAILURE);
    }
209
  }
210
211

  if ( !dirname && !fileName)
212
  {
213
214
215
216
    cout << "Error: you must specify a directory of DICOM data or a .vti file or a .mha!" << endl;
    cout << endl;
    PrintUsage();
    exit(EXIT_FAILURE);
217
  }
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238

  // Create the renderer, render window and interactor
  vtkRenderer *renderer = vtkRenderer::New();
  vtkRenderWindow *renWin = vtkRenderWindow::New();
  renWin->AddRenderer(renderer);

  // Connect it all. Note that funny arithematic on the
  // SetDesiredUpdateRate - the vtkRenderWindow divides it
  // allocated time across all renderers, and the renderer
  // divides it time across all props. If clip is
  // true then there are two props
  vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
  iren->SetRenderWindow(renWin);
  iren->SetDesiredUpdateRate(frameRate / (1+clip) );

  iren->GetInteractorStyle()->SetDefaultRenderer(renderer);

  // Read the data
  vtkAlgorithm *reader=0;
  vtkImageData *input=0;
  if(dirname)
239
  {
240
241
242
243
244
    vtkDICOMImageReader *dicomReader = vtkDICOMImageReader::New();
    dicomReader->SetDirectoryName(dirname);
    dicomReader->Update();
    input=dicomReader->GetOutput();
    reader=dicomReader;
245
  }
246
  else if ( fileType == VTI_FILETYPE )
247
  {
248
249
250
251
252
    vtkXMLImageDataReader *xmlReader = vtkXMLImageDataReader::New();
    xmlReader->SetFileName(fileName);
    xmlReader->Update();
    input=xmlReader->GetOutput();
    reader=xmlReader;
253
  }
254
  else if ( fileType == MHA_FILETYPE )
255
  {
256
257
258
259
260
    vtkMetaImageReader *metaReader = vtkMetaImageReader::New();
    metaReader->SetFileName(fileName);
    metaReader->Update();
    input=metaReader->GetOutput();
    reader=metaReader;
261
  }
262
  else
263
  {
264
265
    cout << "Error! Not VTI or MHA!" << endl;
    exit(EXIT_FAILURE);
266
  }
267
268
269
270
271
272
273

  // Verify that we actually have a volume
  int dim[3];
  input->GetDimensions(dim);
  if ( dim[0] < 2 ||
       dim[1] < 2 ||
       dim[2] < 2 )
274
  {
275
276
    cout << "Error loading data!" << endl;
    exit(EXIT_FAILURE);
277
  }
278
279
280

  vtkImageResample *resample = vtkImageResample::New();
  if ( reductionFactor < 1.0 )
281
  {
282
283
284
285
    resample->SetInputConnection( reader->GetOutputPort() );
    resample->SetAxisMagnificationFactor(0, reductionFactor);
    resample->SetAxisMagnificationFactor(1, reductionFactor);
    resample->SetAxisMagnificationFactor(2, reductionFactor);
286
  }
287
288
289
290
291
292

  // Create our volume and mapper
  vtkVolume *volume = vtkVolume::New();
  vtkFixedPointVolumeRayCastMapper *mapper = vtkFixedPointVolumeRayCastMapper::New();

  if ( reductionFactor < 1.0 )
293
  {
294
    mapper->SetInputConnection( resample->GetOutputPort() );
295
  }
296
  else
297
  {
298
    mapper->SetInputConnection( reader->GetOutputPort() );
299
  }
300
301
302
303
304


  // Set the sample distance on the ray to be 1/2 the average spacing
  double spacing[3];
  if ( reductionFactor < 1.0 )
305
  {
306
    resample->GetOutput()->GetSpacing(spacing);
307
  }
308
  else
309
  {
310
    input->GetSpacing(spacing);
311
  }
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334

//  mapper->SetSampleDistance( (spacing[0]+spacing[1]+spacing[2])/6.0 );
//  mapper->SetMaximumImageSampleDistance(10.0);


  // Create our transfer function
  vtkColorTransferFunction *colorFun = vtkColorTransferFunction::New();
  vtkPiecewiseFunction *opacityFun = vtkPiecewiseFunction::New();

  // Create the property and attach the transfer functions
  vtkVolumeProperty *property = vtkVolumeProperty::New();
  property->SetIndependentComponents(independentComponents);
  property->SetColor( colorFun );
  property->SetScalarOpacity( opacityFun );
  property->SetInterpolationTypeToLinear();

  // connect up the volume to the property and the mapper
  volume->SetProperty( property );
  volume->SetMapper( mapper );

  // Depending on the blend type selected as a command line option,
  // adjust the transfer function
  switch ( blendType )
335
  {
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
    // MIP
    // Create an opacity ramp from the window and level values.
    // Color is white. Blending is MIP.
    case 0:
      colorFun->AddRGBSegment(0.0, 1.0, 1.0, 1.0, 255.0, 1.0, 1.0, 1.0 );
      opacityFun->AddSegment( opacityLevel - 0.5*opacityWindow, 0.0,
                              opacityLevel + 0.5*opacityWindow, 1.0 );
      mapper->SetBlendModeToMaximumIntensity();
      break;

    // CompositeRamp
    // Create a ramp from the window and level values. Use compositing
    // without shading. Color is a ramp from black to white.
    case 1:
      colorFun->AddRGBSegment(opacityLevel - 0.5*opacityWindow, 0.0, 0.0, 0.0,
                              opacityLevel + 0.5*opacityWindow, 1.0, 1.0, 1.0 );
      opacityFun->AddSegment( opacityLevel - 0.5*opacityWindow, 0.0,
                              opacityLevel + 0.5*opacityWindow, 1.0 );
      mapper->SetBlendModeToComposite();
      property->ShadeOff();
      break;

    // CompositeShadeRamp
    // Create a ramp from the window and level values. Use compositing
    // with shading. Color is white.
    case 2:
      colorFun->AddRGBSegment(0.0, 1.0, 1.0, 1.0, 255.0, 1.0, 1.0, 1.0 );
      opacityFun->AddSegment( opacityLevel - 0.5*opacityWindow, 0.0,
                              opacityLevel + 0.5*opacityWindow, 1.0 );
      mapper->SetBlendModeToComposite();
      property->ShadeOn();
      break;

    // CT_Skin
    // Use compositing and functions set to highlight skin in CT data
    // Not for use on RGB data
    case 3:
      colorFun->AddRGBPoint( -3024, 0, 0, 0, 0.5, 0.0 );
      colorFun->AddRGBPoint( -1000, .62, .36, .18, 0.5, 0.0 );
      colorFun->AddRGBPoint( -500, .88, .60, .29, 0.33, 0.45 );
      colorFun->AddRGBPoint( 3071, .83, .66, 1, 0.5, 0.0 );

      opacityFun->AddPoint(-3024, 0, 0.5, 0.0 );
      opacityFun->AddPoint(-1000, 0, 0.5, 0.0 );
      opacityFun->AddPoint(-500, 1.0, 0.33, 0.45 );
      opacityFun->AddPoint(3071, 1.0, 0.5, 0.0);

      mapper->SetBlendModeToComposite();
      property->ShadeOn();
      property->SetAmbient(0.1);
      property->SetDiffuse(0.9);
      property->SetSpecular(0.2);
      property->SetSpecularPower(10.0);
      property->SetScalarOpacityUnitDistance(0.8919);
      break;

    // CT_Bone
    // Use compositing and functions set to highlight bone in CT data
    // Not for use on RGB data
    case 4:
      colorFun->AddRGBPoint( -3024, 0, 0, 0, 0.5, 0.0 );
      colorFun->AddRGBPoint( -16, 0.73, 0.25, 0.30, 0.49, .61 );
      colorFun->AddRGBPoint( 641, .90, .82, .56, .5, 0.0 );
      colorFun->AddRGBPoint( 3071, 1, 1, 1, .5, 0.0 );

      opacityFun->AddPoint(-3024, 0, 0.5, 0.0 );
      opacityFun->AddPoint(-16, 0, .49, .61 );
      opacityFun->AddPoint(641, .72, .5, 0.0 );
      opacityFun->AddPoint(3071, .71, 0.5, 0.0);

      mapper->SetBlendModeToComposite();
      property->ShadeOn();
      property->SetAmbient(0.1);
      property->SetDiffuse(0.9);
      property->SetSpecular(0.2);
      property->SetSpecularPower(10.0);
      property->SetScalarOpacityUnitDistance(0.8919);
      break;

    // CT_Muscle
    // Use compositing and functions set to highlight muscle in CT data
    // Not for use on RGB data
    case 5:
      colorFun->AddRGBPoint( -3024, 0, 0, 0, 0.5, 0.0 );
      colorFun->AddRGBPoint( -155, .55, .25, .15, 0.5, .92 );
      colorFun->AddRGBPoint( 217, .88, .60, .29, 0.33, 0.45 );
      colorFun->AddRGBPoint( 420, 1, .94, .95, 0.5, 0.0 );
      colorFun->AddRGBPoint( 3071, .83, .66, 1, 0.5, 0.0 );

      opacityFun->AddPoint(-3024, 0, 0.5, 0.0 );
      opacityFun->AddPoint(-155, 0, 0.5, 0.92 );
      opacityFun->AddPoint(217, .68, 0.33, 0.45 );
      opacityFun->AddPoint(420,.83, 0.5, 0.0);
      opacityFun->AddPoint(3071, .80, 0.5, 0.0);

      mapper->SetBlendModeToComposite();
      property->ShadeOn();
      property->SetAmbient(0.1);
      property->SetDiffuse(0.9);
      property->SetSpecular(0.2);
      property->SetSpecularPower(10.0);
      property->SetScalarOpacityUnitDistance(0.8919);
      break;

    // RGB_Composite
    // Use compositing and functions set to highlight red/green/blue regions
    // in RGB data. Not for use on single component data
    case 6:
      opacityFun->AddPoint(0, 0.0);
      opacityFun->AddPoint(5.0, 0.0);
      opacityFun->AddPoint(30.0, 0.05);
      opacityFun->AddPoint(31.0, 0.0);
      opacityFun->AddPoint(90.0, 0.0);
      opacityFun->AddPoint(100.0, 0.3);
      opacityFun->AddPoint(110.0, 0.0);
      opacityFun->AddPoint(190.0, 0.0);
      opacityFun->AddPoint(200.0, 0.4);
      opacityFun->AddPoint(210.0, 0.0);
      opacityFun->AddPoint(245.0, 0.0);
      opacityFun->AddPoint(255.0, 0.5);

      mapper->SetBlendModeToComposite();
      property->ShadeOff();
      property->SetScalarOpacityUnitDistance(1.0);
      break;
    default:
       vtkGenericWarningMacro("Unknown blend type.");
       break;
464
  }
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493

  // Set the default window size
  renWin->SetSize(600,600);
  renWin->Render();

  // Add the volume to the scene
  renderer->AddVolume( volume );

  renderer->ResetCamera();

  // interact with data
  renWin->Render();

  iren->Start();

  opacityFun->Delete();
  colorFun->Delete();
  property->Delete();

  volume->Delete();
  mapper->Delete();
  reader->Delete();
  resample->Delete();
  renderer->Delete();
  renWin->Delete();
  iren->Delete();

  return 0;
}