# reslicesphere.cxx

    /*=========================================================================

      Program: GDCM (Grassroots DICOM). A DICOM library
      Module:  $URL$

      Copyright (c) 2006-2010 Mathieu Malaterre
      All rights reserved.
      See Copyright.txt or http://gdcm.sourceforge.net/Copyright.html 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.

    =========================================================================*/
    //
    // Load a DICOM series.
    // Position a sphere within the volume.
    // Allow the user to change between Axial, Sagittal, Coronal, and
    //     Oblique view of the images and move through the slices.
    // The display should show the resliced image and the cross section
    // of the sphere intersecting that plane.
    //

    /*
    from  Scott Johnson /Scott Johnson neuwave com/
    to  VTK /vtkusers vtk.org/
    date  Tue, May 11, 2010 at 7:01 PM
    */
    #include <strstream>
    #include <string>

    #include <vtkDICOMImageReader.h>
    #include <vtkStringArray.h>
    #include <vtkDirectory.h>
    #include <vtkImageThreshold.h>
    #include <vtkImageShiftScale.h>
    #include <vtkImageReslice.h>
    #include <vtkRenderWindowInteractor.h>
    #include <vtkImageViewer2.h>
    #include <vtkSphereSource.h>
    #include <vtkPolyDataMapper.h>
    #include <vtkPlane.h>
    #include <vtkCutter.h>
    #include <vtkActor.h>
    #include <vtkCommand.h>
    #include <vtkSmartPointer.h>
    #include <vtkMatrix4x4.h>
    #include <vtkInteractorObserver.h>
    #include <vtkProperty.h>
    #include <vtkRenderer.h>
    #include <vtkImageData.h>
    #include <vtkImageActor.h>
    #include "vtkTransformPolyDataFilter.h"
    #include <vtkCamera.h>
    #include <vtkMath.h>
    #include <vtkTransform.h>
    #include <vtkTextActor.h>
    #include <vtkActor2D.h>
    #include <vtkPolyDataMapper2D.h>
    #include <vtkProperty2D.h>
    #include <vtkGDCMImageReader.h>
    #include <vtkImageChangeInformation.h>

    #include "gdcmDirectory.h"
    #include "gdcmTesting.h"
    #include "gdcmIPPSorter.h"

    // Change to match the path to find Raw_0.vti or provide
    // the parameter when starting ResliceSphere.

    const double sphereCenter[3]={74, 219, 70};

    // Angles (0, 0, 0)
    const double AxialMatrix[]    = { 1.0,  0.0,  0.0,  0.0,
                                      0.0,  1.0,  0.0,  0.0,
                                      0.0,  0.0,  1.0,  0.0,
                                      0.0,  0.0,  0.0,  1.0 };
    // Angles (0, 90, 0)
    const double SagittalMatrix[] = { 0.0,  0.0,  1.0,  0.0,
                                      0.0,  1.0,  0.0,  0.0,
                                     -1.0,  0.0,  0.0,  0.0,
                                      0.0,  0.0,  0.0,  1.0 };
    // Angles (-90, 0, 0)
    const double CoronalMatrix[]  = { 1.0,  0.0,  0.0,  0.0,
                                      0.0,  0.0,  1.0,  0.0,
                                      0.0, -1.0,  0.0,  0.0,
                                      0.0,  0.0,  0.0,  1.0 };
    // Angles (0, 90, 31)
    const double ObliqueMatrix[] =  { 0.0, -0.515038, 0.857167, 0.0,
                                      0.0,  0.857167, 0.515038, 0.0,
                                     -1.0,  0.0,      0.0,      0.0,
                                      0.0,  0.0,      0.0,      1.0 };

    class ResliceRender;

    // Class to handle key press events.
    class KeyCallback : public vtkCommand
    {
    public:
        static KeyCallback* New()
        {
            return new KeyCallback();
        }

        void Execute(vtkObject* caller, unsigned long eventId, void *calldata);
        void SetCallbackData(ResliceRender* reslice);

    protected:
        ResliceRender* _reslice;
    };

    class ResliceRender
    {
    public:
        typedef enum _ORIENTATION
        {
            AXIAL = 0,
            SAGITTAL = 1,
            CORONAL = 2,
            OBLIQUE = 3
        } ORIENTATION;

        ResliceRender()
        {
            _orientation=AXIAL;
        }

        ~ResliceRender()
        {
            _transform->Delete();
            _reader->Delete();
            _reslice->Delete();
            _interactor->Delete();
            _imageViewer->Delete();

            _sphere->Delete();
            _sphereMapper->Delete();
            _sphereActor->Delete();

            _plane->Delete();
            _cutter->Delete();
            _polyTransform->Delete();
            _ROIMapper->Delete();
            _ROIActor->Delete();

            _annotation->Delete();
        }

        void CreatePipeline(const char* fileName)
        {
            vtkProperty2D* props;

            //_reader=vtkXMLImageDataReader::New();
            //_reader->SetFileName(fileName);
            //_reader->Update();

            //_reader=qzDICOMImageReader::New();
            _reader=vtkGDCMImageReader::New();

          //vtkDirectory *d = vtkDirectory::New();
          //d->Open(fileName);
          //d->Print( std::cout );
          gdcm::Directory d;
          d.Load(fileName);
          gdcm::Directory::FilenamesType const &files = d.GetFilenames();

      gdcm::IPPSorter s;
      s.SetComputeZSpacing( true );
      s.SetZSpacingTolerance( 1e-3 );
      bool b = s.Sort( files );
      if( !b )
        {
        std::cerr << "Failed to sort:" << fileName << std::endl;
        //return ;
        }
      //std::cout << "Sorting succeeded:" << std::endl;
      //s.Print( std::cout );

      //std::cout << "Found z-spacing:" << std::endl;
      //std::cout << s.GetZSpacing() << std::endl;
      double ippzspacing = s.GetZSpacing();

      const std::vector<std::string> & sorted = s.GetFilenames();
      vtkStringArray *vtkfiles = vtkStringArray::New();
      std::vector< std::string >::const_iterator it = sorted.begin();
      for( ; it != sorted.end(); ++it)
        {
        const std::string &f = *it;
        vtkfiles->InsertNextValue( f.c_str() );
        }

            //_reader->SetDirectoryName(fileName);
            //_reader->SetFileNames( d->GetFiles() );
            _reader->SetFileNames( vtkfiles );
            _reader->Update();

      const vtkFloatingPointType *spacing = _reader->GetOutput()->GetSpacing();

      vtkImageChangeInformation *v16 = vtkImageChangeInformation::New();
      v16->SetInput( _reader->GetOutput() );
      v16->SetOutputSpacing( spacing[0], spacing[1], ippzspacing );
      v16->Update();

            _threshold=vtkImageThreshold::New();
            _threshold->ThresholdByUpper(-3024.0);
            _threshold->ReplaceOutOn();
            _threshold->SetOutValue(0.0);
            _threshold->SetInputConnection(v16->GetOutputPort());

            _shift=vtkImageShiftScale::New();
            _shift->SetShift(0);
            _shift->SetScale(1);
            _shift->SetInputConnection(_threshold->GetOutputPort());

            // Initialize the reslice with an axial orientation.
            vtkSmartPointer<vtkMatrix4x4> matrix =
                vtkSmartPointer<vtkMatrix4x4>::New();
            matrix->Identity();

            _transform = vtkTransform::New();
            _transform->SetMatrix(matrix);

            _reslice = vtkImageReslice::New();
            _reslice->SetOutputDimensionality(3);

            // PROBLEM:
            // The original intent was to connect the same transform
            // to the vtkImageReslice and vtkTransformPolyDataFilter,
            // but the resulting reslices appear different using the
            // vtkTransform as opposed to explicitly setting the
            // reslice axes via SetResliceAxes.  Also, if the vtkTransform
            // is connected and orientated other than axial, the extents
            // don't seem to update resulting in VTK believing the slice
            // is out of range.

            //_reslice->SetResliceTransform(_transform);
            _reslice->SetResliceAxes(matrix);
            //_reslice->SetInputConnection(_reader->GetOutputPort());
            _reslice->SetInputConnection(_shift->GetOutputPort());

            // Create the sphere target shape.
            _sphere=vtkSphereSource::New();
            _sphere->SetRadius(7.0);
            _sphere->SetThetaResolution(16);
            _sphere->SetPhiResolution(16);
            _sphere->SetCenter(sphereCenter[0], sphereCenter[1], sphereCenter[2]);

            _sphereMapper=vtkPolyDataMapper::New();
            _sphereMapper->SetInputConnection(_sphere->GetOutputPort());

            _sphereActor=vtkActor::New();
            _sphereActor->SetMapper(_sphereMapper);
            _sphereActor->PickableOff();
            _sphereActor->GetProperty()->SetColor(1.0, 0.0, 0.0);
            _sphereActor->GetProperty()->SetEdgeColor(1.0, 0.0, 0.0);
            _sphereActor->GetProperty()->SetDiffuseColor(1.0, 0.0, 0.0);
            _sphereActor->SetVisibility(true);

            // Create the cutting pipeline.
            // This plane will be positioned in the original image coordinate system.
            _plane = vtkPlane::New();
            _plane->SetNormal(0.0, 0.0, 1.0);

            _cutter = vtkCutter::New();
            _cutter->SetInputConnection(_sphere->GetOutputPort());
            _cutter->SetCutFunction(_plane);
            _cutter->GenerateCutScalarsOn();
            _cutter->SetValue(0, 0.5);

            // The transform attached to _polyTransform should move the cut
            // ROI into the resliced coordinate system, which should be the
            // same as the coordinate system of the resliced images.
            // PROBLEM:  It doesn't.
            _polyTransform = vtkTransformPolyDataFilter::New();
            _polyTransform->SetTransform(_transform);
            _polyTransform->SetInputConnection(_cutter->GetOutputPort());

            _ROIMapper = vtkPolyDataMapper2D::New();
            _ROIMapper->SetInputConnection(_polyTransform->GetOutputPort());

        vtkCoordinate* coordinate = vtkCoordinate::New();
        coordinate->SetCoordinateSystemToWorld();
        _ROIMapper->SetTransformCoordinate(coordinate);

            _ROIActor = vtkActor2D::New();
            _ROIActor->SetMapper(_ROIMapper);

            // Make sure the cut can be seen, especially the edges.
            props=_ROIActor->GetProperty();
            props->SetLineWidth(2);
            props->SetOpacity(1.0);
    //        props->EdgeVisibilityOn();
    //        props->SetDiffuse(0.8);
    //        props->SetSpecular(0.3);
    //        props->SetSpecularPower(20);
    //        props->SetRepresentationToSurface();
    //        props->SetDiffuseColor(1.0, 0.0, 0.0);
    //        props->SetEdgeColor(1.0, 0.0, 0.0);
            props->SetColor(1.0, 0.0, 0.0);

            _interactor = vtkRenderWindowInteractor::New();

            // Create the image viewer and add the actor with the cut ROI.
            _imageViewer = vtkImageViewer2::New();
            _imageViewer->SetupInteractor(_interactor);
            _imageViewer->SetSize(400, 400);
            _imageViewer->SetColorWindow(1024);
            _imageViewer->SetColorLevel(800);
            _imageViewer->SetInputConnection(_reslice->GetOutputPort());
            _imageViewer->GetImageActor()->SetOpacity(0.5);

            _annotation = vtkTextActor::New();
            _annotation->SetTextScaleModeToViewport();
            _imageViewer->GetRenderer()->AddActor(_annotation);

            // Add the cut shape actor to the renderer.
            _imageViewer->GetRenderer()->AddActor(_ROIActor);

            // Set up the key handler.
            vtkSmartPointer<KeyCallback> callback = vtkSmartPointer<KeyCallback>::New();
            callback->SetCallbackData(this);
            _interactor->AddObserver(vtkCommand::KeyPressEvent, callback);

            _interactor->Initialize();
        }

        void Start()
        {
            _interactor->Start();
        }

        void ResetOrientation()
        {
            vtkSmartPointer<vtkMatrix4x4> matrix =
                vtkSmartPointer<vtkMatrix4x4>::New();
            matrix->Identity();

            SetOrientation(matrix);
        }

        // Make sure the orientation of the vtkImageReslice and
        // vtkTransform are in sync.
        void SetOrientation(vtkMatrix4x4* matrix)
        {
            _reslice->SetResliceAxes(matrix);
            _reslice->Update();

        vtkMatrix4x4* inverse = vtkMatrix4x4::New();
        vtkMatrix4x4::Invert(matrix, inverse);

            _transform->SetMatrix(inverse);
            _transform->Update();
        }

        // Set the current slice of the current view.
        void SetSlice(int slice)
        {
            std::strstream posString;

            double    center[3];
            double    spacing[3];
            double    origin[3];
            double    point[4];
            double    newPoint[4];

            vtkImageData* imageData;
            int newSlice;

            // Try to make sure the extents of the reslice are updated.
            // PROBLEM:  It doesn't seem to work when changing the orientation.
            imageData=vtkImageData::SafeDownCast(_reslice->GetOutput());
            imageData->UpdateInformation();

            // Let vtkImageViewer2 handle the slice limits.
            _imageViewer->SetSlice(slice);
            newSlice=GetSlice();

            imageData->GetCenter(center);
            imageData->GetSpacing(spacing);
            imageData->GetOrigin(origin);

            // Compute the position of the center of the slice based on the
            // spacing of the slices.  The resliced axis will always
            // be the "Z" axis.
            point[0]=center[0];
            point[1]=center[1];
            point[2]=(newSlice * spacing[2]) + origin[2];
            point[3]=1.0;

            // Convert the coordinate from the reslice coordinate system to the
            // original image coordinate system.
            // PROBLEM:  Logically this seems like it should have been multiplied
            // by the inverse to translate from the resliced coordinate system to
            // the original coordinate system.  However, multiplying by the inverse
            // sticks the plane in the wrong place completely.  Using the original
            // matrix at least gets the Z coordinate right.
            vtkMatrix4x4* matrix=_reslice->GetResliceAxes();
            vtkSmartPointer<vtkMatrix4x4> inverse =
                vtkSmartPointer<vtkMatrix4x4>::New();
            vtkMatrix4x4::Invert(matrix, inverse);

            matrix->MultiplyPoint(point, newPoint);
            _plane->SetOrigin(newPoint[0], newPoint[1], newPoint[2]);

            // Annotate the image.
            posString << "Position: (" << newPoint[0] << ", " << newPoint[1]
                      << ", " << newPoint[2] << ")  Slice: " << newSlice;
            _annotation->SetInput(posString.str());

            _imageViewer->Render();
        }

        int GetSlice()
        {
            return _imageViewer->GetSlice();
        }

        // Set the orientation of the view.
        void SetOrientation(ResliceRender::ORIENTATION orientation)
        {
            vtkCamera* camera=_imageViewer->GetRenderer()->GetActiveCamera();

            double spacing[3];
            double origin[3];
            double point[4];
            double newPoint[4];
            double initialPosition;
            double xDirCosine[3];
            double yDirCosine[3];
            double zDirCosine[3];
            double normal[3];

            vtkImageData* imageData;

            vtkSmartPointer<vtkMatrix4x4> matrix =
                vtkSmartPointer<vtkMatrix4x4>::New();

            _orientation=orientation;

            // Reset ViewUp
            camera->SetViewUp(0.0, 1.0, 0.0);

            // Compute the cut plane position to the input coordinate system.
            imageData=vtkImageData::SafeDownCast(_reslice->GetInput());
            imageData->UpdateInformation();
            imageData->GetSpacing(spacing);
            imageData->GetOrigin(origin);

            point[0]=origin[0];
            point[1]=origin[1];
            point[2]=origin[2];
            point[3]=1.0;

            switch (_orientation)
            {
            case AXIAL:
                matrix->DeepCopy(AxialMatrix);
                initialPosition=sphereCenter[2];
                break;

            case CORONAL:
                matrix->DeepCopy(CoronalMatrix);
                initialPosition=sphereCenter[1];
                break;

            case SAGITTAL:
                matrix->DeepCopy(SagittalMatrix);
                initialPosition=sphereCenter[0];
                break;

            case OBLIQUE:
                matrix->DeepCopy(ObliqueMatrix);
                initialPosition=sphereCenter[2];
                break;
            }

            // Move the origin from the original image coordinate system to the
            // resliced image coordinate system.
            matrix->MultiplyPoint(point, newPoint);
            matrix->SetElement(0, 3, newPoint[0]);
            matrix->SetElement(1, 3, newPoint[1]);
            matrix->SetElement(2, 3, newPoint[2]);

            ResetOrientation();
            SetOrientation(matrix);

            // Compute the cutting plane normal and set it.
            // PROBLEM:  If the transformation is connected rather than
            // using SetResliceAxes, the Direction Cosines do not reflect
            // the orientation of the vtkImageReslice.
            _reslice->GetResliceAxesDirectionCosines(xDirCosine, yDirCosine,
                                                     zDirCosine);
            vtkMath::Cross(xDirCosine, yDirCosine, normal);
            _plane->SetNormal(normal);

            // Set the extents and spacing of the reslice to account for
            // all of the data.
            _reslice->SetOutputExtentToDefault();
            _reslice->SetOutputSpacing(spacing[0], spacing[0], spacing[0]);

            // Force the vtkImageViewer2 to update.
            // PROBLEM:  The whole extent does not seem to be set in time
            // for the first render.  This results in an error because the
            // slice is positioned outside the old bounds.
            _imageViewer->SetInput(NULL);
            _imageViewer->SetInputConnection(_reslice->GetOutputPort());

            _imageViewer->GetRenderer()->ResetCameraClippingRange();
            _imageViewer->GetRenderer()->ResetCamera();

            // Set the initial slice to be at the center of the sphere.
            // Divide by the spacing because this will be undone in SetSlice.
            SetSlice(initialPosition / spacing[0]);
        }

        vtkRenderWindowInteractor* GetInteractor()
        {
            return _interactor;
        }

    protected:
        ORIENTATION                 _orientation;

        //qzDICOMImageReader*        _reader;
        vtkGDCMImageReader*        _reader;
        vtkImageThreshold*          _threshold;
        vtkImageShiftScale*         _shift;
        vtkImageReslice*            _reslice;
        vtkRenderWindowInteractor*  _interactor;
        vtkImageViewer2*            _imageViewer;

        vtkSphereSource*            _sphere;
        vtkPolyDataMapper*          _sphereMapper;
        vtkActor*                   _sphereActor;

        vtkPlane*                   _plane;
        vtkCutter*                  _cutter;
        vtkTransform*               _transform;
        vtkTransformPolyDataFilter* _polyTransform;
        vtkPolyDataMapper2D*          _ROIMapper;
        vtkActor2D*                   _ROIActor;

        vtkTextActor*               _annotation;
    };

    // Catch KeyPress events.
    // Up Arrow   - increases the slice
    // Down Arrow - decreases the slice
    // 'A'        - sets the view to Axial
    // 'S'        - sets the view to Sagittal
    // 'C'        - sets the view to Coronal
    // 'O'        - set the view to Oblique

    void KeyCallback::Execute(vtkObject* caller, unsigned long eventId, void *calldata)
    {
        std::string sym=_reslice->GetInteractor()->GetKeySym();

        if (!sym.compare("Up"))
        {
            _reslice->SetSlice(_reslice->GetSlice() + 1);
        }
        else if (!sym.compare("Down"))
        {
            _reslice->SetSlice(_reslice->GetSlice() - 1);
        }
        else if ((!sym.compare("A")) || (!sym.compare("a")))
        {
            _reslice->SetOrientation(ResliceRender::AXIAL);
        }
        else if ((!sym.compare("C")) || (!sym.compare("c")))
        {
            _reslice->SetOrientation(ResliceRender::CORONAL);
        }
        else if ((!sym.compare("S")) || (!sym.compare("s")))
        {
            _reslice->SetOrientation(ResliceRender::SAGITTAL);
        }
        else if ((!sym.compare("O")) || (!sym.compare("o")))
        {
            _reslice->SetOrientation(ResliceRender::OBLIQUE);
        }
    }

    void KeyCallback::SetCallbackData(ResliceRender* reslice)
    {
        _reslice=reslice;
    }

    // Usage:  ResliceSphere [fileName]
    int main(int argc, char *argv[])
    {
        ResliceRender render;

        if (argc == 1)
          {
          const char *root = gdcm::Testing::GetDataExtraRoot();
          std::string dir3 = root;
          dir3 += "/gdcmSampleData/ForSeriesTesting/Dentist/images/";
          render.CreatePipeline(dir3.c_str());
          }
        else
          {
          render.CreatePipeline(argv[1]);
          }

        render.SetOrientation(ResliceRender::AXIAL);
        render.Start();

        return EXIT_SUCCESS;
    }