[vtkusers] Correct coordinate transformations for DICOM data set.

[David Gobbi]

Hi V,

I found the vtkDICOMImageReader to be so lacking in functionality
that I wrote my own: http://dgobbi.github.io/vtk-dicom/

And, of course, many other people use ITK to read their DICOMs
(though you have to be very careful or you will lose coord system
fidelity when importing the image from ITK to VTK).  Finally, there
is also the vtkGDCMImageReader.


But if you are using the vtkDICOMImageReader, here are a few tips:

1) The reader ignores FileLowerLeft, so turning it on does nothing

2) The reader always applies a vertical flip to the images (nearly all
VTK image readers flip images when they read them).

3) The reader uses the "Image Position (Patient)" to sort the slices,
but it reverses the order.  It has to do this, because it has already
done a "y" flip, and two flips in total are necessary to avoid having
the volume become a mirror image of itself.

4) When you call reader->GetImagePositionPatient(), it returns the
information from the last slice that it read.  Since it reads the
slices in reverse order (see #3), this actually means that it returns
the information from the first slice.


So, putting everything together, here is the code that I used to read
images with the vtkDICOMImageReader (before I wrote my own reader):

  vtkSmartPointer<vtkDICOMImageReader> reader =
    vtkSmartPointer<vtkDICOMImageReader>::New();
  reader->SetDirectoryName(directoryName);

  // flip the image in Y and Z directions
  vtkSmartPointer<vtkImageReslice> flip =
    vtkSmartPointer<vtkImageReslice>::New();
  flip->SetInputConnection(reader->GetOutputPort());
  flip->SetResliceAxesDirectionCosines(
    1,0,0, 0,-1,0, 0,0,-1);

  flip->Update();
  vtkSmartPointer<vtkImageData> image = flip->GetOutput();

  // generate the matrix
  float *position = reader->GetImagePositionPatient();
  float *orientation = reader->GetImageOrientationPatient();
  float *xdir = &orientation[0];
  float *ydir = &orientation[3];
  float zdir[3];
  vtkMath::Cross(xdir, ydir, zdir);

  vtkSmartPointer<vtkMatrix4x4> matrix =
    vtkSmartPointer<vtkMatrix4x4>::New();
  for (int i = 0; i < 3; i++)
    {
    matrix->Element[i][0] = xdir[i];
    matrix->Element[i][1] = ydir[i];
    matrix->Element[i][2] = zdir[i];
    matrix->Element[i][3] = position[i];
    }


Then you can attach the matrix to the vtkVolume with
SetUserMatrix(matrix).  If you have points that are already in
the DICOM patient coordinate system, then you don't need to
apply any transformations to those points.  I'm not going to give
detailed info about how to display points in VTK (you can probably
find wiki examples), but my advice is to use vtkGlyph3D with
vtkSphereSource.  Just search for "VTK Glyph3D" and you'll
find a suitable example.

 - David


On Wed, Aug 20, 2014 at 9:55 AM, vipulraikar <vpai at g.clemson.edu> wrote:
> Coordinate transformations for DICOM data set.
>
> Hello everybody,
>
> Before I begin to explain my roadblock\problem, let me briefly put forward
> what I am trying to achieve. I am using VTK to handle the visualization in
> my work with Image-Guidance (IGS). At the current time I am trying to
> register a tracked tool (already done) to a data-set and render it correctly
> in 3D space on screen.
>
> To begin with I have extensively gone through the posts here on
> understanding how to visualize the 3D space in DICOM coordinate frame
> (Patient\Scanner\Reference).  To begin with, I am generating a vtkMatrix4X4
> with direction cosines and the origin from the dicom header of the first
> image in the dataset (LowerLeft is ON). I have then applied this as a user
> matrix to my vtkVolume, which according to my understanding from all the
> wonderful explanations here (Especially Dr Gobbi), will apply the DICOM
> coordinate frame to the rendering (for the lack of better way to word this).
>
> My questions are two fold. First: How to I now render a point, which is
> registered to this coordinate frame, correctly (in the form of a vtkActor)?
> Second: How to I move this object\actor in real-time?
>
> I have tried and failed to correctly do so by applying the same user
> transform to the actor. I am still quite new to this. I would very much
> appreciate any help I can get to move forward in the right direction.
>
> best regards,
>
> V