[Insight-users] Problems using the ithDeformableMesh3DFilter filter...need help!
Yan Pingkun
engp1734 at nus.edu.sg
Tue, 9 Mar 2004 11:26:29 +0800
Hi, Julien,
I have built your code and tried to run it. The problem is not caused by =
the itkDeformableMesh3DFilter. If you insert an exception:
try=20
{
gradientMapFilter->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cout << "Exception Caught !" << std::endl;
std::cout << excep << std::endl;
} =20
You may find that the image file 'sphere.img' couldn't be read into =
program. Then I tried to use RawImageIO but faild. I also tried to =
include a 'sphere.mhd' header file but it doesn't work either... The =
error is: Could not create IO object for file sphere.mhd
Hope others can solve this problem.
Regards,
Pingkun
-----Original Message-----
From: insight-users-admin at itk.org [mailto:insight-users-admin at itk.org] =
On Behalf Of Julien Mercenier
Sent: Sunday, March 07, 2004 8:10 AM
To: Luis Ibanez
Cc: insight-users at itk.org
Subject: Re: [Insight-users] Problems using the =
ithDeformableMesh3DFilter filter...need help!
Hi Luis,
before all : thank you for all your advices...
as other students working on the brain (registration, fusion,
vizualisation) with me
You really do a great job...
I've been working on the same problem for more than 5 days. I still =
have the same problem (even trying to tune all the parameters).
If someone could launch or read the code joined here...it would be =
so kind !!!!!!!!!!!!!
The image is a little sphere in *.img (*.hdr) format.
The initial mesh I create (with a lot of work) by hand is a cube. (49 =
nodes and 96 triangles).
The description of the bug I encounter is below.
As a reminder :
The little description of my system : - athlon1400MHz, 256 DDRAM, =
Windows XP
-ALL_BUILD - =
Win32Debug for ITK, VTK and all the programs I write
If someone finds the problem, please tell me...
A million thanks...
Julien
P.S : in the pipeline, there is a point I'm really not sure to be true. =
It is where I define the initial mesh of te DeformableMesh3DFilter "
deformableModelFilter->SetInput(mesh); //mesh is the mesh I creayed by =
hand (smartpointer)
deformableModelFilter->SetGradient(gradientMapFilter->GetOutput()); =
"
----- Original Message -----=20
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 9:51 PM
Subject: Re: [Insight-users] Problems using the =
ithDeformableMesh3DFilter filter...need help!
> 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_1DemonsRegistrationFil=
ter.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.=20
> > (The segmentation is already done...) Therefore, I'd like to use=20
> > forces derived from the deformed brain image (first a simple=20
> > gradient, then balloon and then gradient vector flow) so the the=20
> > 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=20
> >>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
> >>> -=20
> >>>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=20
> >>>create the initial mesh myself.
> >>> This mesh is the result coming from the program Isosurf giving=20
> >>>a surface mesh in 3D.
> >>> For the calculation of the force field, I use here the same=20
> >>>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) =3D> this is thus
> >>> a file of 8Mo). The pixel distances are 0.98, 0.98 and 1.56=20
> >>>(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::C=
ell
> > 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::Ce=
llT
> > 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=20
> >>>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=20
> >>>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=20
> >>>=3D " Unhandled exception in def_surf.exe : 0x0000005 : Access=20
> >>>Violation" in vxl/vnl/vnl_matrix_fixed.h indicating the " //get=20
> >>>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=20
> >>>'J:\WINDOWS\system32\MSVCP60D.DLL'
> >>>Loaded symbols for 'J:\WINDOWS\system32\MSVCRTD.DLL' Loaded=20
> >>>'J:\WINDOWS\system32\user32.dll', no matching symbolic information=20
> >>>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=20
> >>>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**************************=20
> >>>**************************************************************/
> >>>
> >>>
> >>>/*************************************************************
> >>>----------------
> >>>Several steps =3D=3D>
> >>>----------------
> >>> 1) use the file to create the mesh
> >>> 2) computing the force field (gradient)
> >>> 3) mesh and forxce field into the DeformableMeshFilter=20
> >>>***************************************************************/
> >>>
> >>>#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 =3D 0;
> >>> int nbtri =3D 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 =3D 0;
> >>>
> >>> /*opening of the file containing the informations about the mesh
> >>>to create (*.off)*/ fichier =3D fopen("entry_deformable.off","r");
> >>> fscanf(fichier," %d %d\n",&nbnoeuds, &nbtri);
> >>> std::cout<<"nbnoeuds =3D"<<nbnoeuds<<'\n';
> >>> std::cout<<"nbtri =3D"<<nbtri<<'\n';
> >>>
> >>> /*nodes coordinates reading*/
> >>> for(int i=3D0; 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=3D0; 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 =3D 3;
> >>> const unsigned int MaxTopologicalDimension =3D 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 =3D MeshType::New();
> >>>
> >>> /*points insertion*/
> >>> typedef MeshType::PointType PointType;
> >>> PointType point;
> >>>
> >>> const unsigned int numberOfPoints =3D nbnoeuds;
> >>>
> >>> for(unsigned int id=3D0; id<numberOfPoints; id++)
> >>> {
> >>> point[0] =3D noeudsAr[id][0];
> >>> point[1] =3D noeudsAr[id][1];
> >>> point[2] =3D noeudsAr[id][2];
> >>> mesh->SetPoint( id, point );
> >>> }
> >>>
> >>> /*cells inclusion*/
> >>> CellType::CellAutoPointer triangle;
> >>> const unsigned int numberOfCells =3D nbtri;
> >>> for(unsigned int cellId=3D0; cellId<numberOfCells; cellId++) { =20
> >>> triangle.TakeOwnership( new TriangleType ); =20
> >>> triangle->SetPointId( 0, triAr[cellId][0] ); // first point =20
> >>> triangle->SetPointId( 1, triAr[cellId][1] ); // second point =20
> >>> triangle->SetPointId( 2, triAr[cellId][2]); // third point
> >>> mesh->SetCell( cellId, triangle ); // insert the cell
> >>> }
> >>> /*just a little check*/
> >>> std::cout << "Points =3D " << mesh->GetNumberOfPoints() <<=20
> >>> std::endl; std::cout << "Cells =3D " << mesh->GetNumberOfCells() =
> >>> << std::endl;
> >>>
> >>> /****************************************************************
> >>> gradient force field computing
> >>> *****************************************************************/
> >>>
> >>> const unsigned int Dimension =3D 3;
> >>> typedef double PixelType2;
> >>> typedef itk::Image<PixelType2, Dimension> ImageType;
> >>>
> >>>
> >>> typedef itk::CovariantVector< double, Dimension >
> >>> GradientPixelType; typedef itk::Image< GradientPixelType,=20
> >>> Dimension > GradientImageType;
> >>>
> >>>
> >>> typedef itk::GradientRecursiveGaussianImageFilter<ImageType,
> >>>GradientImageType>
> >>> GradientFilterType;
> >>> typedef
> >>>itk::GradientMagnitudeRecursiveGaussianImageFilter<ImageType,ImageT
> >>>ype>
> >>> GradientMagnitudeFilterType;
> >>>
> >>> typedef itk::DeformableMesh3DFilter<MeshType,MeshType>
> >>>DeformableFilterType;
> >>>
> >>> typedef itk::ImageFileReader< ImageType > ReaderType;
> >>> ReaderType::Pointer imageReader =3D ReaderType::New();
> >>>
> >>> imageReader->SetFileName("brain.img");
> >>>
> >>>
> >>> GradientMagnitudeFilterType::Pointer gradientMagnitudeFilter
> >>> =3D GradientMagnitudeFilterType::New();
> >>>
> >>> gradientMagnitudeFilter->SetInput( imageReader->GetOutput() );
> >>> gradientMagnitudeFilter->SetSigma( 1.0 );
> >>>
> >>> GradientFilterType::Pointer gradientMapFilter =3D
> >
> > GradientFilterType::New();
> >
> >>> gradientMapFilter->SetInput(
> >>> gradientMapFilter->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 =3D
> >>> DeformableFilterType::New();
> >>> deformableModelFilter->SetInput(mesh);
> >>> deformableModelFilter->SetGradient(gradientMapFilter->GetOutput())
> >>> deformableModelFilter->;
> >>>
> >>>
> >>> typedef itk::CovariantVector<double, 2> double2DVector;
> >>> typedef itk::CovariantVector<double, 3> double3DVector;
> >>>
> >>> double2DVector stiffness;
> >>> stiffness[0] =3D 0.0001;
> >>> stiffness[1] =3D 0.1;
> >>>
> >>> double3DVector scale;
> >>> scale[0] =3D 0.98;
> >>> scale[1] =3D 0.98;
> >>> scale[2] =3D 1.56;
> >>>
> >>> deformableModelFilter->SetStiffness( stiffness ); SetScale( scale
> >>> deformableModelFilter->);
> >>>
> >>> deformableModelFilter->SetGradientMagnitude( 0.8 ); SetTimeStep(
> >>> deformableModelFilter->0.01 ); 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=E8ge
> >>
> >>
> >>
> >>_______________________________________________
> >>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
> >
>
>
>
>