# gdcmorthoplanes.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.

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

    #include "vtkActor.h"
    #include "vtkCamera.h"
    #include "vtkMatrix4x4.h"
    #include "vtkTransform.h"
    #include "vtkAssembly.h"
    #include "vtkCellPicker.h"
    #include "vtkCommand.h"
    #include "vtkImageActor.h"
    #include "vtkImageMapToColors.h"
    #include "vtkImageOrthoPlanes.h"
    #include "vtkImagePlaneWidget.h"
    #include "vtkImageReader.h"
    #include "vtkInteractorEventRecorder.h"
    #include "vtkLookupTable.h"
    #include "vtkOutlineFilter.h"
    #include "vtkPolyDataMapper.h"
    #include "vtkProperty.h"
    #include "vtkRenderWindow.h"
    #include "vtkRenderWindowInteractor.h"
    #include "vtkRenderer.h"
    #include "vtkVolume16Reader.h"
    #include "vtkImageData.h"
    #include "vtkImageChangeInformation.h"
    #include "vtkOrientationMarkerWidget.h"
    #include "vtkAnnotatedCubeActor.h"
    #include "vtkAxesActor.h"
    #include "vtkCaptionActor2D.h"
    #include "vtkTextProperty.h"
    #include "vtkPropAssembly.h"

    #include "vtkGDCMImageReader.h"
    #include "vtkGDCMImageWriter.h"
    #include "vtkStringArray.h"

    #include "gdcmSystem.h"
    #include "gdcmDirectory.h"
    #include "gdcmIPPSorter.h"

    #ifndef vtkFloatingPointType
    #define vtkFloatingPointType float
    #endif

    //----------------------------------------------------------------------------
    class vtkOrthoPlanesCallback : public vtkCommand
    {
    public:
      static vtkOrthoPlanesCallback *New()
      { return new vtkOrthoPlanesCallback; }

      void Execute( vtkObject *caller, unsigned long vtkNotUsed( event ),
                    void *callData )
      {
        vtkImagePlaneWidget* self =
          reinterpret_cast< vtkImagePlaneWidget* >( caller );
        if(!self) return;

        double* wl = static_cast<double*>( callData );

        if ( self == this->WidgetX )
          {
          this->WidgetY->SetWindowLevel(wl[0],wl[1],1);
          this->WidgetZ->SetWindowLevel(wl[0],wl[1],1);
          }
        else if( self == this->WidgetY )
          {
          this->WidgetX->SetWindowLevel(wl[0],wl[1],1);
          this->WidgetZ->SetWindowLevel(wl[0],wl[1],1);
          }
        else if (self == this->WidgetZ)
          {
          this->WidgetX->SetWindowLevel(wl[0],wl[1],1);
          this->WidgetY->SetWindowLevel(wl[0],wl[1],1);
          }
      }

      vtkOrthoPlanesCallback():WidgetX( 0 ), WidgetY( 0 ), WidgetZ ( 0 ) {}

      vtkImagePlaneWidget* WidgetX;
      vtkImagePlaneWidget* WidgetY;
      vtkImagePlaneWidget* WidgetZ;
    };

    int main( int argc, char *argv[] )
    {
      //char* fname = vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/headsq/quarter");

      //vtkVolume16Reader* v16 =  vtkVolume16Reader::New();
      //  v16->SetDataDimensions( 64, 64);
      //  v16->SetDataByteOrderToLittleEndian();
      //  v16->SetImageRange( 1, 93);
      //  v16->SetDataSpacing( 3.2, 3.2, 1.5);
      //  v16->SetFilePrefix( fname );
      //  v16->SetDataMask( 0x7fff);
      //  v16->Update();
      std::vector<std::string> filenames;
      if( argc < 2 )
        {
        std::cerr << argv[0] << " filename1.dcm [filename2.dcm ...]\n";
        return 1;
        }
      else
        {
        // Is it a single directory ? If so loop over all files contained in it:
        const char *filename = argv[1];
        if( argc == 2 && gdcm::System::FileIsDirectory( filename ) )
          {
          std::cout << "Loading directory: " << filename << std::endl;
          bool recursive = false;
          gdcm::Directory d;
          d.Load(filename, recursive);
          gdcm::Directory::FilenamesType const &files = d.GetFilenames();
          for( gdcm::Directory::FilenamesType::const_iterator it = files.begin(); it != files.end(); ++it )
            {
            filenames.push_back( it->c_str() );
            }
          }
        else // list of files passed directly on the cmd line:
            // discard non-existing or directory
          {
          for(int i=1; i < argc; ++i)
            {
            filename = argv[i];
            if( gdcm::System::FileExists( filename ) )
              {
              if( gdcm::System::FileIsDirectory( filename ) )
                {
                std::cerr << "Discarding directory: " << filename << std::endl;
                }
              else
                {
                filenames.push_back( filename );
                }
              }
            else
              {
              std::cerr << "Discarding non existing file: " << filename << std::endl;
              }
            }
          }
        //names->Print( std::cout );
        }

      //gdcm::Trace::DebugOn();
      //gdcm::Trace::WarningOn();
      gdcm::IPPSorter s;
      s.SetComputeZSpacing( true );
      s.SetZSpacingTolerance( 1e-3 );
      bool b = s.Sort( filenames );
      if( !b )
        {
        //std::cerr << "Failed to sort:" << directory << std::endl;
        return 1;
        }
      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 *files = vtkStringArray::New();
      std::vector< std::string >::const_iterator it = sorted.begin();
      for( ; it != sorted.end(); ++it)
        {
        const std::string &f = *it;
        files->InsertNextValue( f.c_str() );
        }

      //delete[] fname;
      vtkGDCMImageReader * reader = vtkGDCMImageReader::New();
      //reader->SetFileLowerLeft( 1 );
      reader->SetFileNames( files );
      reader->Update(); // important
      //reader->GetOutput()->Print( std::cout );
      //vtkFloatingPointType range[2];
      //reader->GetOutput()->GetScalarRange(range);
      //std::cout << "Range: " << range[0] << " " << range[1] << std::endl;

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

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

    #if 0
        vtkGDCMImageWriter *writer = vtkGDCMImageWriter::New();
        writer->SetInput( v16->GetOutput() );
        writer->SetFileLowerLeft( reader->GetFileLowerLeft() );
        writer->SetDirectionCosines( reader->GetDirectionCosines() );
        writer->SetImageFormat( reader->GetImageFormat() );
        writer->SetFileDimensionality( 3); //reader->GetFileDimensionality() );
        writer->SetMedicalImageProperties( reader->GetMedicalImageProperties() );
        writer->SetShift( reader->GetShift() );
        writer->SetScale( reader->GetScale() );
        writer->SetFileName( "out.dcm" );
        writer->Write();
    #endif

      files->Delete();

      vtkOutlineFilter* outline = vtkOutlineFilter::New();
        outline->SetInputConnection(v16->GetOutputPort());

      vtkPolyDataMapper* outlineMapper = vtkPolyDataMapper::New();
        outlineMapper->SetInputConnection(outline->GetOutputPort());

      vtkActor* outlineActor =  vtkActor::New();
        outlineActor->SetMapper( outlineMapper);

      vtkRenderer* ren1 = vtkRenderer::New();
      vtkRenderer* ren2 = vtkRenderer::New();

      vtkRenderWindow* renWin = vtkRenderWindow::New();
        renWin->AddRenderer(ren2);
        renWin->AddRenderer(ren1);

      vtkRenderWindowInteractor* iren = vtkRenderWindowInteractor::New();
        iren->SetRenderWindow(renWin);

      vtkCellPicker* picker = vtkCellPicker::New();
        picker->SetTolerance(0.005);

      vtkProperty* ipwProp = vtkProperty::New();
       //assign default props to the ipw's texture plane actor

      vtkImagePlaneWidget* planeWidgetX = vtkImagePlaneWidget::New();
        planeWidgetX->SetInteractor( iren);
        planeWidgetX->SetKeyPressActivationValue('x');
        planeWidgetX->SetPicker(picker);
        planeWidgetX->RestrictPlaneToVolumeOn();
        planeWidgetX->GetPlaneProperty()->SetColor(1,0,0);
        planeWidgetX->SetTexturePlaneProperty(ipwProp);
        planeWidgetX->TextureInterpolateOff();
        planeWidgetX->SetResliceInterpolateToNearestNeighbour();
        planeWidgetX->SetInput(v16->GetOutput());
        planeWidgetX->SetPlaneOrientationToXAxes();
        //planeWidgetX->SetSliceIndex(32);
        planeWidgetX->DisplayTextOn();
        planeWidgetX->On();
        planeWidgetX->InteractionOff();
        planeWidgetX->InteractionOn();

      vtkImagePlaneWidget* planeWidgetY = vtkImagePlaneWidget::New();
        planeWidgetY->SetInteractor( iren);
        planeWidgetY->SetKeyPressActivationValue('y');
        planeWidgetY->SetPicker(picker);
        planeWidgetY->GetPlaneProperty()->SetColor(1,1,0);
        planeWidgetY->SetTexturePlaneProperty(ipwProp);
        planeWidgetY->TextureInterpolateOn();
        planeWidgetY->SetResliceInterpolateToLinear();
        planeWidgetY->SetInput(v16->GetOutput());
        planeWidgetY->SetPlaneOrientationToYAxes();
        //planeWidgetY->SetSlicePosition(102.4);
        planeWidgetY->SetLookupTable( planeWidgetX->GetLookupTable());
        planeWidgetY->DisplayTextOn();
        planeWidgetY->UpdatePlacement();
        planeWidgetY->On();

      vtkImagePlaneWidget* planeWidgetZ = vtkImagePlaneWidget::New();
        planeWidgetZ->SetInteractor( iren);
        planeWidgetZ->SetKeyPressActivationValue('z');
        planeWidgetZ->SetPicker(picker);
        planeWidgetZ->GetPlaneProperty()->SetColor(0,0,1);
        planeWidgetZ->SetTexturePlaneProperty(ipwProp);
        planeWidgetZ->TextureInterpolateOn();
        planeWidgetZ->SetResliceInterpolateToCubic();
        planeWidgetZ->SetInput(v16->GetOutput());
        planeWidgetZ->SetPlaneOrientationToZAxes();
        //planeWidgetZ->SetSliceIndex(25);
        planeWidgetZ->SetLookupTable( planeWidgetX->GetLookupTable());
        planeWidgetZ->DisplayTextOn();
        planeWidgetZ->On();

      vtkImageOrthoPlanes *orthoPlanes = vtkImageOrthoPlanes::New();
        orthoPlanes->SetPlane(0, planeWidgetX);
        orthoPlanes->SetPlane(1, planeWidgetY);
        orthoPlanes->SetPlane(2, planeWidgetZ);
        orthoPlanes->ResetPlanes();

       vtkOrthoPlanesCallback* cbk = vtkOrthoPlanesCallback::New();
       cbk->WidgetX = planeWidgetX;
       cbk->WidgetY = planeWidgetY;
       cbk->WidgetZ = planeWidgetZ;
       planeWidgetX->AddObserver( vtkCommand::EndWindowLevelEvent, cbk );
       planeWidgetY->AddObserver( vtkCommand::EndWindowLevelEvent, cbk );
       planeWidgetZ->AddObserver( vtkCommand::EndWindowLevelEvent, cbk );
       cbk->Delete();

      double wl[2];
      planeWidgetZ->GetWindowLevel(wl);

      // Add a 2D image to test the GetReslice method
      //
      vtkImageMapToColors* colorMap = vtkImageMapToColors::New();
        colorMap->PassAlphaToOutputOff();
        colorMap->SetActiveComponent(0);
        colorMap->SetOutputFormatToLuminance();
        colorMap->SetInput(planeWidgetZ->GetResliceOutput());
        colorMap->SetLookupTable(planeWidgetX->GetLookupTable());

      vtkImageActor* imageActor = vtkImageActor::New();
        imageActor->PickableOff();
        imageActor->SetInput(colorMap->GetOutput());

      // Add the actors
      //
      ren1->AddActor( outlineActor);
      ren2->AddActor( imageActor);

      ren1->SetBackground( 0.1, 0.1, 0.2);
      ren2->SetBackground( 0.2, 0.1, 0.2);

      renWin->SetSize( 600, 350);

      ren1->SetViewport(0,0,0.58333,1);
      ren2->SetViewport(0.58333,0,1,1);

      // Set the actors' postions
      //
      renWin->Render();
      //iren->SetEventPosition( 175,175);
      //iren->SetKeyCode('r');
      //iren->InvokeEvent(vtkCommand::CharEvent,NULL);
      //iren->SetEventPosition( 475,175);
      //iren->SetKeyCode('r');
      //iren->InvokeEvent(vtkCommand::CharEvent,NULL);
      //renWin->Render();

      //ren1->GetActiveCamera()->Elevation(110);
      //ren1->GetActiveCamera()->SetViewUp(0, 0, -1);
      //ren1->GetActiveCamera()->Azimuth(45);
      //ren1->GetActiveCamera()->Dolly(1.15);
      ren1->ResetCameraClippingRange();

      vtkAnnotatedCubeActor* cube = vtkAnnotatedCubeActor::New();
      cube->SetXPlusFaceText ( "R" );
      cube->SetXMinusFaceText( "L" );
      cube->SetYPlusFaceText ( "A" );
      cube->SetYMinusFaceText( "P" );
      cube->SetZPlusFaceText ( "H" );
      cube->SetZMinusFaceText( "F" );
      cube->SetFaceTextScale( 0.666667 );

      vtkAxesActor* axes2 = vtkAxesActor::New();

      vtkMatrix4x4 *invert = vtkMatrix4x4::New();
      invert->DeepCopy( reader->GetDirectionCosines() );
      invert->Invert();

      // simulate a left-handed coordinate system
      //
      vtkTransform *transform = vtkTransform::New();
      transform->Identity();
      //transform->RotateY(90);
      transform->Concatenate(invert);
      axes2->SetShaftTypeToCylinder();
      axes2->SetUserTransform( transform );
      cube->GetAssembly()->SetUserTransform( transform );

      axes2->SetTotalLength( 1.5, 1.5, 1.5 );
      axes2->SetCylinderRadius( 0.500 * axes2->GetCylinderRadius() );
      axes2->SetConeRadius    ( 1.025 * axes2->GetConeRadius() );
      axes2->SetSphereRadius  ( 1.500 * axes2->GetSphereRadius() );

      vtkTextProperty* tprop = axes2->GetXAxisCaptionActor2D()->
        GetCaptionTextProperty();
      tprop->ItalicOn();
      tprop->ShadowOn();
      tprop->SetFontFamilyToTimes();

      axes2->GetYAxisCaptionActor2D()->GetCaptionTextProperty()->ShallowCopy( tprop );
      axes2->GetZAxisCaptionActor2D()->GetCaptionTextProperty()->ShallowCopy( tprop );

      vtkPropAssembly* assembly = vtkPropAssembly::New();
      assembly->AddPart( axes2 );
      assembly->AddPart( cube );

      vtkOrientationMarkerWidget* widget = vtkOrientationMarkerWidget::New();
      widget->SetOutlineColor( 0.9300, 0.5700, 0.1300 );
      widget->SetOrientationMarker( assembly );
      widget->SetInteractor( iren );
      widget->SetViewport( 0.0, 0.0, 0.4, 0.4 );
      widget->SetEnabled( 1 );
      widget->InteractiveOff();
      widget->InteractiveOn();

      // Playback recorded events
      //
      //vtkInteractorEventRecorder *recorder = vtkInteractorEventRecorder::New();
      //recorder->SetInteractor(iren);
      //recorder->ReadFromInputStringOn();
      //recorder->SetInputString(IOPeventLog);

      // Interact with data
      // Render the image
      //
      iren->Initialize();
      renWin->Render();

      // Test SetKeyPressActivationValue for one of the widgets
      //
      //iren->SetKeyCode('z');
      //iren->InvokeEvent(vtkCommand::CharEvent,NULL);
      //iren->SetKeyCode('z');
      //iren->InvokeEvent(vtkCommand::CharEvent,NULL);

      //int retVal = vtkRegressionTestImage( renWin );
      //
      //if ( retVal == vtkRegressionTester::DO_INTERACTOR)
        {
        iren->Start();
        }

      // Clean up
      //
      //recorder->Off();
      //recorder->Delete();

      ipwProp->Delete();
      orthoPlanes->Delete();
      planeWidgetX->Delete();
      planeWidgetY->Delete();
      planeWidgetZ->Delete();
      colorMap->Delete();
      imageActor->Delete();
      picker->Delete();
      outlineActor->Delete();
      outlineMapper->Delete();
      outline->Delete();
      iren->Delete();
      renWin->Delete();
      ren1->Delete();
      ren2->Delete();
      v16->Delete();
      reader->Delete();

      return 0;
    }