[vtkusers] Using vtkImageReslice

[Audrius Stundzia]

Hi fellow vtker's,

I would appreciate any insight on using vtkImageReslice to extract 2D image 
slices from a 3D.

The catch (for me, anyways) is that to get the slices I want I first have to 
perform a rotation R(1), followed by a translation T and a second rotation 
R(2), ie., R(2) T R(1)

I've tested the rotation + translation tranform by extracting a surface and 
comparing it center and orientation to a reference grid and it checks out 
okay.

Step 1: Here I create and fill my initial vtkImageData object with origin at 
(0,0,0)

( nVoxels[X] = nVoxels[Y] and nVoxels[Z] < nVoxels[X] and

  voxelSize[X] = voxelSize[Y] and voxelSize[Z] > voxelSize[X] )

   imageDataVolume = vtkImageData::New();
   imageDataVolume->PrepareForNewData();
   imageDataVolume->SetDimensions( nVoxels[X]
                                 , nVoxels[Y]
                                 , nVoxels[Z] );
   imageDataVolume->SetOrigin( -firstCenter[X]
                             , -firstCenter[Y]
                             , -firstCenter[Z] );
   imageDataVolume->SetSpacing( voxelSize[X]
                              , voxelSize[Y]
                              , voxelSize[Z] );
   imageDataVolume->SetExtent( 0
                             , nVoxels[X] - 1
                             , 0
                             , nVoxels[Y] - 1
                             , 0
                             , nVoxels[Z] - 1 );

   imageDataVolume->SetScalarType( VTK_FLOAT);
   imageDataVolume->SetNumberOfScalarComponents( 1);
   imageDataVolume->AllocateScalars();
   imageDataVolume->GetPointData()->SetScalars( imageVolumeScalars);
   imageDataVolume->Update();

Step 2: First I reslice along the 1st rotation axes

   imageDataFirstRotationVolumeReslice = vtkImageReslice::New();
   imageDataFirstRotationVolumeReslice->SetInput( imageDataVolume);
   imageDataFirstRotationVolumeReslice
   ->SetResliceAxes( firstRoatationMatrix);
   imageDataFirstRotationVolumeReslice->WrapOff();
   imageDataFirstRotationVolumeReslice->MirrorOff();
        imageDataFirstRotationVolumeReslice
   ->SetInterpolationModeToLinear();
   imageDataFirstRotationVolumeReslice->InterpolateOn();
   imageDataFirstRotationVolumeReslice->OptimizationOn();
   imageDataFirstRotationVolumeReslice->ReleaseDataFlagOn();
   imageDataFirstRotationVolumeReslice->SetBackgroundColor( 0.0
                                                          , 0.0
                                                          , 0.0
                                                          , 0.0 );

   imageDataFirstRotationVolumeReslice
   ->OutputAlwaysCenteredOnInputOff();

   imageDataFirstRotationVolumeReslice
   ->SetOutputSpacing( voxelSize[X]
                     , voxelSize[Y]
                     , voxelSize[Z] );
                                                                
imageDataFirstRotationVolumeReslice
   ->SetOutputOrigin( -firstCenter[X]
                    , -firstCenter[Y]
                    , -firstCenter[Z] );

   imageDataFirstRotationVolumeReslice
   ->SetOutputExtent( 0
                    , ( nVoxels[X] - 1)
                    , 0
                    , ( nVoxels[Y] - 1)
                    , 0
                    , ( nVoxels[Z] - 1) );

   imageDataFirstRotationVolumeReslice->Update();

Question: as the SetOutputOrigin arguements are specified as (float) am I 
correct in assuming that one uses physical coordinate (x,y,z) rather than 
voxel coordinates (i,j,k)?

Step 3: Now I reslice along the 2nd rotation axes, taking the translation T 
into account by changing the output origin function call SetOutputOrigin.

Questions: Do I need to first update the origin of the input for the second 
reslice? Ie.,

imageDataFirstRotationVolumeReslice->GetOutput()
->SetOrigin( -secondCenter[X]
           , -secondCenter[Y]
           , -secondCenter[Z] );

to take the translation T into proper account?

   imageDataSecondRotationVolumeReslice = vtkImageReslice::New();
   imageDataSecondRotationVolumeReslice
   ->SetInput( imageDataFirstRotationVolumeReslice->GetOutput() );
   imageDataSecondRotationVolumeReslice
   ->SetResliceAxes( secondRotationMatrix);

   imageDataSecondRotationVolumeReslice
   -> same set of function calls as in the first reslice;

   imageDataSecondRotationVolumeReslice
   ->OutputAlwaysCenteredOnInputOff();

   imageDataSecondRotationVolumeReslice
   ->SetOutputSpacing( voxelSize[X]
                     , voxelSize[Y]
                     , voxelSize[Z] );

   imageDataSecondRotationVolumeReslice
   ->SetOutputOrigin( -secondCenter[X]
                    , -secondCenter[Y]
                    , -secondCenter[Z] );

   imageDataSecondRotationVolumeReslice
   ->SetOutputExtent( 0
                    , ( nVoxels[X] - 1)
                    , 0
                    , ( nVoxels[Y] - 1)
                    , 0
                    , ( nVoxels[Z] - 1) );

   imageDataSecondRotationVolumeReslice->Update();

Step 4: At this stage my reslicing should be done and I can go ahead and 
extract my 2D image slices along the axes defined by an identity matrix. The 
relevant parts

   imageXYPlaneReslice[iPlane]
   ->SetInput( imageDataSecondRotationVolumeReslice->GetOutput() );
   imageXYPlaneReslice[iPlane]->SetResliceAxes( identityMatrix);

   int zIndex
   = (int) ( ( secondCentre[Z] + GetPlaneZ()[iPlane] )
             / voxelSize[Z] );

   imageXYPlaneReslice[iPlane]->SetOutputExtent( 0
                                               , ( nVoxels[X] - 1)
                                               , 0
                                               , ( nVoxels[Y] - 1)
                                               , zIndex
                                               , zIndex );

Three problems:

1. Step 4 is slooooow. I assume this is what is meant by the  comment in the 
online manual, "Warning: This filter is very inefficient if the output X 
dimension is 1."

However, my extract plane voxel thickness needs to be 1 as I then need to 
texture map the output.

Should I use vtkExtractVOI instead? Any yabuts?

2. The translation was not done properly. It's off in all three axes.

3. The corresponding XZ and YZ plane images look like they've been    left 
out in the rain and have had their colours run along the Z axis.

Any insight would be most appreciated. Thanks.

Audrius




_________________________________________________________________________
Get Your Private, Free E-mail from MSN Hotmail at http://www.hotmail.com.

Share information about yourself, create your own public profile at 
http://profiles.msn.com.