]> Creatis software - clitk.git/blobdiff - vv/vvSlicer.cxx
merge master
[clitk.git] / vv / vvSlicer.cxx
index 2255cce74c5b8a719ec671caaddbea425636eb4b..9809b06ddc771b574d4de2923c2a750ad11e1f67 100644 (file)
@@ -3,7 +3,7 @@
 
   Authors belong to:
   - University of LYON              http://www.universite-lyon.fr/
-  - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
+  - Léon Bérard cancer center       http://www.centreleonberard.fr
   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
 
   This software is distributed WITHOUT ANY WARRANTY; without even
@@ -14,7 +14,7 @@
 
   - BSD        See included LICENSE.txt file
   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-  ======================================================================-====*/
+  ===========================================================================**/
 
 #include "vvSlicer.h"
 #include "vvImage.h"
 #include <vtkImageAccumulate.h>
 #include <vtkImageReslice.h>
 
+// template <class T, unsigned int dim>
+// void print_vector(const char* pmsg, T* pvec)
+// {
+//   std::cout << pmsg << ": ";
+//   for (unsigned int i = 0; i < dim; i++)
+//     std::cout << pvec[i] << " ";
+//   std::cout << std::endl;
+// }
+
+
 vtkCxxRevisionMacro(vvSlicer, "DummyRevision");
 vtkStandardNewMacro(vvSlicer);
 
@@ -309,7 +319,18 @@ void vvSlicer::SetImage(vvImage::Pointer image)
 {
   if (image->GetVTKImages().size()) {
     mImage = image;
-    this->Superclass::SetInput(image->GetTransformedVTKImages()[0]);
+
+    if (!mImageReslice) {
+      mImageReslice = vtkSmartPointer<vtkImageReslice>::New();
+      mImageReslice->SetInterpolationModeToLinear();
+      mImageReslice->AutoCropOutputOn();
+      mImageReslice->SetBackgroundColor(-1000,-1000,-1000,1);
+    }
+    mImageReslice->SetResliceTransform(mImage->GetTransform());
+    mImageReslice->SetInput(0, mImage->GetFirstVTKImageData());
+    mImageReslice->UpdateInformation();
+
+    this->Superclass::SetInput(mImageReslice->GetOutput());
 
     int extent[6];
     this->GetInput()->GetWholeExtent(extent);
@@ -322,8 +343,8 @@ void vvSlicer::SetImage(vvImage::Pointer image)
     // Make sure that the required part image has been computed
     extent[SliceOrientation*2] = Slice;
     extent[SliceOrientation*2+1] = Slice;
-    image->GetTransformedVTKImages()[0]->SetUpdateExtent(extent);
-    image->GetTransformedVTKImages()[0]->Update();
+    mImageReslice->GetOutput()->SetUpdateExtent(extent);
+    mImageReslice->GetOutput()->Update();
 
     this->UpdateDisplayExtent();
 
@@ -339,10 +360,19 @@ void vvSlicer::SetOverlay(vvImage::Pointer overlay)
 {
   if (overlay->GetVTKImages().size()) {
     mOverlay = overlay;
+    
+    if (!mOverlayReslice) {
+      mOverlayReslice = vtkSmartPointer<vtkImageReslice>::New();
+      mOverlayReslice->SetInterpolationModeToLinear();
+      mOverlayReslice->AutoCropOutputOn();
+      mOverlayReslice->SetBackgroundColor(-1000,-1000,-1000,1);
+    }
+    mOverlayReslice->SetResliceTransform(mOverlay->GetTransform());
+    mOverlayReslice->SetInput(0, mOverlay->GetFirstVTKImageData());
 
     if (!mOverlayMapper)
       mOverlayMapper = vtkSmartPointer<vtkImageMapToWindowLevelColors>::New();
-    mOverlayMapper->SetInput(overlay->GetTransformedVTKImages()[0]);
+    mOverlayMapper->SetInput(mOverlayReslice->GetOutput());
 
     if (!mOverlayActor) {
       mOverlayActor = vtkSmartPointer<vvBlendImageActor>::New();
@@ -350,7 +380,6 @@ void vvSlicer::SetOverlay(vvImage::Pointer overlay)
       mOverlayActor->SetPickable(0);
       mOverlayActor->SetVisibility(false);
       mOverlayActor->SetOpacity(0.5);
-      this->UpdateDisplayExtent();
     }
 
     //stupid but necessary : the Overlay need to be rendered before fusion
@@ -361,8 +390,9 @@ void vvSlicer::SetOverlay(vvImage::Pointer overlay)
     } else
       this->GetRenderer()->AddActor(mOverlayActor);
 
-    //Synchronize slice
-    SetTSlice(mCurrentTSlice);
+    //Synchronize orientation and slice
+    this->SetSliceOrientation(this->SliceOrientation);
+    this->SetTSlice(mCurrentTSlice);
   }
 }
 //------------------------------------------------------------------------------
@@ -374,9 +404,18 @@ void vvSlicer::SetFusion(vvImage::Pointer fusion)
   if (fusion->GetVTKImages().size()) {
     mFusion = fusion;
 
+    if (!mFusionReslice) {
+      mFusionReslice = vtkSmartPointer<vtkImageReslice>::New();
+      mFusionReslice->SetInterpolationModeToLinear();
+      mFusionReslice->AutoCropOutputOn();
+      mFusionReslice->SetBackgroundColor(-1000,-1000,-1000,1);
+    }
+    mFusionReslice->SetResliceTransform(mFusion->GetTransform());
+    mFusionReslice->SetInput(0, mFusion->GetFirstVTKImageData());
+
     if (!mFusionMapper)
       mFusionMapper = vtkSmartPointer<vtkImageMapToWindowLevelColors>::New();
-    mFusionMapper->SetInput(fusion->GetTransformedVTKImages()[0]);
+    mFusionMapper->SetInput(mFusionReslice->GetOutput());
 
     if (!mFusionActor) {
       mFusionActor = vtkSmartPointer<vtkImageActor>::New();
@@ -384,12 +423,12 @@ void vvSlicer::SetFusion(vvImage::Pointer fusion)
       mFusionActor->SetPickable(0);
       mFusionActor->SetVisibility(false);
       mFusionActor->SetOpacity(0.7);
-      this->UpdateDisplayExtent();
       this->GetRenderer()->AddActor(mFusionActor);
     }
 
-    //Synchronize slice
-    SetTSlice(mCurrentTSlice);
+    //Synchronize orientation and slice
+    this->SetSliceOrientation(this->SliceOrientation);
+    this->SetTSlice(mCurrentTSlice);
   }
 }
 //------------------------------------------------------------------------------
@@ -425,7 +464,7 @@ void vvSlicer::SetVF(vvImage::Pointer vf)
       mVOIFilter = vtkSmartPointer<vtkExtractVOI>::New();
       mVOIFilter->SetSampleRate(mSubSampling,mSubSampling,mSubSampling);
     }
-    mVOIFilter->SetInput(vf->GetTransformedVTKImages()[0]);
+    mVOIFilter->SetInput(vf->GetFirstVTKImageData());
     mAAFilter->SetInput(mVOIFilter->GetOutput());
     ///This tells VTK to use the scalar (pixel) data of the image to draw the little arrows
     mAAFilter->Assign(vtkDataSetAttributes::SCALARS, vtkDataSetAttributes::VECTORS, vtkAssignAttribute::POINT_DATA);
@@ -620,18 +659,18 @@ void vvSlicer::SetTSlice(int t)
   if (mCurrentTSlice == t) return;
 
   mCurrentTSlice = t;
-  this->SetInput(mImage->GetTransformedVTKImages()[t]);
+  mImageReslice->SetInput( mImage->GetVTKImages()[mCurrentTSlice] );
   if (mVF && mVFActor->GetVisibility()) {
     if (mVF->GetVTKImages().size() > (unsigned int)mCurrentTSlice)
-      mVOIFilter->SetInput(mVF->GetTransformedVTKImages()[mCurrentTSlice]);
+      mVOIFilter->SetInput(mVF->GetVTKImages()[mCurrentTSlice]);
   }
   if (mOverlay && mOverlayActor->GetVisibility()) {
-    if (mOverlay->GetTransformedVTKImages().size() > (unsigned int)mCurrentTSlice)
-      mOverlayMapper->SetInput(mOverlay->GetTransformedVTKImages()[mCurrentTSlice]);
+    if (mOverlay->GetVTKImages().size() > (unsigned int)mCurrentTSlice)
+      mOverlayReslice->SetInput( mOverlay->GetVTKImages()[mCurrentTSlice] );
   }
   if (mFusion && mFusionActor->GetVisibility()) {
     if (mFusion->GetVTKImages().size() > (unsigned int)mCurrentTSlice)
-      mFusionMapper->SetInput(mFusion->GetTransformedVTKImages()[mCurrentTSlice]);
+      mFusionReslice->SetInput( mFusion->GetVTKImages()[mCurrentTSlice]);
   }
   if (mSurfaceCutActors.size() > 0)
     for (std::vector<vvMeshActor*>::iterator i=mSurfaceCutActors.begin();
@@ -649,7 +688,6 @@ int vvSlicer::GetTSlice()
 }
 //------------------------------------------------------------------------------
 
-
 //------------------------------------------------------------------------------
 void vvSlicer::SetSliceOrientation(int orientation)
 {
@@ -657,16 +695,22 @@ void vvSlicer::SetSliceOrientation(int orientation)
   int extent[6];
   this->GetInput()->GetWholeExtent(extent);
   if (extent[5]-extent[4] <= 2)
-    orientation=2;
+    orientation = vtkImageViewer2::SLICE_ORIENTATION_XY;
 
   if (orientation < vtkImageViewer2::SLICE_ORIENTATION_YZ ||
       orientation > vtkImageViewer2::SLICE_ORIENTATION_XY) {
     vtkErrorMacro("Error - invalid slice orientation " << orientation);
     return;
   }
-
+  
   this->SliceOrientation = orientation;
 
+  if(mFusion)
+    AdjustResliceToSliceOrientation(mFusionReslice);
+
+  if(mOverlay)
+    AdjustResliceToSliceOrientation(mOverlayReslice);
+
   // Update the viewer
   int *range = this->GetSliceRange();
   if (range)
@@ -690,6 +734,41 @@ void vvSlicer::SetSliceOrientation(int orientation)
 }
 //----------------------------------------------------------------------------
 
+//------------------------------------------------------------------------------
+// This function ensures that we sample the slices of a vtkImageReslice filter
+// in the direction of the slicer (SliceOrientation) similarly as mImageReslice.
+// In other words, we change the grid of the reslice in the same way as the grid
+// of the displayed image in the slicing direction.
+void vvSlicer::AdjustResliceToSliceOrientation(vtkImageReslice *reslice)
+{
+  // Reset autocrop
+  double origin[3] = {VTK_DOUBLE_MAX, VTK_DOUBLE_MAX, VTK_DOUBLE_MAX};
+  double spacing[3] = {VTK_DOUBLE_MAX, VTK_DOUBLE_MAX, VTK_DOUBLE_MAX};
+  reslice->SetOutputOrigin(origin);
+  reslice->SetOutputSpacing(spacing);
+  reslice->GetOutput()->UpdateInformation();
+  reslice->GetOutput()->GetOrigin(origin);
+  reslice->GetOutput()->GetSpacing(spacing);
+
+  // Use similar spacing as the image in the direction SliceOrientation
+  spacing[this->SliceOrientation] = mImageReslice->GetOutput()->GetSpacing()[this->SliceOrientation];
+
+  // Modify origin to be on the image grid in the direction SliceOrientation in 3 steps
+  // Step 1: from world coordinates to image coordinates
+  origin[this->SliceOrientation] -= mImageReslice->GetOutput()->GetOrigin()[this->SliceOrientation];
+  origin[this->SliceOrientation] /= mImageReslice->GetOutput()->GetSpacing()[this->SliceOrientation];
+  // Step 2: round to superior grid positionInc
+  origin[this->SliceOrientation] = itk::Math::Ceil<double>(origin[this->SliceOrientation]);
+  // Step 3: back to world coordinates
+  origin[this->SliceOrientation] *= mImageReslice->GetOutput()->GetSpacing()[this->SliceOrientation];
+  origin[this->SliceOrientation] += mImageReslice->GetOutput()->GetOrigin()[this->SliceOrientation];
+
+  // Set new spacing and origin
+  reslice->SetOutputOrigin(origin);
+  reslice->SetOutputSpacing(spacing);
+  reslice->UpdateInformation();
+}
+//------------------------------------------------------------------------------
 
 //----------------------------------------------------------------------------
 int * vvSlicer::GetExtent()
@@ -710,6 +789,7 @@ int vvSlicer::GetOrientation()
 }
 //----------------------------------------------------------------------------
 
+
 //----------------------------------------------------------------------------
 void vvSlicer::UpdateDisplayExtent()
 {
@@ -735,59 +815,70 @@ void vvSlicer::UpdateDisplayExtent()
   // Image actor
   this->ImageActor->SetDisplayExtent(w_ext);
   
-  // Position vector
-  double position[3] = {0.,0.,0.};
-  double positionInc = (Renderer->GetActiveCamera()->GetPosition()[this->SliceOrientation] > this->Slice)?10:-10;
-  position[this->SliceOrientation] += positionInc;
-  
   // Overlay image actor
   if (mOverlay && mOverlayActor->GetVisibility()) {
     int overExtent[6];
-    mOverlay->GetTransformedVTKImages()[0]->UpdateInformation();
-    this->ConvertImageToImageDisplayExtent(input, w_ext, mOverlay->GetTransformedVTKImages()[0], overExtent);
+    mOverlayReslice->GetOutput()->UpdateInformation();
+    this->ConvertImageToImageDisplayExtent(input, w_ext, mOverlayReslice->GetOutput(), overExtent);
     ClipDisplayedExtent(overExtent, mOverlayMapper->GetInput()->GetWholeExtent());
     mOverlayActor->SetDisplayExtent( overExtent );
-    mOverlayActor->SetPosition(position);
   }
-  position[this->SliceOrientation] += positionInc;
 
   // Fusion image actor
   if (mFusion && mFusionActor->GetVisibility()) {
     int fusExtent[6];
-    mFusion->GetTransformedVTKImages()[0]->UpdateInformation();
-    this->ConvertImageToImageDisplayExtent(input, w_ext, mFusion->GetTransformedVTKImages()[0], fusExtent);
+    mFusionReslice->GetOutput()->UpdateInformation();
+    this->ConvertImageToImageDisplayExtent(input, w_ext, mFusionReslice->GetOutput(), fusExtent);
     ClipDisplayedExtent(fusExtent, mFusionMapper->GetInput()->GetWholeExtent());
     mFusionActor->SetDisplayExtent(fusExtent);
-    mFusionActor->SetPosition(position);
   }
-  position[this->SliceOrientation] += positionInc;
 
   // Vector field actor
+  double* camera = Renderer->GetActiveCamera()->GetPosition();
+  double* image_bounds = ImageActor->GetBounds();
+  double position[3] = {0, 0, 0};
+  position[this->SliceOrientation] = image_bounds[this->SliceOrientation*2]; 
+
+  //print_vector<double, 6>("camera", camera);
+  //print_vector<double, 6>("image_bounds", image_bounds);
+  //print_vector<double, 3>("position", position);
+
+  // find where to place the VF actor. to deal with
+  // z-buffer issues, the VF is placed right in front of the image,
+  // subject to a small offset. the position actually depends on the
+  // the location of the camera relative to the image. 
+  double offset = 1;
+  if (camera[this->SliceOrientation] < image_bounds[this->SliceOrientation*2])
+    offset = -1;
+  
   if (mVF && mVFActor->GetVisibility()) {
     int vfExtent[6];
-    mVF->GetTransformedVTKImages()[0]->UpdateInformation();
-    this->ConvertImageToImageDisplayExtent(input, w_ext, mVF->GetTransformedVTKImages()[0], vfExtent);
+    mVF->GetVTKImages()[0]->UpdateInformation();
+    this->ConvertImageToImageDisplayExtent(input, w_ext, mVF->GetVTKImages()[0], vfExtent);
     ClipDisplayedExtent(vfExtent, mVOIFilter->GetInput()->GetWholeExtent());
     mVOIFilter->SetVOI(vfExtent);
     int orientation[3] = {1,1,1};
     orientation[this->SliceOrientation] = 0;
     mGlyphFilter->SetOrientation(orientation[0], orientation[1], orientation[2]);
     mVFMapper->Update();
+
+    position[this->SliceOrientation] += offset;
     mVFActor->SetPosition(position);
   }
-  position[this->SliceOrientation] += positionInc;
-
+  
   // Landmarks actor
   if (mLandActor) {
     if (mClipBox) {
       double bounds [6];
       for(unsigned int i=0; i<6; i++)
         bounds[i] = ImageActor->GetBounds()[i];
-      bounds[ this->SliceOrientation*2   ] = ImageActor->GetBounds()[ this->SliceOrientation*2  ]-fabs(0.5/this->GetInput()->GetSpacing()[this->SliceOrientation]);
-      bounds[ this->SliceOrientation*2+1 ] = ImageActor->GetBounds()[ this->SliceOrientation*2+1 ]+fabs(0.5/this->GetInput()->GetSpacing()[this->SliceOrientation]);
+      bounds[ this->SliceOrientation*2   ] = ImageActor->GetBounds()[ this->SliceOrientation*2  ]-fabs(this->GetInput()->GetSpacing()[this->SliceOrientation]);
+      bounds[ this->SliceOrientation*2+1 ] = ImageActor->GetBounds()[ this->SliceOrientation*2+1 ]+fabs(this->GetInput()->GetSpacing()[this->SliceOrientation]);
       mClipBox->SetBounds(bounds);
       UpdateLandmarks();
     }
+    
+    position[this->SliceOrientation] = offset;
     mLandActor->SetPosition(position);
   }
 
@@ -811,6 +902,7 @@ void vvSlicer::UpdateDisplayExtent()
       }
     }
   }
+  
 }
 //----------------------------------------------------------------------------
 
@@ -827,7 +919,7 @@ void vvSlicer::ConvertImageToImageDisplayExtent(vtkImageData *sourceImage, const
     dExtents[i] = (dExtents[i]- targetImage->GetOrigin()[i/2]) / targetImage->GetSpacing()[i/2];
     
     // Round to nearest
-    targetExtent[i] = itk::Math::Round(dExtents[i]);
+    targetExtent[i] = itk::Math::Round<double>(dExtents[i]);
   }
 }
 //----------------------------------------------------------------------------
@@ -1081,7 +1173,7 @@ void vvSlicer::GetExtremasAroundMousePointer(double & min, double & max)
 //----------------------------------------------------------------------------
 
 //----------------------------------------------------------------------------
-double vvSlicer::GetScalarComponentAsDouble(vtkImageData *image, int X, double Y, double Z, int &ix, int &iy, int &iz, int component)
+double vvSlicer::GetScalarComponentAsDouble(vtkImageData *image, double X, double Y, double Z, int &ix, int &iy, int &iz, int component)
 {
   ix = lrint(X);
   iy = lrint(Y);
@@ -1092,7 +1184,7 @@ double vvSlicer::GetScalarComponentAsDouble(vtkImageData *image, int X, double Y
       iy > image->GetWholeExtent()[3] ||
       iz < image->GetWholeExtent()[4] ||
       iz > image->GetWholeExtent()[5] )
-    return sqrt(-1.);
+    return std::numeric_limits<double>::quiet_NaN();
 
   image->SetUpdateExtent(ix, ix, iy, iy, iz, iz);
   image->Update();
@@ -1109,13 +1201,7 @@ void vvSlicer::Render()
   } else legend->SetVisibility(0);
 
   if (ca->GetVisibility()) {
-    std::string worldPos = "";
-    std::stringstream world1;
-    std::stringstream world2;
-    std::stringstream world3;
-    world1 << (int)mCurrent[0];
-    world2 << (int)mCurrent[1];
-    world3 << (int)mCurrent[2];
+    std::stringstream worldPos;
     double X = (mCurrent[0] - this->GetInput()->GetOrigin()[0])/this->GetInput()->GetSpacing()[0];
     double Y = (mCurrent[1] - this->GetInput()->GetOrigin()[1])/this->GetInput()->GetSpacing()[1];
     double Z = (mCurrent[2] - this->GetInput()->GetOrigin()[2])/this->GetInput()->GetSpacing()[2];
@@ -1138,24 +1224,19 @@ void vvSlicer::Render()
       int ix, iy, iz;
       double value = this->GetScalarComponentAsDouble(this->GetInput(), X, Y, Z, ix, iy, iz);
 
-      std::stringstream pixel1;
-      std::stringstream pixel2;
-      std::stringstream pixel3;
-      std::stringstream temps;
-      pixel1 << ix;
-      pixel2 << iy;
-      pixel3 << iz;
-      temps << mCurrentTSlice;
-
-      std::stringstream val;
-      val << value;
-      worldPos += "data value : " + val.str();
-      worldPos += "\n mm : " + world1.str() + " " + world2.str() + " " +
-        world3.str() + " " + temps.str();
-      worldPos += "\n pixel : " + pixel1.str() + " " + pixel2.str() + " " +
-        pixel3.str() + " " + temps.str();
+      worldPos << "data value : " << value << std::endl;
+      worldPos << "mm : " << lrint(mCurrent[0]) << ' '
+                          << lrint(mCurrent[1]) << ' '
+                          << lrint(mCurrent[2]) << ' '
+                          << mCurrentTSlice
+                          << std::endl;
+      worldPos << "pixel : " << ix << ' '
+                             << iy << ' '
+                             << iz << ' '
+                             << mCurrentTSlice
+                             << std::endl;
     }
-    ca->SetText(1,worldPos.c_str());
+    ca->SetText(1,worldPos.str().c_str());
   }
 
   if (pdmA->GetVisibility()) {