[vtkusers] Correct coordinate transformations for DICOM data set.

[David Gobbi]

About the code that I included in my previous email, I was
sure that I would forget something, and I did!

5) The "Origin" of the vtkImageData that you get from the
DICOM reader _must_ be (0,0,0) or else the matrix will not
work.  The matrix already has the position information from
the DICOM header.  The VTK "Origin" is not the same as the
DICOM "Position", so you cannot get "Image Position (Patient)"
from the DICOM header and then use that as your Origin!

Likewise, if you ever import oriented ITK images into VTK,
the ITK Origin and the VTK Origin are not the same thing!
The ITK Origin is defined with respect to the oriented image
(I.e. with respect to a rotated coordinate system), while the
VTK Origin is defined with respect to an un-oriented image
(i.e. with respect to a coordinate system where X, Y, and Z
are aligned with the column, row, and slice directions).

I hope this clarifies things, rather than just confusing them
further!

 - David



On Wed, Aug 20, 2014 at 10:58 AM, David Gobbi <david.gobbi at gmail.com> wrote:
> 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