[Insight-users] Active-contour_RG

yasser salman yass71 at yahoo.com
Mon May 23 16:36:06 EDT 2005


Hi All, 
The attachment files describes the   using
geodesicactive contour, Connectedthreshold ,
confidence connected to extract tumor in 3D and
visualize the output using VTK ,  i'm using a data
from IBSR "Another subject with tumor, multiple scans
over time series 32i,45i,....", i have the following
questions about the output: 
I think suitable filter is geodesicactive contour not
the other two filters, I try all of them at different
seed points and geodesicAC pass while the other fails,
For example for series 47i at seed 160, 135, 21,
connected at threshold 70-100, and confidence f=2 
both filter fails , while active at beta=8, alpha =
-.3 was pass, is it convenient ? Or I have some
mistakes?
2- At series 32i I'm using beta=2 (try and error) and
using beta = 8 at 45, 47,… try and error, what is the
visual criteria  to decide the beta or alpha values
even high or low ?
3- if i use variable  "Vol_cm" to convert  no of
pixels to cm3 the type should be  unsigned long  but i
need it to be float how can i did it type casting or
what?
4- when I changed   "sizeOfObjects =
rcfilter->GetSizeOfObjectsInPixels();" TO 
sizeOfObjects = rcfilter->GetSizeOfObjectsInPixels();
"D:\InsightToolkit-1.6.0\Code\BasicFilters\
itkRelabelComponentImageFilter.h (129): error C2440:
'return' : cannot convert from 'const class
std::vector<float,class std::allocator<float> >' to
'const class std::vector<unsigned long,class
std::allocator<unsigned long> > &'Reason: cannot
convert from 'const class std::vector<float,class
std::allocator<float> >' to 'const class
std::vector<unsigned long,class
std::allocator<unsigned long> >'
 No constructor could take the source type, or
constructor overload resolution was ambiguous 
D:\InsightToolkit-1.6.0\Code\BasicFilters\itkRelabelComponentImageFilter.h(129)
: while compiling class-template member function
'const class std::vector<unsigned long,class
std::allocator<unsigned long> > &__thiscall
itk::RelabelComponentIm
ageFilter<class itk::Image<unsigned short,3>,class
itk::Image<unsigned short,3>
>::GetSizeOfObjectsInPhysicalUnits(void) const'
Error executing cl.exe."
5- The active contour filter takes long time for
calculation "about 2.5 min. per volume at Dual
processor work station with 2G ram.
6- The visualization output of the active contour
differs from the output from connected and confidence
i try it in many series specially when using surface
rendering instead of volume rendering. 
I think I have many mistakes and any help will be
appreciated. 
Regard's
Yasser Salman.



		
__________________________________ 
Do you Yahoo!? 
Yahoo! Small Business - Try our new Resources site
http://smallbusiness.yahoo.com/resources/
-------------- next part --------------
#include "itkImageFileReader.h"
#include "itkImage.h"

#include "itkCastImageFilter.h"
#include "itkConfidenceConnectedImageFilter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkConnectedThresholdImageFilter.h"
#include <itkImageMomentsCalculator.h> 
#include "itkRescaleIntensityImageFilter.h"

#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"

#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkContourFilter.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkImageActor.h"
#include "vtkImageShiftScale.h"
#include "vtkMarchingCubes.h"

#include "vtkProperty.h"
#include "vtkPolyDataNormals.h"
#include "vtkContourFilter.h"
#include "vtkImageViewer.h"
#include "vtkImageReader2.h"
#include "vtkVolumeRayCastFunction.h"
#include "vtkVolumeRayCastCompositeFunction.h"
#include "vtkVolumeRayCastMapper.h"
#include "vtkPiecewiseFunction.h"
#include "vtkVolumeProperty.h" 
#include "vtkColorTransferFunction.h"
#include "vtkImageViewer2.h"
#include "vtkPolyData.h"
#include "vtkImageData.h"
#include "vtkDICOMImageReader.h"
#include "vtkImageFlip.h"
#include "vtkContourFilter.h"
#include "vtkPolyDataMapper.h"
#include "itkDICOMImageIO2.h"
#include "itkImageSeriesReader.h"
#include "itkDICOMSeriesFileNames.h"
#include "vtkWin32OpenGLRenderWindow.h"
#include "vtkWin32RenderWindowInteractor.h"
#include "vtkCamera.h"
#include "vtkImageGaussianSmooth.h"
#include "vtkWindowToImageFilter.h"
#include "vtkBMPWriter.h"
#include "vtkImageWriter.h"
#include "vtkStructuredPointsWriter.h"
#include "vtkPolyDataWriter.h"
#include "vtkMCubesWriter.h"
#include "itkNumericSeriesFileNames.h"
#include "itkAnalyzeImageIO.h"
#include "itkMetaImageIO.h"
#include "itkRelabelComponentImageFilter.h"

#include "itkImageFileWriter.h"




#include "itkImageToVTKImageFilter.h"
#include "itkVTKImageToImageFilter.h"

int main( int argc, char ** argv )
{


  typedef unsigned short   InputPixelType;
  typedef float        InternalPixelType;
  typedef unsigned short   SegmentedPixelType;

  typedef itk::Image< InputPixelType, 3 >         InputImageType;
  typedef itk::Image< InternalPixelType, 3>       InternalImageType;
  typedef itk::Image< SegmentedPixelType, 3 >     SegmentedImageType;
  typedef itk::ImageMomentsCalculator<SegmentedImageType> ImageMomentsType;
typedef unsigned char OutputPixelType;
  typedef itk::Image< OutputPixelType, 3 > OutputImageType;
  typedef   itk::CastImageFilter< 
                 InputImageType, 
                 InternalImageType >     CastImageFilterType;

  typedef   itk::CurvatureFlowImageFilter< 
                 InternalImageType, 
                 InternalImageType >     CurvatureFlowImageFilterType;

  typedef   itk::ConfidenceConnectedImageFilter< 
                         InternalImageType, 
                         SegmentedImageType >     
ConfidenceConnectedImageFilterType;

  typedef itk::ConnectedThresholdImageFilter< 
			InternalImageType, SegmentedImageType > ConnectedFilterType;
  
  typedef itk::ImageToVTKImageFilter< SegmentedImageType >  
ITK2VTKConnectorFilterType;

  typedef itk::VTKImageToImageFilter< InputImageType     >  
VTK2ITKConnectorFilterType;


typedef itk::BinaryThresholdImageFilter< 
                        InternalImageType, 
                        SegmentedImageType /*char*/   >    ThresholdingFilterType; 

typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter< 
                               InternalImageType, 
                               InternalImageType >  GradientFilterType;
typedef  itk::FastMarchingImageFilter< 
                              InternalImageType, 
                              InternalImageType >    FastMarchingFilterType;
 typedef   itk::SigmoidImageFilter<                               
                               InternalImageType, 
                               InternalImageType >  SigmoidFilterType;
typedef  itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType, 
                InternalImageType >    GeodesicActiveContourFilterType;

ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
//SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();
SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
GeodesicActiveContourFilterType::Pointer geodesicActiveContour = GeodesicActiveContourFilterType::New();


  /*  vtkImageReader2 *Reader=vtkImageReader2::New(); 
	vtkImageViewer2 *viewer = vtkImageViewer2::New();     
	Reader->SetFilePrefix("D:/ImageData/isbr/126_26.img/126_26_");
	Reader->SetFilePattern("%s%.d.img");
	Reader->SetDataScalarTypeToUnsignedShort();
  Reader->SetDataByteOrderToBigEndian();
  Reader->SetDataExtent(0,255,0,255,1,28);//
   Reader->SetDataSpacing(1,1,1.7);
    Reader->Update();										  

	viewer->SetInput(Reader->GetOutput());				  
int VolData_Images = viewer->GetWholeZMax ();			  


for (int i=1;i<VolData_Images;i++)
{
	
	viewer->SetZSlice(i);		  
    viewer->SetSize(512,512);
	viewer->SetPosition(500,0);
	viewer->SetColorWindow(512);
	viewer->SetColorLevel(256);
	viewer->Render();

// TODO: Add your command handler code here
}*/	

  //----------------------------------
  //           ITK Pipeline
  //----------------------------------

 // READ ANALAYZE IMAGES..... 
  typedef itk::ImageFileReader< InputImageType >     ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
  
  typedef  itk::ImageFileWriter<  SegmentedImageType  > WriterType;

  WriterType::Pointer writer = WriterType::New();
  reader->SetFileName(  argv[1]  );
  try
    {
	 
    reader->Update(); 
    }
  catch (itk::ExceptionObject &ex)
    {
    std::cout << ex << std::endl;
    return EXIT_FAILURE;
    }
  

  CastImageFilterType::Pointer cast = CastImageFilterType::New();

  CurvatureFlowImageFilterType::Pointer smoothing = 
CurvatureFlowImageFilterType::New();

  ConfidenceConnectedImageFilterType::Pointer confidence = 
ConfidenceConnectedImageFilterType::New();

  ConnectedFilterType::Pointer connectedThreshold = 
	  ConnectedFilterType::New();

  VTK2ITKConnectorFilterType::Pointer VTK2ITKconnector = 
VTK2ITKConnectorFilterType::New();
  ImageMomentsType::Pointer calculator=ImageMomentsType::New();
  

  cast->SetInput( reader->GetOutput() );

  smoothing->SetInput(  cast->GetOutput() );
  smoothing->SetTimeStep( 0.0625 );
  smoothing->SetNumberOfIterations( 2 );
   thresholder->SetLowerThreshold( -1000.0 );
  thresholder->SetUpperThreshold(     0.0 );

  thresholder->SetOutsideValue(  0  );
  thresholder->SetInsideValue(  255 );
  sigmoid->SetOutputMinimum(  0.0  );
  sigmoid->SetOutputMaximum(  1.0  );
    const double propagationScaling = atof( argv[10] );

  gradientMagnitude->SetInput( smoothing->GetOutput() );
  sigmoid->SetInput( gradientMagnitude->GetOutput() );

geodesicActiveContour->SetInput(  fastMarching->GetOutput() );
geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
thresholder->SetInput( geodesicActiveContour->GetOutput() );
const double sigma = atof( argv[6] );
gradientMagnitude->SetSigma(  sigma  );
const double alpha = atof( argv[7] );
const double beta  =  atof( argv[8] );
sigmoid->SetAlpha( alpha );
sigmoid->SetBeta(  beta  );
 typedef FastMarchingFilterType::NodeContainer  NodeContainer;
 typedef FastMarchingFilterType::NodeType       NodeType;
  NodeContainer::Pointer seeds = NodeContainer::New();
  InternalImageType::IndexType  seedPosition;

  seedPosition[0] = atoi( argv[3] ); //135
  seedPosition[1] = atoi( argv[4] );
  seedPosition[2]=atoi( argv[5] );
  const double initialDistance = atof( argv[9] );
  NodeType node;
  const double seedValue = - initialDistance;
  node.SetValue( seedValue );
  node.SetIndex( seedPosition );
  seeds->Initialize();
  seeds->InsertElement( 0, node );
   std::cout << "Using seed = " << seedPosition << std::endl;
  fastMarching->SetTrialPoints(  seeds  );
  fastMarching->SetSpeedConstant( 1.0 );
  fastMarching->SetOutputSize(reader->GetOutput()->GetBufferedRegion().GetSize() );


  std::cout << "Using seed = " << seeds << std::endl;
  

  typedef itk::RelabelComponentImageFilter<SegmentedImageType,SegmentedImageType> RelabelImageFilter;
RelabelImageFilter::Pointer rcfilter= RelabelImageFilter::New();
std::vector<unsigned long>::const_iterator it;
rcfilter->SetInput(thresholder->GetOutput());  
rcfilter->Update();
unsigned long nObjects = rcfilter->GetNumberOfObjects();
	const std::vector<unsigned long> sizeOfObjects = rcfilter->GetSizeOfObjectsInPixels();
it =  sizeOfObjects.begin();
  cout<<"No. of pixels ="<<*it<<"\n";

typedef itk::RescaleIntensityImageFilter<SegmentedImageType,SegmentedImageType>   RescaleType;

RescaleType:: Pointer rescale= RescaleType::New();
rescale->SetInput(rcfilter->GetOutput()); // Case "confidence" replace "connectedThreshold"  
rescale->SetOutputMinimum(0);
rescale->SetOutputMaximum(255);
rescale->Update();
   writer->SetInput( rescale->GetOutput() );
  writer->SetFileName( argv[2] );
writer->Update();
  //

ITK2VTKConnectorFilterType::Pointer ITK2VTKconnector = 
ITK2VTKConnectorFilterType::New();

  ITK2VTKconnector->SetInput(rescale->GetOutput());
  ITK2VTKconnector->GetImporter()->SetDataScalarTypeToUnsignedChar();
  
  //----------------------------------
  //           VTK Pipeline
  //----------------------------------

  vtkRenderer * renderer = vtkRenderer::New();
  vtkWin32OpenGLRenderWindow * renderWindow = vtkWin32OpenGLRenderWindow::New();
  vtkWin32RenderWindowInteractor * renderWindowInteractor = vtkWin32RenderWindowInteractor::New();
vtkImageActor* imageactor = vtkImageActor::New();
vtkActor *actor = vtkActor::New();
  renderWindow->AddRenderer( renderer );
  renderWindow->SetInteractor( renderWindowInteractor );
 vtkMarchingCubes *mc = vtkMarchingCubes::New();
//viewer->SetInput( ITK2VTKconnector->GetOutput() );
 vtkImageGaussianSmooth  *gaussian=vtkImageGaussianSmooth::New();
 gaussian->SetStandardDeviations (1, 1, 1);
     gaussian-> SetRadiusFactors (1.0, 1.0, 1.0);
gaussian ->SetInput(ITK2VTKconnector->GetOutput());

/*	mc->SetInput( gaussian->GetOutput() );  
	mc->SetValue(0,127);

	vtkPolyDataMapper* mcmap=vtkPolyDataMapper::New();
	mcmap->SetInput(mc->GetOutput());
	mcmap->ScalarVisibilityOff();
	vtkActor* mcactor=vtkActor::New();
	mcactor->SetMapper(mcmap); 
	renderer->AddActor(mcactor);//*/
	
	vtkPiecewiseFunction *oTFun =vtkPiecewiseFunction::New();
	oTFun->AddSegment(80,0.0,255,1);
	vtkPiecewiseFunction *gTFun =vtkPiecewiseFunction::New();
	gTFun->AddSegment(0,1.0,255,1.0);

	vtkVolumeProperty *volProperty = vtkVolumeProperty::New();
	volProperty->SetColor(gTFun);

	volProperty->SetScalarOpacity(oTFun);
	volProperty->SetInterpolationTypeToLinear();
	volProperty->ShadeOn();
	vtkVolumeRayCastCompositeFunction *compositeFunction= vtkVolumeRayCastCompositeFunction::New();
	vtkVolumeRayCastMapper *volMapper = vtkVolumeRayCastMapper::New();
    volMapper->SetInput(ITK2VTKconnector->GetOutput()); 
	volMapper->SetVolumeRayCastFunction(compositeFunction);

	vtkVolume *vol=vtkVolume::New();
	vol->SetMapper(volMapper);  //in original SetVolumeMapper
	vol->SetProperty(volProperty); // in original SetVolumeProperty
	renderer->AddVolume(vol);
	renderer->SetBackground(1,1,1);


	vtkCamera *aCamera = vtkCamera::New();
     aCamera->SetViewUp (0,  1, 0);
    aCamera->SetPosition (0, 0, 1);
    aCamera->SetFocalPoint (1, 0, 1);
    aCamera->ComputeViewPlaneNormal();
	 renderer->SetActiveCamera(aCamera);
	 renderer->ResetCamera ();
	 aCamera->Azimuth(0);


	 



	 renderWindow->SetSize(320, 320);
	renderWindowInteractor->SetDesiredUpdateRate(3.0);
	
	renderer->SetBackground( 0.1, 0.2, 0.4 );
  renderWindow->Render();
 renderWindowInteractor->Start();
 
  return 0;

}


-------------- next part --------------
#include "itkImageFileReader.h"
#include "itkImage.h"

#include "itkCastImageFilter.h"
#include "itkConfidenceConnectedImageFilter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkConnectedThresholdImageFilter.h"
#include <itkImageMomentsCalculator.h> 
#include "itkRescaleIntensityImageFilter.h"


#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkContourFilter.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkImageActor.h"
#include "vtkImageShiftScale.h"
#include "vtkMarchingCubes.h"

#include "vtkProperty.h"
#include "vtkPolyDataNormals.h"
#include "vtkContourFilter.h"
#include "vtkImageViewer.h"
#include "vtkImageReader2.h"
#include "vtkVolumeRayCastFunction.h"
#include "vtkVolumeRayCastCompositeFunction.h"
#include "vtkVolumeRayCastMapper.h"
#include "vtkPiecewiseFunction.h"
#include "vtkVolumeProperty.h" 
#include "vtkColorTransferFunction.h"
#include "vtkImageViewer2.h"
#include "vtkPolyData.h"
#include "vtkImageData.h"
#include "vtkDICOMImageReader.h"
#include "vtkImageFlip.h"
#include "vtkContourFilter.h"
#include "vtkPolyDataMapper.h"
#include "itkDICOMImageIO2.h"
#include "itkImageSeriesReader.h"
#include "itkDICOMSeriesFileNames.h"
#include "vtkWin32OpenGLRenderWindow.h"
#include "vtkWin32RenderWindowInteractor.h"
#include "vtkCamera.h"
#include "vtkImageGaussianSmooth.h"
#include "vtkWindowToImageFilter.h"
#include "vtkBMPWriter.h"
#include "vtkImageWriter.h"
#include "vtkStructuredPointsWriter.h"
#include "vtkPolyDataWriter.h"
#include "vtkMCubesWriter.h"
#include "itkNumericSeriesFileNames.h"
#include "itkAnalyzeImageIO.h"
#include "itkMetaImageIO.h"
#include "itkRelabelComponentImageFilter.h"
#include "vtkLookupTable.h"
#include "vtkImageMapToColors.h"
#include "vtkImageBlend.h"
#include "vtkImageCast.h"







#include "itkImageToVTKImageFilter.h"
#include "itkVTKImageToImageFilter.h"

int main( int argc, char ** argv )
{


  typedef unsigned short   InputPixelType;
  typedef float        InternalPixelType;
  typedef unsigned short   SegmentedPixelType;

  typedef itk::Image< InputPixelType, 3 >         InputImageType;
  typedef itk::Image< InternalPixelType, 3>       InternalImageType;
  typedef itk::Image< SegmentedPixelType, 3 >     SegmentedImageType;
  typedef itk::ImageMomentsCalculator<SegmentedImageType> ImageMomentsType;

  typedef   itk::CastImageFilter< 
                 InputImageType, 
                 InternalImageType >     CastImageFilterType;

  typedef   itk::CurvatureFlowImageFilter< 
                 InternalImageType, 
                 InternalImageType >     CurvatureFlowImageFilterType;

  typedef   itk::ConfidenceConnectedImageFilter< 
                         InternalImageType, 
                         SegmentedImageType >     
ConfidenceConnectedImageFilterType;

  typedef itk::ConnectedThresholdImageFilter< 
			InternalImageType, SegmentedImageType > ConnectedFilterType;
  
  typedef itk::ImageToVTKImageFilter< SegmentedImageType >  
ITK2VTKConnectorFilterType;

  typedef itk::VTKImageToImageFilter< InputImageType     >  
VTK2ITKConnectorFilterType;




    vtkImageReader2 *Reader=vtkImageReader2::New(); 
	vtkImageViewer2 *viewer = vtkImageViewer2::New();     
	Reader->SetFilePrefix("D:/ImageData/img2/32i/536_32_");
	Reader->SetFilePattern("%s%.d.img");
	Reader->SetDataScalarTypeToUnsignedShort();
  Reader->SetDataByteOrderToBigEndian();
  Reader->SetDataExtent(0,255,0,255,1,60);//
   Reader->SetDataSpacing(.9,.9,2);
    Reader->Update();										  

	viewer->SetInput(Reader->GetOutput());				  
int VolData_Images = viewer->GetWholeZMax ();			  


for (int i=1;i<VolData_Images;i++)
{
	
	viewer->SetZSlice(i);		  
    viewer->SetSize(512,512);
	viewer->SetPosition(500,0);
	viewer->SetColorWindow(512);
	viewer->SetColorLevel(256);
	viewer->Render();

// TODO: Add your command handler code here
}	

  //----------------------------------
  //           ITK Pipeline
  //----------------------------------

 // READ ANALAYZE IMAGES..... 
  typedef itk::ImageFileReader< InputImageType >     ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(  argv[1]   );
  try
    {
	 
    reader->Update(); 
    }
  catch (itk::ExceptionObject &ex)
    {
    std::cout << ex << std::endl;
    return EXIT_FAILURE;
    }
  

  CastImageFilterType::Pointer cast = CastImageFilterType::New();

  CurvatureFlowImageFilterType::Pointer smoothing = 
CurvatureFlowImageFilterType::New();

  ConfidenceConnectedImageFilterType::Pointer confidence = 
ConfidenceConnectedImageFilterType::New();

  ConnectedFilterType::Pointer connectedThreshold = 
	  ConnectedFilterType::New();

  VTK2ITKConnectorFilterType::Pointer VTK2ITKconnector = 
VTK2ITKConnectorFilterType::New();
  ImageMomentsType::Pointer calculator=ImageMomentsType::New();
  

  cast->SetInput( reader->GetOutput() );

  smoothing->SetInput(  cast->GetOutput() );
  confidence->SetInput( smoothing->GetOutput() );

  smoothing->SetTimeStep( 0.0625 );
  smoothing->SetNumberOfIterations( 2 );
    connectedThreshold->SetInput( smoothing->GetOutput() );
	 const InternalPixelType lowerThreshold = atoi( argv[2] );
  const InternalPixelType upperThreshold = atoi( argv[3] );
  connectedThreshold->SetLower(  lowerThreshold  );
  connectedThreshold->SetUpper(  upperThreshold  );
  connectedThreshold->SetReplaceValue( 255 );


  typedef ConnectedFilterType::IndexType IndexType;
  IndexType seed;
  seed[0] = atoi( argv[4] );
  seed[1] = atoi( argv[5] );
  seed[2] = atoi( argv[6] );

  std::cout << "Using seed = " << seed << std::endl;
  
  connectedThreshold->SetSeed( seed );
  
  confidence->SetSeed( seed );
  confidence->SetMultiplier( 1);
  confidence->SetNumberOfIterations(1);
  confidence->SetInitialNeighborhoodRadius(2);
  confidence->SetReplaceValue( 255 );

  typedef itk::RelabelComponentImageFilter<SegmentedImageType,SegmentedImageType> RelabelImageFilter;
RelabelImageFilter::Pointer rcfilter= RelabelImageFilter::New();
std::vector<unsigned long>::const_iterator it;
rcfilter->SetInput(connectedThreshold->GetOutput());  
rcfilter->Update();
unsigned long nObjects = rcfilter->GetNumberOfObjects();
	const std::vector<unsigned long> sizeOfObjects = rcfilter->GetSizeOfObjectsInPixels();
it =  sizeOfObjects.begin();
  cout<<"No. of pixels ="<<*it<<"\n";
unsigned long Vol_cm;
  Vol_cm=(.9*.9*2*(*it))/1000; 
  cout<<"Volume in cm ="<<Vol_cm<<"\n";


typedef itk::RescaleIntensityImageFilter<SegmentedImageType,SegmentedImageType>   RescaleType;

RescaleType:: Pointer rescale= RescaleType::New();
rescale->SetInput(rcfilter->GetOutput()); // Case "confidence" replace "connectedThreshold"  
rescale->SetOutputMinimum(0);
rescale->SetOutputMaximum(255);
rescale->Update();


  //

ITK2VTKConnectorFilterType::Pointer ITK2VTKconnector = 
ITK2VTKConnectorFilterType::New();

  ITK2VTKconnector->SetInput(rescale->GetOutput());
  ITK2VTKconnector->GetImporter()->SetDataScalarTypeToUnsignedChar();
  
  //----------------------------------
  //           VTK Pipeline
  //----------------------------------

  vtkRenderer * renderer = vtkRenderer::New();
  vtkWin32OpenGLRenderWindow * renderWindow = vtkWin32OpenGLRenderWindow::New();
  vtkWin32RenderWindowInteractor * renderWindowInteractor = vtkWin32RenderWindowInteractor::New();
vtkImageActor* imageactor = vtkImageActor::New();
vtkActor *actor = vtkActor::New();
  renderWindow->AddRenderer( renderer );
  renderWindow->SetInteractor( renderWindowInteractor );
 vtkMarchingCubes *mc = vtkMarchingCubes::New();
//viewer->SetInput( ITK2VTKconnector->GetOutput() );
 vtkImageGaussianSmooth  *gaussian=vtkImageGaussianSmooth::New();
 gaussian->SetStandardDeviations (1, 1, 1);
     gaussian-> SetRadiusFactors (1.0, 1.0, 1.0);
gaussian ->SetInput(ITK2VTKconnector->GetOutput());

/*	mc->SetInput( gaussian->GetOutput() );  
	mc->SetValue(0,127);

	vtkPolyDataMapper* mcmap=vtkPolyDataMapper::New();
	mcmap->SetInput(mc->GetOutput());
	mcmap->ScalarVisibilityOff();
	vtkActor* mcactor=vtkActor::New();
	mcactor->SetMapper(mcmap); 
	renderer->AddActor(mcactor);//*/
	

	vtkPiecewiseFunction *oTFun =vtkPiecewiseFunction::New();
	oTFun->AddSegment(80,0.0,255,1);
	vtkPiecewiseFunction *gTFun =vtkPiecewiseFunction::New();
	gTFun->AddSegment(0,1.0,255,1.0);

	vtkVolumeProperty *volProperty = vtkVolumeProperty::New();
	volProperty->SetColor(gTFun);
	vtkLookupTable *hueLut0 = vtkLookupTable::New();
 hueLut0->SetTableRange (0,2000);
 hueLut0->SetSaturationRange (0,1);
    hueLut0->SetValueRange (0, 5);

	hueLut0->SetHueRange (.6,.6);
	
	vtkLookupTable *hueLut1 = vtkLookupTable::New();
	hueLut1->SetTableRange (0,2000);
    hueLut1->SetSaturationRange (1,1);
    hueLut1->SetValueRange (0, 5);
	hueLut1->SetHueRange (1, 1);

  vtkImageMapToColors *alpha0 = vtkImageMapToColors::New();
    alpha0->SetInput(Reader->GetOutput());
    alpha0->SetLookupTable(hueLut0);
  vtkImageMapToColors *alpha1 = vtkImageMapToColors::New();
    alpha1->SetInput(gaussian->GetOutput());
   alpha1->SetLookupTable(hueLut1);
   
   vtkImageCast *cast0=vtkImageCast::New();
   cast0->SetInput(alpha0->GetOutput());
   cast0->SetOutputScalarTypeToUnsignedShort();
   vtkImageCast *cast1=vtkImageCast::New();
   cast1->SetInput(alpha1->GetOutput());
   cast1->SetOutputScalarTypeToUnsignedShort();
   
  vtkImageBlend *blend=vtkImageBlend::New();
    blend-> SetInput(0,Reader->GetOutput());
    blend-> SetInput(1,gaussian->GetOutput());// end of 1
   blend->SetOpacity(0,.9);
	blend->SetOpacity(1,.7);
	volProperty->SetScalarOpacity(oTFun);
	volProperty->SetInterpolationTypeToLinear();
	volProperty->ShadeOn();
	vtkVolumeRayCastCompositeFunction *compositeFunction= vtkVolumeRayCastCompositeFunction::New();
	vtkVolumeRayCastMapper *volMapper = vtkVolumeRayCastMapper::New();
    volMapper->SetInput(blend->GetOutput()); //if tumor only gaussian
	volMapper->SetVolumeRayCastFunction(compositeFunction);

	vtkVolume *vol=vtkVolume::New();
	vol->SetMapper(volMapper);  //in original SetVolumeMapper
	vol->SetProperty(volProperty); // in original SetVolumeProperty
	

	
	
	renderer->AddVolume(vol);
	renderer->SetBackground(1,1,1);


	vtkCamera *aCamera = vtkCamera::New();
     aCamera->SetViewUp (0,  1, 0);
    aCamera->SetPosition (0, 0, 1);
    aCamera->SetFocalPoint (1, 0, 1);
    aCamera->ComputeViewPlaneNormal();
	 renderer->SetActiveCamera(aCamera);
	 renderer->ResetCamera ();
	 aCamera->Azimuth(0);


	 



	 renderWindow->SetSize(320, 320);
	renderWindowInteractor->SetDesiredUpdateRate(3.0);
	
	renderer->SetBackground( 0.1, 0.2, 0.4 );
  renderWindow->Render();
 renderWindowInteractor->Start();
 
  return 0;
}




More information about the Insight-users mailing list