[Insight-users] Problems using the ithDeformableMesh3DFilter filter...need help!
Julien Mercenier
itk_julienmercenier at hotmail.com
Thu, 4 Mar 2004 15:58:03 +0100
This is a multi-part message in MIME format.
------=_NextPart_000_0033_01C40201.7A16E480
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
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) =3D> this is thus=20
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 :=20
J:\Program Files\Microsoft Visual Studio\VC98\INCLUDE\xmemory(64) : =
warning C4786: =
'std::allocator<itk::SmartPointer<itk::CellInterfaceVisitor<double,itk::C=
ellTraitsInfo<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=
llTraitsInfo<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 =3D " 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 =3D=3D>
----------------=20
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>=20
#include <iostream>=20
#include "itkMesh.h" =20
#include "itkDefaultStaticMeshTraits.h"
#include "itkTriangleCell.h"
#include "itkDeformableMesh3DFilter.h"=20
#include "itkGradientRecursiveGaussianImageFilter.h"=20
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"=20
#include "itkImage.h"
#include "itkCovariantVector.h"//will contain the force field
#include "itkImageFileReader.h"=20
int main()
{
=20
int nbnoeuds =3D 0;
int nbtri =3D 0;
double noeudsAr[3648][3]; //3648 is the number of nodes=20
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]);
}
=20
/*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; =20
typedef double InterpolationWeightType;
typedef itk::DefaultStaticMeshTraits<=20
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++)
{
triangle.TakeOwnership( new TriangleType );
triangle->SetPointId( 0, triAr[cellId][0] ); // first point=20
triangle->SetPointId( 1, triAr[cellId][1] ); // second point
triangle->SetPointId( 2, triAr[cellId][2]); // third point=20
mesh->SetCell( cellId, triangle ); // insert the cell=20
} =20
/*just a little check*/
std::cout << "Points =3D " << mesh->GetNumberOfPoints() << 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, 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 =3D ReaderType::New();
imageReader->SetFileName("brain.img");
GradientMagnitudeFilterType::Pointer gradientMagnitudeFilter
=3D GradientMagnitudeFilterType::New();
gradientMagnitudeFilter->SetInput( imageReader->GetOutput() );=20
gradientMagnitudeFilter->SetSigma( 1.0 );
GradientFilterType::Pointer gradientMapFilter =3D =
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 =3D=20
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] =3D 0.0001;
stiffness[1] =3D 0.1;
double3DVector scale;
scale[0] =3D 0.98;
scale[1] =3D 0.98;=20
scale[2] =3D 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=20
{
deformableModelFilter->Update(); //I think it bugs here
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception Caught !" << std::endl;
std::cerr << excep << std::endl;
}
=20
/**************************************************
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=20
the same coordinates) has this form :
for the nodes (coordinates x, y, z) :=20
96.47731 79.87897 0.00311=20
103.69083 79.68633 0.00311=20
91.25018 82.59742 0.00311=20
for the cells (node1, node 2, node 3)
1598 1594 1718=20
1595 1600 1719=20
1599 1595 1719=20
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
------=_NextPart_000_0033_01C40201.7A16E480
Content-Type: text/html;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=3DContent-Type content=3D"text/html; =
charset=3Diso-8859-1">
<META content=3D"MSHTML 6.00.2800.1400" name=3DGENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=3D#ffffff>
<DIV><FONT face=3DArial size=3D2>Hi everybody,</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2> Here is a little =
description of=20
my system : - athlon1400MHz, 256 DDRAM, Windows XP</FONT></DIV>
<DIV><FONT face=3DArial=20
size=3D2> &nbs=
p;  =
; =
&=
nbsp; &n=
bsp; =20
- ALL_BUILD - Win32Debug for ITK, VTK and all the programs I =
write</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2> I want to =
characterize the=20
displacements of the brain surface in 3D. </FONT></DIV>
<DIV><FONT face=3DArial size=3D2> To do this, I use =
the ITK=20
User-Guide exemple using the ithDeformableMesh3DFilter filter (page=20
379).</FONT></DIV>
<DIV><FONT face=3DArial size=3D2> But, except of using =
the itkBinaryMask3DMeshSource filter, I create the initial mesh=20
myself.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2> This mesh is the =
result coming=20
from the program Isosurf giving a surface mesh in 3D.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2> For the calculation =
of the force=20
field, I use here the same method as in the ITK User-guide =
example.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2> The binary image =
used for this=20
calculation is a deformed brain, in *.img format (*.hdr as =
header)</FONT></DIV>
<DIV><FONT face=3DArial size=3D2> (the dimension of =
this image is=20
256*256*128 with 00 (not in the brain) or 01(in the brain) =3D> this =
is thus=20
</FONT></DIV>
<DIV><FONT face=3DArial size=3D2> a file of 8Mo). The =
pixel=20
distances are 0.98, 0.98 and 1.56 (mm).</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2> There are several =
things I don't=20
understant or I can't solve by myself.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2> - <STRONG>At=20
compiling-time, there are 56 warnings (all the sames but for =
different=20
variables).</STRONG> Here is an example :</FONT> </DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>J:\Program Files\Microsoft Visual=20
Studio\VC98\INCLUDE\xmemory(64) : warning C4786:=20
'std::allocator<itk::SmartPointer<itk::CellInterfaceVisitor<doub=
le,itk::CellTraitsInfo<3,double,double,unsigned=20
long,unsigned long,unsigned=20
long,itk::Point<double,3>,<BR>itk::VectorContainer<unsigned=20
long,itk::Point<double,3> >,std::set<unsigned=20
long,std::less<unsigned long>,std::allocator<unsigned long> =
>=20
> > > >' : identifier was truncated to '255' characters in =
the debug=20
information<BR> J:\Program=20
Files\Microsoft Visual Studio\VC98\INCLUDE\xmemory(64) : while compiling =
class-template member function 'void __thiscall=20
std::allocator<itk::SmartPointer<itk::CellInterfaceVisitor<doubl=
e,itk::CellTraitsInfo<3,double,double,unsigne<BR>d=20
long,unsigned long,unsigned=20
long,itk::Point<double,3>,itk::VectorContainer<unsigned=20
long,itk::Point<double,3> >,std::set<unsigned=20
long,std::less<unsigned long>,std::allocator<unsigned long> =
>=20
> > > >::deallocate(void *,unsigned int)'</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2> - <STRONG>At=20
running-time</STRONG></FONT></DIV>
<DIV><FONT face=3DArial size=3D2>The program stops when using the=20
ithDeformableMesh3DFilter filter : there is a bug !! But I don't =
understand=20
why...even when reading the codes.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2><U>Then, am message apperas in a =
window. And, when=20
I click on the debug button, a message appears in the visualC++ window=20
:</U></FONT></DIV>
<DIV><FONT face=3DArial size=3D2>message =3D " Unhandled exception in =
def_surf.exe :=20
0x0000005 : Access Violation"</FONT></DIV>
<DIV><FONT face=3DArial size=3D2>in =
vxl/vnl/vnl_matrix_fixed.h indicating=20
the " //get element " function.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2><U>In the debug window appears the =
messages=20
:</U></FONT></DIV>
<DIV><FONT face=3DArial size=3D2>Loaded symbols for=20
'J:\DEF_SURF\Debug\def_surf.exe'<BR>Loaded =
'J:\WINDOWS\system32\ntdll.dll', no=20
matching symbolic information found.<BR>Loaded=20
'J:\WINDOWS\system32\kernel32.dll', no matching symbolic information=20
found.<BR>Loaded symbols for =
'J:\WINDOWS\system32\MSVCP60D.DLL'<BR>Loaded=20
symbols for 'J:\WINDOWS\system32\MSVCRTD.DLL'<BR>Loaded=20
'J:\WINDOWS\system32\user32.dll', no matching symbolic information=20
found.<BR>Loaded 'J:\WINDOWS\system32\gdi32.dll', no matching symbolic=20
information found.<BR>Loaded 'J:\WINDOWS\system32\advapi32.dll', no =
matching=20
symbolic information found.<BR>Loaded 'J:\WINDOWS\system32\rpcrt4.dll', =
no=20
matching symbolic information found.<BR>Loaded=20
'J:\WINDOWS\system32\version.dll', no matching symbolic information=20
found.<BR>Loaded 'J:\WINDOWS\system32\msvcrt.dll', no matching symbolic=20
information found.<BR>Loaded 'J:\WINDOWS\system32\SHLWAPI.DLL', no =
matching=20
symbolic information found.<BR>Loaded 'J:\WINDOWS\system32\apphelp.dll', =
no=20
matching symbolic information found.<BR>The thread 0xBB0 has exited with =
code 0=20
(0x0).<BR>First-chance exception in def_surf.exe: 0xC0000005: Access=20
Violation.<BR>The program 'J:\DEF_SURF\Debug\def_surf.exe' has exited =
with code=20
0 (0x0).</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2><STRONG>I put my code just here =
</STRONG>(the=20
creation of the mesh don't seem to pose any problem)</FONT></DIV>
<DIV><STRONG><FONT face=3DArial size=3D2></FONT></STRONG> </DIV>
<DIV><FONT face=3DArial=20
size=3D2>/*************************************************************<B=
R>*******active=20
deformable=20
surface**************************<BR>************************************=
**************************/</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV><FONT face=3DArial =
size=3D2>
<DIV><BR>/*************************************************************<B=
R>----------------<BR>Several=20
steps =3D=3D><BR>---------------- <BR> 1) use the =
file to=20
create the mesh<BR> 2) computing the force field=20
(gradient)<BR> 3) mesh and forxce field into the=20
DeformableMeshFilter<BR>*************************************************=
**************/</DIV>
<DIV> </DIV>
<DIV>#include <cstdio> //C-style I/O to read the Isosurf=20
file<BR>#include <stdio.h> </DIV>
<DIV> </DIV>
<DIV>#include <iostream> </DIV>
<DIV> </DIV>
<DIV>#include=20
"itkMesh.h" &n=
bsp; =20
<BR>#include "itkDefaultStaticMeshTraits.h"<BR>#include=20
"itkTriangleCell.h"</DIV>
<DIV> </DIV>
<DIV>#include "itkDeformableMesh3DFilter.h" </DIV>
<DIV> </DIV>
<DIV>#include "itkGradientRecursiveGaussianImageFilter.h" <BR>#include=20
"itkGradientMagnitudeRecursiveGaussianImageFilter.h" </DIV>
<DIV> </DIV>
<DIV>#include "itkImage.h"</DIV>
<DIV> </DIV>
<DIV>#include "itkCovariantVector.h"//will contain the force field</DIV>
<DIV> </DIV>
<DIV>#include "itkImageFileReader.h" </DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV>int main()<BR>{<BR> <BR> int nbnoeuds =3D 0;<BR> int =
nbtri =3D=20
0;</DIV>
<DIV> </DIV>
<DIV> double=20
noeudsAr[3648][3]; =
//3648=20
is the number of nodes <BR> double =
noeudsdefAr[3648][3];<BR> int=20
triAr[7280][3]; &nbs=
p;  =
;=20
//7280 is the number of triangles</DIV>
<DIV> </DIV>
<DIV> FILE* fichier =3D 0;</DIV>
<DIV> </DIV>
<DIV> /*opening of the file containing the informations about the =
mesh to=20
create (*.off)*/<BR> fichier =3D=20
fopen("entry_deformable.off","r");<BR> fscanf(fichier," %d=20
%d\n",&nbnoeuds, &nbtri);<BR> std::cout<<"nbnoeuds=20
=3D"<<nbnoeuds<<'\n';<BR> std::cout<<"nbtri=20
=3D"<<nbtri<<'\n';</DIV>
<DIV> </DIV>
<DIV> /*nodes coordinates reading*/<BR> for(int i=3D0; =
i<nbnoeuds;=20
i++)<BR> {<BR> fscanf(fichier,"%lf %lf %lf",=20
&noeudsAr[i][0], &noeudsAr[i][1], =
&noeudsAr[i][2]);<BR> }</DIV>
<DIV> </DIV>
<DIV> /*nodes in each triangle reading*/<BR> for(int i=3D0; =
i<nbtri;=20
i++)<BR> {<BR> fscanf(fichier,"%i %i %i", =
&triAr[i][0],=20
&triAr[i][1], =
&triAr[i][2]);<BR> }<BR> <BR> /*closing the=20
file (*.off)*/<BR> fclose(fichier);</DIV>
<DIV> </DIV>
<DIV> /*******************************************************<BR>&n=
bsp;mesh=20
creation from the noeudsAr and=20
triAr<BR> ********************************************************/<=
/DIV>
<DIV> </DIV>
<DIV> /*Characterization of the mesh*/<BR> const unsigned int=20
PointDimension =3D 3;<BR> const unsigned int =
MaxTopologicalDimension =3D=20
2;</DIV>
<DIV> </DIV>
<DIV> typedef double PixelType; //must be the same as=20
GradientPixelType<BR> typedef double CellDataType;</DIV>
<DIV> </DIV>
<DIV> typedef double=20
CoordinateType; <BR> =
typedef=20
double InterpolationWeightType;</DIV>
<DIV> </DIV>
<DIV> typedef itk::DefaultStaticMeshTraits<=20
<BR> PixelType, PointDimension,=20
MaxTopologicalDimension,<BR> CoordinateType,=20
InterpolationWeightType, CellDataType > MeshTraits;</DIV>
<DIV> </DIV>
<DIV> typedef itk::Mesh< PixelType, PointDimension, MeshTraits =
>=20
MeshType;</DIV>
<DIV> </DIV>
<DIV> typedef=20
MeshType::CellType &=
nbsp; =20
CellType;<BR> typedef itk::TriangleCell< CellType=20
> TriangleType;</DIV>
<DIV> </DIV>
<DIV> MeshType::Pointer mesh =3D MeshType::New();</DIV>
<DIV> </DIV>
<DIV> /*points insertion*/<BR> typedef MeshType::PointType=20
PointType;<BR> PointType point;</DIV>
<DIV> </DIV>
<DIV> const unsigned int numberOfPoints =3D nbnoeuds;</DIV>
<DIV> </DIV>
<DIV> for(unsigned int id=3D0; id<numberOfPoints;=20
id++)<BR> {<BR> point[0] =3D=20
noeudsAr[id][0];<BR> point[1] =3D=20
noeudsAr[id][1];<BR> point[2] =3D=20
noeudsAr[id][2];<BR> mesh->SetPoint( id, point =
);<BR> }</DIV>
<DIV> </DIV>
<DIV> /*cells inclusion*/<BR> CellType::CellAutoPointer=20
triangle;<BR> const unsigned int numberOfCells =3D=20
nbtri;<BR> for(unsigned int cellId=3D0; cellId<numberOfCells;=20
cellId++)<BR> {<BR> triangle.TakeOwnership( new=20
TriangleType );<BR> triangle->SetPointId( 0,=20
triAr[cellId][0] ); // first point =
<BR> triangle->SetPointId( 1,=20
triAr[cellId][1] ); // second =
point<BR> triangle->SetPointId( 2,=20
triAr[cellId][2]); // third point <BR> mesh->SetCell( =
cellId,=20
triangle ); // insert the cell=20
<BR> } <BR> /*just a little =
check*/<BR> std::cout=20
<< "Points =3D " << mesh->GetNumberOfPoints() <<=20
std::endl;<BR> std::cout << "Cells =3D " <<=20
mesh->GetNumberOfCells() << std::endl;</DIV>
<DIV> </DIV>
<DIV> /*************************************************************=
***<BR> gradient=20
force field=20
computing<BR> ******************************************************=
***********/</DIV>
<DIV> </DIV>
<DIV> const unsigned int =20
Dimension =3D 3;<BR> typedef =20
double =20
PixelType2;<BR> typedef itk::Image<PixelType2,=20
Dimension> ImageType;</DIV>
<DIV> </DIV>
<DIV><BR> typedef itk::CovariantVector< double, Dimension =
> =20
GradientPixelType;<BR> typedef itk::Image< GradientPixelType, =
Dimension=20
> GradientImageType;</DIV>
<DIV> </DIV>
<DIV><BR> typedef =
itk::GradientRecursiveGaussianImageFilter<ImageType,=20
GradientImageType><BR> GradientFilterType;<BR> typedef=
=20
itk::GradientMagnitudeRecursiveGaussianImageFilter<ImageType,ImageType=
><BR> GradientMagnitudeFilterType;</DIV>
<DIV> </DIV>
<DIV> typedef =
itk::DeformableMesh3DFilter<MeshType,MeshType> =20
DeformableFilterType;</DIV>
<DIV> </DIV>
<DIV> typedef itk::ImageFileReader<=20
ImageType > =20
ReaderType;<BR> ReaderType::Pointer &nb=
sp;=20
imageReader =3D ReaderType::New();</DIV>
<DIV> </DIV>
<DIV> imageReader->SetFileName("brain.img");</DIV>
<DIV> </DIV>
<DIV><BR> GradientMagnitudeFilterType::Pointer =20
gradientMagnitudeFilter<BR> &nbs=
p; =3D=20
GradientMagnitudeFilterType::New();</DIV>
<DIV> </DIV>
<DIV> gradientMagnitudeFilter->SetInput( =
imageReader->GetOutput() );=20
<BR> gradientMagnitudeFilter->SetSigma( 1.0 );</DIV>
<DIV> </DIV>
<DIV> GradientFilterType::Pointer gradientMapFilter =3D=20
GradientFilterType::New();</DIV>
<DIV> </DIV>
<DIV> gradientMapFilter->SetInput(=20
gradientMagnitudeFilter->GetOutput());<BR> gradientMapFilter->=
SetSigma(=20
1.0 );</DIV>
<DIV><BR> std::cout << "Before map filter"<< =
std::endl;</DIV>
<DIV><BR> gradientMapFilter->Update();</DIV>
<DIV><BR> std::cout << "After map filter"<< =
std::endl;</DIV>
<DIV> </DIV>
<DIV> std::cout << "The gradient map created!" <<=20
std::endl;</DIV>
<DIV> </DIV>
<DIV> DeformableFilterType::Pointer deformableModelFilter =3D=20
<BR> &nb=
sp; &nbs=
p; =20
DeformableFilterType::New();<BR> deformableModelFilter->SetInput(=
mesh);<BR> deformableModelFilter->SetGradient(gradientMapFilter-&=
gt;GetOutput());</DIV>
<DIV> </DIV>
<DIV><BR> typedef itk::CovariantVector<double,=20
2> =20
double2DVector;<BR> typedef itk::CovariantVector<double,=20
3> =20
double3DVector;</DIV>
<DIV> </DIV>
<DIV> double2DVector stiffness;<BR> stiffness[0] =3D=20
0.0001;<BR> stiffness[1] =3D 0.1;</DIV>
<DIV> </DIV>
<DIV> double3DVector scale;<BR> scale[0] =3D =
0.98;<BR> scale[1] =3D=20
0.98; <BR> scale[2] =3D 1.56;</DIV>
<DIV> </DIV>
<DIV> deformableModelFilter->SetStiffness( stiffness=20
);<BR> deformableModelFilter->SetScale( scale );</DIV>
<DIV> </DIV>
<DIV> deformableModelFilter->SetGradientMagnitude( 0.8=20
);<BR> deformableModelFilter->SetTimeStep( 0.01=20
);<BR> deformableModelFilter->SetStepThreshold( 5 );</DIV>
<DIV> </DIV>
<DIV> std::cout << "Deformable mesh fitting...";</DIV>
<DIV> </DIV>
<DIV> try=20
<BR> {<BR> deformableModelFilter->Update(); =
; =20
//I think it bugs here<BR> }<BR> catch( =
itk::ExceptionObject=20
& excep )<BR> {<BR> std::cerr << =
"Exception=20
Caught !" << std::endl;<BR> std::cerr << excep =
<<=20
std::endl;<BR> }<BR> <BR> /*************************=
*************************</DIV>
<DIV>I'd like ti use the output of the last filter too...but it must =
work :=20
(</DIV>
<DIV>***************************************************/</DIV>
<DIV> return EXIT_SUCCESS;<BR>}</DIV>
<DIV> </DIV>
<DIV><STRONG>The file used (from ITK, but having been treated to erase =
nodes=20
having </STRONG></DIV>
<DIV><STRONG>the same coordinates) has this form :</STRONG></DIV>
<DIV> </DIV>
<DIV>for the nodes (coordinates x, y, z) : <FONT size=3D2>
<P>96.47731 79.87897 0.00311 </P>
<P>103.69083 79.68633 0.00311 </P>
<P>91.25018 82.59742 0.00311 </P></FONT></DIV>
<DIV>for the cells (node1, node 2, node 3)</DIV>
<DIV><FONT size=3D2>
<P>1598 1594 1718 </P>
<P>1595 1600 1719 </P>
<P>1599 1595 1719 </P></FONT></DIV>
<DIV>If someone knows what to do, please tell me...you'll save my life =
(or at=20
least a good part of it)</DIV>
<DIV> </DIV>
<DIV>Thanks,</DIV>
<DIV> </DIV>
<DIV> Julien Mercenier, BELGIUM, University of=20
Li=E8ge</FONT></DIV></BODY></HTML>
------=_NextPart_000_0033_01C40201.7A16E480--