[Insight-users] Calculate exact surface

[THOMAS Diego]

Hi luis,

Thanks for your help,

I have some problems to create the simplex mesh,
I use the itk TriangleMeshToSimplexMesh Filter after the BinaryMask3DMeshSource.
I have no error message and no segmentation fault during the execution of the algorithm,
but when I update the TriangleMeshToSimplexMesh Filter the calculation doesn't end.
When I create the mesh with itk RegularSphereMeshSource instead of BinaryMask3DMeshSource, it works,
but when I use BinaryMask3DMeshSource I can't transform my mesh.

What am I missing?

Thanks in advance,
Diego

Here is my code:
#include "SurfaceTool.h"
#include "ImageTool.h"
 
#include <math.h>
#include <iostream>

#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkAntiAliasBinaryImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"
#include "itkBinaryMask3DMeshSource.h"
#include "itkMesh.h"
#include "itkSpatialObjectWriter.h"
#include "itkMeshSpatialObject.h"
#include "itkDefaultDynamicMeshTraits.h"
#include "itkSimplexMeshVolumeCalculator.h"
#include "itkTriangleMeshToSimplexMeshFilter.h"
#include "itkSimplexMesh.h"
#include "itkTriangleCell.h"
#include "itkImage.h"

typedef itk::Image<unsigned short, 3> Image3D;
typedef itk::Image<float, 3> ImageReal3D;
typedef itk::CastImageFilter< Image3D, ImageReal3D> CastToRealFilterType;
typedef itk::RescaleIntensityImageFilter<ImageReal3D, Image3D > RescaleFilter;
typedef itk::AntiAliasBinaryImageFilter<ImageReal3D, ImageReal3D> AntiAliasFilterType;
typedef itk::ImageRegionConstIterator<Image3D> ConstIteratorType;
typedef itk::ImageRegionIterator<Image3D> IteratorType;
typedef itk::DefaultDynamicMeshTraits< double,3,3, double, double, double> MeshTrait;
typedef itk::Mesh<double,3,MeshTrait> MeshType;
typedef itk::SimplexMesh<double, 3, MeshTrait> SimplexMeshType;
typedef itk::BinaryMask3DMeshSource< Image3D, MeshType > MeshSourceType;
typedef itk::TriangleMeshToSimplexMeshFilter<MeshType, SimplexMeshType> TriangleToSimplexFilter;
typedef itk::SimplexMeshVolumeCalculator<SimplexMeshType> CalculatorType;



Image3D::Pointer SurfaceTool::GetSurface( Image3D::Pointer data)
{
	//Instantiation of Caster and rescale filter
	CastToRealFilterType::Pointer toReal = CastToRealFilterType::New();
	RescaleFilter::Pointer rescale = RescaleFilter::New();
	rescale->SetOutputMinimum(   0 );
  	rescale->SetOutputMaximum( 255 );    
	
	AntiAliasFilterType::Pointer antiAliasFilter = AntiAliasFilterType::New();    
	
	toReal->SetInput( data );  
	antiAliasFilter->SetInput( toReal->GetOutput() );
	
	antiAliasFilter->SetMaximumRMSError( 0.01 );
	antiAliasFilter->SetNumberOfIterations( 50 );
	antiAliasFilter->SetNumberOfLayers( 2 ); 
	//antiAliasFilter->SetIsoSurfaceValue (0.0);
	
	rescale->SetInput( antiAliasFilter->GetOutput() );
	try 
    {
    	rescale->Update();
    	std::cout << "Completed in " << antiAliasFilter->GetNumberOfIterations() << std::endl;	//->Completed in 50
    	std::cout << "Iso surface value " << antiAliasFilter->GetIsoSurfaceValue () << std::endl;	//->Iso surface value 0.5
    	std::cout << "Last RMS change value was: " << antiAliasFilter->GetRMSChange() << std::endl;	//Last RMS change value was: 0.00946563 
    }
  	catch( itk::ExceptionObject & err ) 
    { 
    	std::cout << "ExceptionObject caught !" << std::endl; 
    	std::cout << err << std::endl; 
    	return data;
    } 
    
    float min = rescale->GetInputMinimum();
    float max = rescale->GetInputMaximum();
    
    double isovalue = (double)( antiAliasFilter->GetIsoSurfaceValue () - min )*(255.0)/(max-min);
    
     
    Image3D::Pointer res = rescale->GetOutput();
	
  	IteratorType outputItTmp(res, res->GetLargestPossibleRegion());
  	for (outputItTmp.GoToBegin(); !outputItTmp.IsAtEnd(); ++outputItTmp)
    {
     	if (outputItTmp.Get() >= 127)
     	{
    		outputItTmp.Set(1);
     	} else {
			outputItTmp.Set(0);
		}
    }
    
    MeshSourceType::Pointer meshSource = MeshSourceType::New();
    meshSource->SetObjectValue( 1 );
    meshSource->SetInput( res );
    try
	{
		meshSource->Update();
	}
	catch( itk::ExceptionObject & exp )
	{
		std::cerr << "Exception thrown during Update() " << std::endl;
		std::cerr << exp << std::endl;
		return data;
	}
	
	std::cout << "Nodes = " << meshSource->GetNumberOfNodes() << std::endl;	//->Nodes = 2003
	std::cout << "Cells = " << meshSource->GetNumberOfCells() << std::endl; //->Cells = 3904
	
	TriangleToSimplexFilter::Pointer Transform = TriangleToSimplexFilter::New();
	
	Transform->SetInput( meshSource->GetOutput() );
	Transform->Update();
	
	CalculatorType::Pointer Calculator = CalculatorType::New();
	
	Calculator->SetSimplexMesh(Transform->GetOutput());
	Calculator->Compute();
	
	std::cout << "Volume = " << Calculator->GetVolume() << std::endl;
	std::cout << "Area = " << Calculator->GetArea() << std::endl;
   
	return data;
} 

-----Message d'origine-----
De : Luis Ibanez [mailto:luis.ibanez at kitware.com]
Envoyé : mardi 28 août 2007 19:03
À : THOMAS Diego
Cc : Insight Users
Objet : Re: [Insight-users] Calculate exact surface



Hi Thomas,

You may want to use the
http://www.itk.org/Insight/Doxygen/html/classitk_1_1SimplexMeshVolumeCalculator.html

The code is in

     Insight/Code/Algorithms/
        itkSimplexMeshVolumeCalculator.h*
        itkSimplexMeshVolumeCalculator.txx


Take a look at the GetArea() method.

You may need to use also the filter:

http://www.itk.org/Insight/Doxygen/html/classitk_1_1TriangleMeshToSimplexMeshFilter.html

in order to get the intermediate simples mesh.



   Regards,


      Luis


-------------------
THOMAS Diego wrote:
> Hi,
> 
> I would like to calculate the exact surface of an object recorded as a 3D image with itk.
> 
> Does anybody have an idea?
> 
> I though about using the itk antialiasing filter, then the itk extractsurface filter and sum the areas of all the triangles of the mesh.
> Does this do the right thing?
> 
> thanks in advance,
> 
> regards
> 
> Diego.
> __________________________
> 
> Ce message (et toutes ses pièces jointes éventuelles) est confidentiel et établi à l'intention exclusive de ses destinataires. Toute utilisation de ce message non conforme à sa destination, toute diffusion ou toute publication, totale ou partielle, est interdite, sauf autorisation expresse. L'IFP décline toute responsabilité au titre de ce message.
> 
> This message and any attachments (the message) are confidential and intended solely for the addressees. Any unauthorised use or dissemination is prohibited. IFP should not be liable for this message.
> 
> Visitez notre site Web / Visit our web site : http://www.ifp.fr
> __________________________
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
> 
__________________________

Ce message (et toutes ses pièces jointes éventuelles) est confidentiel et établi à l'intention exclusive de ses destinataires. Toute utilisation de ce message non conforme à sa destination, toute diffusion ou toute publication, totale ou partielle, est interdite, sauf autorisation expresse. L'IFP décline toute responsabilité au titre de ce message.

This message and any attachments (the message) are confidential and intended solely for the addressees. Any unauthorised use or dissemination is prohibited. IFP should not be liable for this message.

Visitez notre site Web / Visit our web site : http://www.ifp.fr
__________________________