[vtkusers] question on vtkImageData orientation support

Mark Roden mmroden at gmail.com
Fri Nov 5 11:10:07 EDT 2010


The inverted image is actually a pretty tricky problem in vtk when dealing
with DICOM images.  DICOM is very specific about how the orientation vectors
are to be used, and VTK absolutely ignores this standard.  I'm not sure if
it's bugged behavior or not, though.

Right now, there's an option in VTK to flip the Y coordinate, but I don't
remember what it is and it turns out not to be helpful here, because you'll
also need to flip x and z, most likely.

DICOM specifies that each imaging plane has the upper left xy coordinate in
each 2d slice specified.  That coordinate is in tag 0020,0032 (Image
Position (Patient)).  DICOM then specifies that each subsequent pixel in the
x y plane is determined by the spacing from that coordinate (generally
positive).  So, if you have an x extent of -300 and a pixel spacing of 1 on
a 512x512 image, then your x coordinate will span from -300 to 212.  This is
true regardless of the direction cosines in the image.  The direction
cosines tell the user the orientation of the patient as regards to the
scanner, and should be used to determine display, but not the coordinates of
the data.

In VTK, however, that's not true.  If the x direction cosine is -1,0,0, then
your 512x512 image will be from -300 to -812.  That means that all the
coordinates in your RTStruct will be completely wrong, and your RTStruct
will get drawn off in space some where.  That offset problem happens because
the RT contour coordinates are exactly specified to the pixel coordinates I
described in the previous paragraph.

The solution is somewhat complicated, and hindered by the fact that VTK has
some functions in the API that don't work but are apparently there to
confuse people.  You need to avoid PreserveImageExtentOn when you flip, as
it does _absolutely nothing_.

The next problem is that vtk doesn't allow you to flip about any origin
other than 0,0,0.  (Remember, we have to flip the data itself, not the
cameras, because otherwise all the RT information will be misregistered to
the original image, regardless of where you're looking).

This is the code I use in Java to solve this problem:

                IPPSorter sorter = new IPPSorter();
                sorter.SetComputeZSpacing(true);
                sorter.SetZSpacingTolerance(0.000001);
                Boolean sorted = sorter.Sort(theSeriesFiles);
                //sorter.Sort(ft);

                //have to make a vtk string array
                FilenamesType sortedFT = sorter.GetFilenames();
                long theSize = sortedFT.size();
                vtkStringArray sa = new vtkStringArray();
                for (int j = 0; j < theSize; j++) {
                    String theString2 = sortedFT.get(j);
                    sa.InsertNextValue(theString2);
                }
                gdcmReader.SetFileNames(sa);

                gdcmReader.Update();

                //make sure that we're getting spacing and orientation info
                double[] spacing = gdcmReader.GetOutput().GetSpacing();
                spacing[2] = sorter.GetZSpacing();
                gdcmReader.GetOutput().SetSpacing(spacing);
                double[] origin = gdcmReader.GetOutput().GetOrigin();
                double[] bounds = gdcmReader.GetOutput().GetBounds();
                double[] nextOrigin = {0,0,0};//because vtk
                //will only flip about 0,0,0
                vtkImageChangeInformation imageSpacingChange = new
vtkImageChangeInformation();
                imageSpacingChange.SetInput(gdcmReader.GetOutput());
                imageSpacingChange.SetOutputSpacing(spacing[0], spacing[1],
spacing[2]);
                imageSpacingChange.SetOutputOrigin(nextOrigin );
                imageSpacingChange.Update();
                imageData = imageSpacingChange.GetOutput();
                //ok, here's the deal with axes
                //for now, assume all axes are -1 or 1 along the origin
                //and then set them all to 1.
                vtkMatrix4x4 directionCosines =
gdcmReader.GetDirectionCosines();
                //actually have to flip the data, not the camera
                //and the opposite way to itk.
                //but, to make the image display properly, have to flip the
camera
                //anyway? most likely, but depends if this is brain, pelvis,
etc
                boolean willFlip = false;
                for (int j = 0; j < 3; j++){
                    if (directionCosines.GetElement(j, j) < 0) {
                        flip[j] = true;
                        willFlip = true;
                    }
                    else {
                        flip[j] = false;
                    }
                    if (j == 1 ) {
                        flip[j] = !flip[j]; //y is inherently backwards in
vtk
                        willFlip = true;
                    }
                }


                //vtk's pipelining means I can't easily iterate over axes.
                //note the switch that doesn't do anything
                //(that would be PreserveImageExtentOn/Off)
                //require that you translate, flip, translate to get the
                //functionality of this older switch
                if (willFlip){

                    vtk.vtkImageFlip theFlip = new vtk.vtkImageFlip();

theFlip.SetInputConnection(imageSpacingChange.GetOutputPort());
                    theFlip.FlipAboutOriginOn();
                    //gotta find the first axis to flip...
                    boolean keepGoing = false;
                    if (flip[0]){
                        theFlip.SetFilteredAxis(0);
                        keepGoing = flip[1] && flip[2];
                    } else if (flip[1]){
                        theFlip.SetFilteredAxis(1);
                        keepGoing = flip[2];
                    } else {
                        theFlip.SetFilteredAxis(2);
                    }

                    theFlip.Update();
                    imageData = theFlip.GetOutput();
                    if (keepGoing){
                        vtk.vtkImageFlip theFlip2 = new vtk.vtkImageFlip();

theFlip2.SetInputConnection(theFlip.GetOutputPort());
                        theFlip2.FlipAboutOriginOn();
                        //gotta find the first axis to flip...
                        boolean keepGoing2 = false;
                        if (flip[1]){
                            theFlip2.SetFilteredAxis(1);
                            keepGoing = flip[2];
                        } else {
                            theFlip2.SetFilteredAxis(2);
                        }

                        theFlip2.Update();
                        imageData = theFlip2.GetOutput();

                        if (keepGoing2){
                            vtk.vtkImageFlip theFlip3 = new
vtk.vtkImageFlip();

theFlip3.SetInputConnection(theFlip2.GetOutputPort());
                            theFlip3.FlipAboutOriginOn();
                            //gotta find the first axis to flip...
                            theFlip3.SetFilteredAxis(2);

                            theFlip3.Update();
                            imageData = theFlip3.GetOutput();

                        }
                    }
                }
                //have to do this because the flipping above only flips
about 0,0,0
                vtkImageChangeInformation imageSpacingChange2 = new
vtkImageChangeInformation();
                imageSpacingChange2.SetInput(imageData);
                imageSpacingChange2.SetOutputSpacing(spacing[0], spacing[1],
spacing[2]);
                double[] extent = new double[3];
                for (int j = 0; j < 3; j++){
                    extent[j] = bounds[j*2+1]-bounds[j*2];
                }
                imageSpacingChange2.SetOutputOrigin(
                        origin[0] + (flip[0] ? -1 : 0) *extent[0],
                        origin[1] + (flip[1] ? -1 : 0) *extent[1],
                        origin[2] + (flip[2] ? -1 : 0) *extent[2]);
                imageSpacingChange2.Update();
                imageData = imageSpacingChange2.GetOutput();
                imageMedicalProperties =
gdcmReader.GetMedicalImageProperties();

The above code will NOT work if you have direction cosines that are not
unitary (ie, if you have direction cosines that are something like 0.93,
-0.4, etc etc).  That's because you'll have to add some rotational element
in as well.  This rotational problem could become a real issue when dealing
with MR data and associated RT structs, because MR does not require that the
directional cosines be unitary.

I hope that helps,
Mark

-------------
You can use vtkTransform to rotate the RT Dose matric to your desired
orientation.

Jothy

On Fri, Nov 5, 2010 at 2:09 AM, Xiaofeng Z <xf10036 at hotmail.com> wrote:

>  Hi, vtkusers,
>
> I have this algorithm that computes dose-volume histogram out of RT
> structureset and RT Dose.  Right now the algorithm works only when RT Dose
> has a 1,0,0 and 0,1,0 orientation.  I need to use the algorithm for other
RT
> Dose orientations.  I have too options, either to modify the algorithm, or
> to convert RT Dose orientation to 1,0,0 and 0,1,0.  I prefer to former,
but
> only if vtkImageData support multiple orientations.
>
> So my question is, can vtkImageData store pixel data in different
> orientation, and if so, how the orientation information is set and get?
>
> Thanks!
>
> Xiaofeng
>
>
>
>
> ------------------------------
> From: xf10036 at hotmail.com
> To: vtkusers at vtk.org
> Date: Thu, 4 Nov 2010 21:02:08 -0400
> Subject: [vtkusers] question on vtkImageData orientation support
>
>
> Does vtkImageData:SetExt(int ext[]) support inverse orientation, such that
> ext[0]<ext[1], ext[2]<ext[3], etc?
>
> Thanks!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20101105/fa8f2eae/attachment.htm>


More information about the vtkusers mailing list