[Insight-users] Problems using the ithDeformableMesh3DFilter
filter...need help!
Luis Ibanez
luis.ibanez at kitware.com
Thu, 04 Mar 2004 22:36:51 -0500
Hi Julien,
Thanks a lot for providing such complete description
of your system and the problem you are facing.
About your questions:
1) The warning:
"identifier was truncated to '255' characters
in the debug information"
is an annoying known 'feature' of Visual Studio.
The symbol names generated from templated code
(in particular from normal STL) are too long
for the 255 characters limit of in this compiler.
We use to simply disable the warninng using
the pragma:
#ifdef _MSC_VER
#pragma warning ( disable : 4786 )
#endif
Note that this doesn't really solve the problem.
It just ask the compiler to stop complaining
about it, so you can see other more significant
warnings.
2) The set of messages at run time are somehow
related to (1). Given that the symbols were
truncated, when the program is crashing, there
is no way to recover some of the symbols.
The real concern is : why the program is crashing
One of the typical pathological behaviors of this
filter is to generate forces that are too large
for the scale of displacements, making the nodes
in the mesh be propulsed to far away locations.
Note that this particular filter is extremely
sensitive to the force computations. You may
want to tune it in order to generate very small
forces in order to produce small displacements
of the mesh nodes.
In particular you want to tune the parameters
of "stiffnes" and "time step" in order to
control the magnitude of the forces.
Just for the record,
You probably may want to reconsider your approach for
characterizing the displacements of the brain surface...
Are you bulding something like a finite elements model
of the brain ?
If your first goal is to obtain a segmentation of the
brain, you may get better results using Level Sets,
which are a more robuts and reliable implementation
of deformable models.
Regards,
Luis
------------------------
Julien Mercenier wrote:
> Hi everybody,
>
> Here is a little description of my system : - athlon1400MHz, 256
> DDRAM, Windows XP
> -
> ALL_BUILD - Win32Debug for ITK, VTK and all the programs I write
>
>
> I want to characterize the displacements of the brain surface in 3D.
> To do this, I use the ITK User-Guide exemple using the
> ithDeformableMesh3DFilter filter (page 379).
> But, except of using the itkBinaryMask3DMeshSource filter, I create
> the initial mesh myself.
> This mesh is the result coming from the program Isosurf giving a
> surface mesh in 3D.
> For the calculation of the force field, I use here the same method
> as in the ITK User-guide example.
> The binary image used for this calculation is a deformed brain, in
> *.img format (*.hdr as header)
> (the dimension of this image is 256*256*128 with 00 (not in the
> brain) or 01(in the brain) => this is thus
> a file of 8Mo). The pixel distances are 0.98, 0.98 and 1.56 (mm).
>
> There are several things I don't understant or I can't solve by myself.
>
> - *At compiling-time, there are 56 warnings (all the sames but for
> different variables).* Here is an example :
>
> J:\Program Files\Microsoft Visual Studio\VC98\INCLUDE\xmemory(64) :
> warning C4786:
> 'std::allocator<itk::SmartPointer<itk::CellInterfaceVisitor<double,itk::CellTraitsInfo<3,double,double,unsigned
> long,unsigned long,unsigned long,itk::Point<double,3>,
> itk::VectorContainer<unsigned long,itk::Point<double,3>
> >,std::set<unsigned long,std::less<unsigned
> long>,std::allocator<unsigned long> > > > > >' : identifier was
> truncated to '255' characters in the debug information
> J:\Program Files\Microsoft Visual
> Studio\VC98\INCLUDE\xmemory(64) : while compiling class-template member
> function 'void __thiscall
> std::allocator<itk::SmartPointer<itk::CellInterfaceVisitor<double,itk::CellTraitsInfo<3,double,double,unsigne
> d long,unsigned long,unsigned
> long,itk::Point<double,3>,itk::VectorContainer<unsigned
> long,itk::Point<double,3> >,std::set<unsigned long,std::less<unsigned
> long>,std::allocator<unsigned long> > > > > >::deallocate(void
> *,unsigned int)'
>
> - *At running-time*
> The program stops when using the ithDeformableMesh3DFilter filter :
> there is a bug !! But I don't understand why...even when reading the codes.
>
> _Then, am message apperas in a window. And, when I click on the debug
> button, a message appears in the visualC++ window :_
> message = " Unhandled exception in def_surf.exe : 0x0000005 : Access
> Violation"
> in vxl/vnl/vnl_matrix_fixed.h indicating the " //get element " function.
>
> _In the debug window appears the messages :_
> Loaded symbols for 'J:\DEF_SURF\Debug\def_surf.exe'
> Loaded 'J:\WINDOWS\system32\ntdll.dll', no matching symbolic information
> found.
> Loaded 'J:\WINDOWS\system32\kernel32.dll', no matching symbolic
> information found.
> Loaded symbols for 'J:\WINDOWS\system32\MSVCP60D.DLL'
> Loaded symbols for 'J:\WINDOWS\system32\MSVCRTD.DLL'
> Loaded 'J:\WINDOWS\system32\user32.dll', no matching symbolic
> information found.
> Loaded 'J:\WINDOWS\system32\gdi32.dll', no matching symbolic information
> found.
> Loaded 'J:\WINDOWS\system32\advapi32.dll', no matching symbolic
> information found.
> Loaded 'J:\WINDOWS\system32\rpcrt4.dll', no matching symbolic
> information found.
> Loaded 'J:\WINDOWS\system32\version.dll', no matching symbolic
> information found.
> Loaded 'J:\WINDOWS\system32\msvcrt.dll', no matching symbolic
> information found.
> Loaded 'J:\WINDOWS\system32\SHLWAPI.DLL', no matching symbolic
> information found.
> Loaded 'J:\WINDOWS\system32\apphelp.dll', no matching symbolic
> information found.
> The thread 0xBB0 has exited with code 0 (0x0).
> First-chance exception in def_surf.exe: 0xC0000005: Access Violation.
> The program 'J:\DEF_SURF\Debug\def_surf.exe' has exited with code 0 (0x0).
>
> *I put my code just here *(the creation of the mesh don't seem to pose
> any problem)
> **
> /*************************************************************
> *******active deformable surface**************************
> **************************************************************/
>
>
> /*************************************************************
> ----------------
> Several steps ==>
> ----------------
> 1) use the file to create the mesh
> 2) computing the force field (gradient)
> 3) mesh and forxce field into the DeformableMeshFilter
> ***************************************************************/
>
> #include <cstdio> //C-style I/O to read the Isosurf file
> #include <stdio.h>
>
> #include <iostream>
>
> #include "itkMesh.h"
> #include "itkDefaultStaticMeshTraits.h"
> #include "itkTriangleCell.h"
>
> #include "itkDeformableMesh3DFilter.h"
>
> #include "itkGradientRecursiveGaussianImageFilter.h"
> #include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
>
> #include "itkImage.h"
>
> #include "itkCovariantVector.h"//will contain the force field
>
> #include "itkImageFileReader.h"
>
>
>
> int main()
> {
>
> int nbnoeuds = 0;
> int nbtri = 0;
>
> double noeudsAr[3648][3]; //3648 is the number of nodes
> double noeudsdefAr[3648][3];
> int triAr[7280][3]; //7280 is the number of triangles
>
> FILE* fichier = 0;
>
> /*opening of the file containing the informations about the mesh to
> create (*.off)*/
> fichier = fopen("entry_deformable.off","r");
> fscanf(fichier," %d %d\n",&nbnoeuds, &nbtri);
> std::cout<<"nbnoeuds ="<<nbnoeuds<<'\n';
> std::cout<<"nbtri ="<<nbtri<<'\n';
>
> /*nodes coordinates reading*/
> for(int i=0; i<nbnoeuds; i++)
> {
> fscanf(fichier,"%lf %lf %lf", &noeudsAr[i][0], &noeudsAr[i][1],
> &noeudsAr[i][2]);
> }
>
> /*nodes in each triangle reading*/
> for(int i=0; i<nbtri; i++)
> {
> fscanf(fichier,"%i %i %i", &triAr[i][0], &triAr[i][1], &triAr[i][2]);
> }
>
> /*closing the file (*.off)*/
> fclose(fichier);
>
> /*******************************************************
> mesh creation from the noeudsAr and triAr
> ********************************************************/
>
> /*Characterization of the mesh*/
> const unsigned int PointDimension = 3;
> const unsigned int MaxTopologicalDimension = 2;
>
> typedef double PixelType; //must be the same as GradientPixelType
> typedef double CellDataType;
>
> typedef double CoordinateType;
> typedef double InterpolationWeightType;
>
> typedef itk::DefaultStaticMeshTraits<
> PixelType, PointDimension, MaxTopologicalDimension,
> CoordinateType, InterpolationWeightType, CellDataType > MeshTraits;
>
> typedef itk::Mesh< PixelType, PointDimension, MeshTraits > MeshType;
>
> typedef MeshType::CellType CellType;
> typedef itk::TriangleCell< CellType > TriangleType;
>
> MeshType::Pointer mesh = MeshType::New();
>
> /*points insertion*/
> typedef MeshType::PointType PointType;
> PointType point;
>
> const unsigned int numberOfPoints = nbnoeuds;
>
> for(unsigned int id=0; id<numberOfPoints; id++)
> {
> point[0] = noeudsAr[id][0];
> point[1] = noeudsAr[id][1];
> point[2] = noeudsAr[id][2];
> mesh->SetPoint( id, point );
> }
>
> /*cells inclusion*/
> CellType::CellAutoPointer triangle;
> const unsigned int numberOfCells = nbtri;
> for(unsigned int cellId=0; cellId<numberOfCells; cellId++)
> {
> triangle.TakeOwnership( new TriangleType );
> triangle->SetPointId( 0, triAr[cellId][0] ); // first point
> triangle->SetPointId( 1, triAr[cellId][1] ); // second point
> triangle->SetPointId( 2, triAr[cellId][2]); // third point
> mesh->SetCell( cellId, triangle ); // insert the cell
> }
> /*just a little check*/
> std::cout << "Points = " << mesh->GetNumberOfPoints() << std::endl;
> std::cout << "Cells = " << mesh->GetNumberOfCells() << std::endl;
>
> /****************************************************************
> gradient force field computing
> *****************************************************************/
>
> const unsigned int Dimension = 3;
> typedef double PixelType2;
> typedef itk::Image<PixelType2, Dimension> ImageType;
>
>
> typedef itk::CovariantVector< double, Dimension > GradientPixelType;
> typedef itk::Image< GradientPixelType, Dimension > GradientImageType;
>
>
> typedef itk::GradientRecursiveGaussianImageFilter<ImageType,
> GradientImageType>
> GradientFilterType;
> typedef
> itk::GradientMagnitudeRecursiveGaussianImageFilter<ImageType,ImageType>
> GradientMagnitudeFilterType;
>
> typedef itk::DeformableMesh3DFilter<MeshType,MeshType>
> DeformableFilterType;
>
> typedef itk::ImageFileReader< ImageType > ReaderType;
> ReaderType::Pointer imageReader = ReaderType::New();
>
> imageReader->SetFileName("brain.img");
>
>
> GradientMagnitudeFilterType::Pointer gradientMagnitudeFilter
> = GradientMagnitudeFilterType::New();
>
> gradientMagnitudeFilter->SetInput( imageReader->GetOutput() );
> gradientMagnitudeFilter->SetSigma( 1.0 );
>
> GradientFilterType::Pointer gradientMapFilter = GradientFilterType::New();
>
> gradientMapFilter->SetInput( gradientMagnitudeFilter->GetOutput());
> gradientMapFilter->SetSigma( 1.0 );
>
> std::cout << "Before map filter"<< std::endl;
>
> gradientMapFilter->Update();
>
> std::cout << "After map filter"<< std::endl;
>
> std::cout << "The gradient map created!" << std::endl;
>
> DeformableFilterType::Pointer deformableModelFilter =
> DeformableFilterType::New();
> deformableModelFilter->SetInput(mesh);
> deformableModelFilter->SetGradient(gradientMapFilter->GetOutput());
>
>
> typedef itk::CovariantVector<double, 2> double2DVector;
> typedef itk::CovariantVector<double, 3> double3DVector;
>
> double2DVector stiffness;
> stiffness[0] = 0.0001;
> stiffness[1] = 0.1;
>
> double3DVector scale;
> scale[0] = 0.98;
> scale[1] = 0.98;
> scale[2] = 1.56;
>
> deformableModelFilter->SetStiffness( stiffness );
> deformableModelFilter->SetScale( scale );
>
> deformableModelFilter->SetGradientMagnitude( 0.8 );
> deformableModelFilter->SetTimeStep( 0.01 );
> deformableModelFilter->SetStepThreshold( 5 );
>
> std::cout << "Deformable mesh fitting...";
>
> try
> {
> deformableModelFilter->Update(); //I think it bugs here
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception Caught !" << std::endl;
> std::cerr << excep << std::endl;
> }
>
> /**************************************************
> I'd like ti use the output of the last filter too...but it must work : (
> ***************************************************/
> return EXIT_SUCCESS;
> }
>
> *The file used (from ITK, but having been treated to erase nodes having *
> *the same coordinates) has this form :*
>
> for the nodes (coordinates x, y, z) :
>
> 96.47731 79.87897 0.00311
>
> 103.69083 79.68633 0.00311
>
> 91.25018 82.59742 0.00311
>
> for the cells (node1, node 2, node 3)
>
> 1598 1594 1718
>
> 1595 1600 1719
>
> 1599 1595 1719
>
> If someone knows what to do, please tell me...you'll save my life (or at
> least a good part of it)
>
> Thanks,
>
> Julien Mercenier, BELGIUM, University of Liège