[Insight-users] Calculate exact surface
Luis Ibanez
luis.ibanez at kitware.com
Mon Sep 3 18:48:35 EDT 2007
Hi Thomas,
This sounds like a bug in the TriangleMeshToSimplexMesh filter.
Tell us,
the surface that you are extracting with the BinaryMask3DMeshSource
is it a closed surface ?
I would suspect that the TriangleMeshToSimplexMesh
makes that assumption.
Please let us know,
Thanks
BTW: You may find interesting to look at the code of the application:
InsightApplications/
SimplexMeshDeformableSurface
This applications does exactly what you are describing.
Luis
--------------------
THOMAS Diego wrote:
> 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
> __________________________
>
More information about the Insight-users
mailing list