[Insight-users] Problems using the ithDeformableMesh3DFilter
filter...need help!
Luis Ibanez
luis.ibanez at kitware.com
Fri, 05 Mar 2004 15:51:31 -0500
Hi Julien,
1) Given that you have a brain image and its deformed version
you could find the vector deformation field relating the two
images by using the DemonsRegistrationFilter.
http://www.itk.org/Insight/Doxygen/html/classitk_1_1DemonsRegistrationFilter.html
You will find a description of this filter in the SoftwareGuide
http://www.itk.org/ItkSoftwareGuide.pdf
2) The truncation of the identifies is not the reason for
the crash. With the DeformableMesh filter, chances are
that the forces are trying to make it evolve too fast, and
the geometrical consistence of the mesh gets disrupted.
3) The way to avoid big forces is to increase stiffness and
to reduce the step size. You also may want to write to
a file the vector force field that you are computing.
They usually have surprising features.
You may use ParaView for visualizing the Vector Field:
http://www.paraview.org
Note that you can also use the FEM-Deformable model
in order to find a deformation field mapping one brain
into another. This FEM model *will not* be equivalent
to the one that you actually want to create for modeling
the physical deformation, but it will make it possible
for you to compute a vector deformation field.
Regards,
Luis
-------------------
Julien Mercenier wrote:
> First,
> thanks for your answers...
>
> About your questions :
> "Are you bulding something like a finite elements model
> of the brain ?"
>
> Yes, I am.
>
> But, I have to take this sort of initial mesh (from Isosurf)
> to characterize the displacements from a brain to a deformed brain.
> (The segmentation is already done...)
> Therefore, I'd like to use forces derived from the deformed brain image
> (first a simple gradient, then balloon and then gradient vector flow)
> so the the itkDeformableMesh3DFilter filter seems to be the answer to this
> problem.
>
> (1) Is there thus possible to do this by the way I told ?
> Or do I have to reconsider the way I to solve this problem ?
>
> (2) Can the fact that some identifiers were truncated make the program
> crash ?
> If this is the case, what can I do (If you know the answers) ?
>
>
> (3) If I can use this filter, the way to avoid such big forces is
> to raise the stiffness ?
> but for the time step ?
>
>
> THANKS A LOT for providing such answers : )
> ----- Original Message -----
> From: "Luis Ibanez" <luis.ibanez at kitware.com>
> To: "Julien Mercenier" <itk_julienmercenier at hotmail.com>
> Cc: <insight-users at itk.org>
> Sent: Friday, March 05, 2004 4:36 AM
> Subject: Re: [Insight-users] Problems using the ithDeformableMesh3DFilter
> filter...need help!
>
>
>
>>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::Cell
> TraitsInfo<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::CellT
> raitsInfo<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
>>
>>
>>
>>_______________________________________________
>>Insight-users mailing list
>>Insight-users at itk.org
>>http://www.itk.org/mailman/listinfo/insight-users
>>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>