vtkImagePolyDataSeedConnectivity.cxx

Go to the documentation of this file.

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

  Program:   wxMaracas
  Module:    $RCSfile: vtkImagePolyDataSeedConnectivity.cxx,v $
  Language:  C++
  Date:      $Date: 2009/05/14 13:54:57 $
  Version:   $Revision: 1.1 $

  Copyright: (c) 2002, 2003
  License:
  
     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 "vtkImagePolyDataSeedConnectivity.h"

#include <vtkPolyData.h>
#include <vtkImageData.h>
#include <vtkMath.h>
#include <vtkPlane.h>
#include <vtkSphere.h>
#include <vtkImplicitBoolean.h>
#include <vtkImplicitFunctionToImageStencil.h>
#include <vtkImageStencil.h>
#include <vtkObjectFactory.h>
#include <vtkImageSeedConnectivity.h>
#include <vtkImageIterator.h>
#include <vtkImageGaussianSmooth.h>
#include <vtkMarchingCubes.h>
#include <vtkPolyDataConnectivityFilter.h>
#include <vtkCleanPolyData.h>
#include <vtkTriangleFilter.h>
#include <vtkPointData.h>
#include <vtkImageContinuousDilate3D.h>

vtkCxxRevisionMacro(vtkImagePolyDataSeedConnectivity, "$Revision: 1.1 $");
vtkStandardNewMacro(vtkImagePolyDataSeedConnectivity);
vtkCxxSetObjectMacro(vtkImagePolyDataSeedConnectivity, Axis,vtkPolyData);

//----------------------------------------------------------------------------
vtkImagePolyDataSeedConnectivity::vtkImagePolyDataSeedConnectivity()
{
  this->Axis = NULL;
  this->ClipImageData = vtkImageData::New();
  this->OuterMold = vtkPolyData::New();
  this->ThresholdRatio = 0.45;
}
//----------------------------------------------------------------------------
vtkImagePolyDataSeedConnectivity::~vtkImagePolyDataSeedConnectivity()
{
  if(this->Axis)
  {
    this->Axis->Delete();
  }
  this->ClipImageData->Delete();
  this->ClipImageData = NULL;
  this->OuterMold->Delete();
  this->OuterMold = NULL;
}
//Macros definitions
#define PSC_MIN(x,y)      ((x)<(y) ? (x) : (y))
#define PSC_MAX(x,y)      ((x)>(y) ? (x) : (y))
//----------------------------------------------------------------------------
void vtkImagePolyDataSeedConnectivity::Execute()
{
  vtkImageData *input = this->GetInput();
  vtkPolyData *output = this->GetOutput();
  int intervalle = 2;
  int dim[3];

  vtkDebugMacro(<< "Executing polydata seed connectivity...");

  if (input == NULL)
    {
    vtkErrorMacro(<<"Input is NULL");
    return;
    }

  // just deal with volumes
  if ( input->GetDataDimension() != 3 )
    {
    vtkErrorMacro("Bad input: only treats 3D structured point datasets");
    return;
    }

  input->GetDimensions(dim);

  //Allocate a temporary image that will be filled in b&w
  //this is the base of the algothm, from a graysccale input + an axis
  //we produce a b&w image the the marching cubes can contour:

  vtkImageData *imagedata = vtkImageData::New();
  imagedata->SetOrigin( input->GetOrigin() );
  imagedata->SetSpacing( input->GetSpacing() );
  imagedata->SetExtent( input->GetExtent() );
  imagedata->SetScalarTypeToUnsignedChar (); //for connectityseed + reduce mem usage
  imagedata->AllocateScalars();

  /*
  Result: piece of crap !!
  There was some leakage from another vessel that was near enough from the heart.
  */

  vtkImageSeedConnectivity *seed = vtkImageSeedConnectivity::New();
  seed->SetInput( imagedata );
  seed->SetInputConnectValue( 255 ); //If 0 -> inverse video, 255 -> rien !

  // Before starting the loop we should make sure the image has been clipped
  // with the polydata:
  this->ClipImageWithAxis();

  vtkPoints *points = this->Axis->GetPoints();
  vtkDataArray* scalars = this->Axis->GetPointData()->GetScalars();
  
  int numPts = points->GetNumberOfPoints();
  double p1[3], p2[3], p3[3], tuple[2];
  double radius, signalvalue;
  int region[6];
  int i_p2; //allow saving of p2 index
  unsigned short* inSI, *inSIEnd;
  unsigned char*  outSI,*outSIEnd;
    
  for(int i=0; i<numPts; i+=intervalle)
  {
    i_p2 = i+intervalle/2;
    points->GetPoint(i, p1);
    if(numPts <= i+2*intervalle)
    {
      i_p2 = i+intervalle;
      points->GetPoint(numPts - 1, p3);
      i = numPts;       //no more iteration
    }
    else
    {
      points->GetPoint(i+intervalle, p3);
    }
    points->GetPoint(i_p2, p2);
    
    //Now we can add p2 point to the see list:
    //seed->AddSeed( 3, p2); //too bad :(
    seed->AddSeed( (int)p2[0], (int)p2[1], (int)p2[2] );

    //this parameter seems to be ok for most image
    //furthermore it doesn't make any change to make it bigger.
    scalars->GetTuple(i_p2, tuple);
    radius              = tuple[0]; 
        signalvalue = tuple[1];
    radius *= 4;                //multiply for a factor 4 for aneurism case

    region[0] = (int)PSC_MAX(p2[0] - radius, 0);
    region[2] = (int)PSC_MAX(p2[1] - radius, 0);
    region[4] = (int)PSC_MAX(p2[2] - radius, 0);

    region[1] = (int)PSC_MIN(p2[0] + radius, dim[0] - 1);
    region[3] = (int)PSC_MIN(p2[1] + radius, dim[1] - 1);
    region[5] = (int)PSC_MIN(p2[2] + radius, dim[2] - 1);

    // multiply signal by a factor this allow user to accuratley find the best
    // range
    signalvalue *= ThresholdRatio;  

      //Solution: this is too much a trouble: rather do it once on the entry (=first  axis point)

      vtkImageIterator<unsigned short> inIt(this->ClipImageData, region);//image_woheart
      vtkImageIterator<unsigned char> outIt(imagedata, region);
      // vtkImageProgressIterator should be nicer in conjonction with a wxGauge

      // Loop through ouput pixels
      while (!inIt.IsAtEnd())
      {
        inSI = inIt.BeginSpan(); inSIEnd = inIt.EndSpan();
        outSI = outIt.BeginSpan(); outSIEnd = outIt.EndSpan();

        // Pixel operation
        while (inSI != inSIEnd)
          *outSI++ = (signalvalue <= *inSI++) ? 255 : 0;

        inIt.NextSpan(); outIt.NextSpan();
      }
   }
   
  //NEW: introduced by Maciek/Marcela: use a vtkImageGaussian to get rid 
  //of the crenelage/steps introduced by MarchingCubes
  
  vtkImageGaussianSmooth *gauss = vtkImageGaussianSmooth::New();
  gauss->SetInput( seed->GetOutput() );
  gauss->SetDimensionality(3);
  //gauss->SetStandardDeviation(1.0);
  //gauss->SetRadiusFactor(3.);

  vtkMarchingCubes *contour = vtkMarchingCubes::New();
  contour->SetInput ( gauss->GetOutput() );
  contour->SetNumberOfContours( 1 );
  contour->SetValue(0, 128);
  contour->ComputeGradientsOn ();
  contour->ComputeScalarsOn ();

  vtkCleanPolyData *clean = vtkCleanPolyData::New();
  clean->SetInput ( contour->GetOutput() );

  vtkTriangleFilter *tf = vtkTriangleFilter::New();
  tf->SetInput( clean->GetOutput() );
  tf->Update();

  //Need to do a deep copy as we'll do a ShallowCopy afterwards
  output->DeepCopy( tf->GetOutput() );

  //output->Squeeze();  //usefull ?

  // Now we do a dilatation of seed->GetOutput() in order to have inner (=output)
  // and outer mold (=OuterMold):
  //First apply a morpho math: dilatation of 2mm:
  vtkImageContinuousDilate3D *dilate = vtkImageContinuousDilate3D ::New();
  dilate->SetInput ( seed->GetOutput() );
  //TODO: determine real kernel based on scale = cm/pixel
  //See: http://public.kitware.com/pipermail/vtkusers/2000-September/004210.html
  //[Spacing is x and y is the distance between the samples. The units of measure
  //are your choice. Often medical images use millimeters. Let's say your medical
  //data resolution is 512 x 512 and your study has a 250 mm field of view. Then
  //the spacing in x and y would be 250/512 or .488 mm.]
  //therefore we have here:

  int kernelsize[3];    
  double spacing[3];
  seed->GetOutput()->GetSpacing( spacing );  //beware input should be Update()
  for(int j=0; j<3; j++)
  {
    //2mm thickness, we assume spacing units is in mm !!
    kernelsize[j] = 3; //(int) (2. / spacing[i]);
    //std::cout << "kernel=" << kernelsize[i] << ",et:" << (int)(2./0.3334) << ",";
    //beware if spacing=0.3334 -> 2/0.3334 = 5 !!!! bad bad +/-0.5 to see
  }
  dilate->SetKernelSize (kernelsize[0], kernelsize[1], kernelsize[2]);

  //VTK filters rules !!!
  gauss->SetInput( dilate->GetOutput());
  tf->Update();
  // shallow copy result into output.
  // Could be re-written, as we know there is no vert neither strips...
  OuterMold->ShallowCopy( tf->GetOutput() );
  
  //OuterMold->Squeeze();  //usefull ?

  //Delete temp stuff:
  imagedata             -> Delete();
  seed                  -> Delete();
  gauss                 -> Delete();
  contour               -> Delete();
  clean                 -> Delete();
  tf                    -> Delete();
  dilate                -> Delete();
  
  //TODO: if this filter is used many times when working on different image size
  //I guess we need to find something brighter for this->ClipImageData
}
//end of macros use
#undef PSC_MAX
#undef PSC_MIN
//----------------------------------------------------------------------------
//We need to add a function that from the raw data is able to clip with
//the begining and ending of polydata to produce sharp edge.
void vtkImagePolyDataSeedConnectivity::ClipImageWithAxis()
{
  //no input need we'll work with member data:
  //Now we only need the region growing part:
  double normal[3];

  this->Axis->Update();
  vtkPoints *points = this->Axis->GetPoints();
  int numPts = points->GetNumberOfPoints();
  double p1[3];
  double p2[3];
  
  points->GetPoint( 0, p2 );
  points->GetPoint( 1, p1 );
  
  normal[0] = p1[0] - p2[0];
  normal[1] = p1[1] - p2[1];
  normal[2] = p1[2] - p2[2];
  vtkMath::Normalize( normal );

  vtkPlane *plane1 = vtkPlane::New();
  plane1->SetOrigin( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
  plane1->SetNormal( normal[0], normal[1], normal[2]);

  vtkDataArray* scalars = this->Axis->GetPointData()->GetScalars();
  double tuple[2];
  scalars->GetTuple(0, tuple);
  double radius = tuple[0];

   vtkSphere *sphere1 = vtkSphere::New();
   sphere1->SetCenter( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
   sphere1->SetRadius( radius + 1);

   vtkImplicitBoolean *disk1 = vtkImplicitBoolean::New();
   disk1->SetOperationTypeToIntersection();
   disk1->AddFunction(sphere1);
   disk1->AddFunction(plane1);

  //ax->getControlPoint(numPts - 1, p2, p2+3);
  points->GetPoint( numPts - 1, p2);
  //ax->getControlPoint(numPts - 2, p1, p1+3); //assume it exists
  points->GetPoint( numPts - 2, p1);
  normal[0] = p1[0] - p2[0];
  normal[1] = p1[1] - p2[1];
  normal[2] = p1[2] - p2[2];
  vtkMath::Normalize( normal );

  vtkPlane *plane2 = vtkPlane::New();
  plane2->SetOrigin( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
  plane2->SetNormal( normal[0], normal[1], normal[2]);

  vtkSphere *sphere2 = vtkSphere::New();
  sphere2->SetCenter( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
  sphere2->SetRadius( radius + 1 );

  vtkImplicitBoolean *disk2 = vtkImplicitBoolean::New();
  disk2->SetOperationTypeToIntersection();
  disk2->AddFunction(sphere2);
  disk2->AddFunction(plane2);

  vtkImplicitBoolean *disk = vtkImplicitBoolean::New();
  disk->SetOperationTypeToUnion();
  disk->AddFunction(disk1);
  disk->AddFunction(disk2);

  vtkImplicitFunctionToImageStencil *functionToStencil = vtkImplicitFunctionToImageStencil::New ();
  functionToStencil->SetInput ( disk );

  vtkImageStencil *stencil = vtkImageStencil::New();
  stencil->SetInput( this->GetInput() );
  stencil->SetStencil( functionToStencil->GetOutput() );
  stencil->ReverseStencilOn();
  stencil->SetBackgroundValue (0);
  //don't forget to call ->Update before using an iterator !!!
  stencil->Update();
  
  this->ClipImageData->ShallowCopy( stencil->GetOutput() );
  
  //Delete stuff:
  plane1                        -> Delete();
  sphere1                       -> Delete();
  disk1                         -> Delete();
  plane2                        -> Delete();
  sphere2                       -> Delete();
  disk2                         -> Delete();
  disk                          -> Delete();
  functionToStencil     -> Delete();
  stencil                       -> Delete();
}
//----------------------------------------------------------------------------
void vtkImagePolyDataSeedConnectivity::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  if ( this->Axis )
    {
    os << indent << "Axis of " << this->Axis->GetNumberOfPoints()
       << "points defined\n";
    }
  else
    {
    os << indent << "Axis not defined\n";
    }
}