[Insight-users] itkDeformableMesh3DFilter: Suggestion for code changes
Thomas Boettger
t . boettger at dkfz-heidelberg . de
Tue, 05 Aug 2003 16:28:00 +0200
This is a multi-part message in MIME format.
--------------060502080703090308040000
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit
Hello everybody,
I already posted some questions about the itkDeformableMesh3DFilter
yesterday. Unfortunately nobody could tell me on which scientific paper
the underlying algorithm is based, which I'd still like to know. (Maybe
its the wrong mailing list here and someone could forward the mail to
the developers list.)
Well I was playing around with the filter, trying to understand the code
and modified the sources so the filter produced my expected results.
I'd really like know more about the algorithm and if necessary discuss
the suggested changes with the responsible developers, because I'm
interested in using this part of ITK as well.
I attached the part of the source code which I modified. The new changes
are made in method ComputeDt() and the changes described yesterday (see
attached mail file from yesterday) are made in GenerateData().
Thank you for your help.
Thomas Boettger
--
Dipl.-Inform. Thomas Boettger
Deutsches Krebsforschungszentrum (German Cancer Research Center)
Div. Medical and Biological Informatics H0100 Tel: (+49) 6221-42 2328
Im Neuenheimer Feld 280 Fax: (+49) 6221-42 2345
D-69120 Heidelberg e-mail: t . boettger at dkfz . de
Germany http://www . dkfz . de/mbi/people/thomasb . shtml
--------------060502080703090308040000
Content-Type: text/plain;
name="modifications_in_itkDeformableMesh3DFilter.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="modifications_in_itkDeformableMesh3DFilter.txt"
/* Compute the derivatives using d'- Kd = f. */
template <typename TInputMesh, typename TOutputMesh>
void
DeformableMesh3DFilter<TInputMesh, TOutputMesh>
::ComputeDt()
{
int i;
const unsigned long *tp;
InputCellsContainerPointer myCells = this->GetInput(0)->GetCells();
InputCellsContainerIterator cells = myCells->Begin();
InputPointsContainerPointer myForces = m_Forces->GetPoints();
InputPointsContainerIterator forces = myForces->Begin();
InputPointsContainerPointer myDerives = m_Derives->GetPoints();
InputPointsContainerIterator derives = myDerives->Begin();
double p = 1.0;
i = 0;
InputPointType v1, v2, v3;
InputPointType* v1_pt;
InputPointType* v2_pt;
InputPointType* v3_pt;
v1_pt = &v1;
v2_pt = &v2;
v3_pt = &v3;
while( cells != myCells->End() ) {
tp = cells.Value()->GetPointIds();
++cells;
m_Displacements->GetPoint (tp[0], v1_pt);
m_Displacements->GetPoint (tp[1], v2_pt);
m_Displacements->GetPoint (tp[2], v3_pt);
v1[0] *= m_K[i]->get(0, 0)*p;
v1[1] *= m_K[i]->get(0, 0)*p;
v1[2] *= m_K[i]->get(0, 0)*p;
v2[0] *= m_K[i]->get(0, 1)*p;
v2[1] *= m_K[i]->get(0, 1)*p;
v2[2] *= m_K[i]->get(0, 1)*p;
v3[0] *= m_K[i]->get(0, 2)*p;
v3[1] *= m_K[i]->get(0, 2)*p;
v3[2] *= m_K[i]->get(0, 2)*p;
v1[0] += v2[0]+v3[0];
v1[1] += v2[1]+v3[1];
v1[2] += v2[2]+v3[2];
m_Forces->GetPoint (tp[0], v2_pt);
/***********************************************************************
Thought we want to solve: d' - K*d = f --> += ?
Shouldn't we use an addition instead of the substraction?
**********************************************************************/
//v2[0] -= v1[0];
//v2[1] -= v1[1];
//v2[2] -= v1[2];
v2[0] += v1[0];
v2[1] += v1[1];
v2[2] += v1[2];
/***********************************************************************
shouldn't the computed derives be set on the derives member???
When the force for this point is needed a second time (which can
actually happen, because a point can be member of many cells)
- in my opinion - an algorithmic error will occur.
**********************************************************************/
//m_Forces->SetPoint (tp[0], v2);
m_Derives->SetPoint (tp[1], v2);
m_Displacements->GetPoint (tp[0], v1_pt);
m_Displacements->GetPoint (tp[1], v2_pt);
m_Displacements->GetPoint (tp[2], v3_pt);
v1[0] *= m_K[i]->get(1, 0)*p;
v1[1] *= m_K[i]->get(1, 0)*p;
v1[2] *= m_K[i]->get(1, 0)*p;
v2[0] *= m_K[i]->get(1, 1)*p;
v2[1] *= m_K[i]->get(1, 1)*p;
v2[2] *= m_K[i]->get(1, 1)*p;
v3[0] *= m_K[i]->get(1, 2)*p;
v3[1] *= m_K[i]->get(1, 2)*p;
v3[2] *= m_K[i]->get(1, 2)*p;
v1[0] += v2[0]+v3[0];
v1[1] += v2[1]+v3[1];
v1[2] += v2[2]+v3[2];
m_Forces->GetPoint (tp[1], v2_pt);
// ************ same as above ***************
//v2[0] -= v1[0];
//v2[1] -= v1[1];
//v2[2] -= v1[2];
v2[0] += v1[0];
v2[1] += v1[1];
v2[2] += v1[2];
// ************ same as above ***************
//m_Forces->SetPoint (tp[0], v2);
m_Derives->SetPoint (tp[1], v2);
m_Displacements->GetPoint (tp[0], v1_pt);
m_Displacements->GetPoint (tp[1], v2_pt);
m_Displacements->GetPoint (tp[2], v3_pt);
v1[0] *= m_K[i]->get(2, 0)*p;
v1[1] *= m_K[i]->get(2, 0)*p;
v1[2] *= m_K[i]->get(2, 0)*p;
v2[0] *= m_K[i]->get(2, 1)*p;
v2[1] *= m_K[i]->get(2, 1)*p;
v2[2] *= m_K[i]->get(2, 1)*p;
v3[0] *= m_K[i]->get(2, 2)*p;
v3[1] *= m_K[i]->get(2, 2)*p;
v3[2] *= m_K[i]->get(2, 2)*p;
v1[0] += v2[0]+v3[0];
v1[1] += v2[1]+v3[1];
v1[2] += v2[2]+v3[2];
m_Forces->GetPoint (tp[2], v2_pt);
// ************ same as above ***************
//v2[0] -= v1[0];
//v2[1] -= v1[1];
//v2[2] -= v1[2];
v2[0] += v1[0];
v2[1] += v1[1];
v2[2] += v1[2];
// ************ same as above ***************
//m_Forces->SetPoint (tp[2], v2);
m_Derives->SetPoint (tp[2], v2);
++i;
}
/***********************************************************************
This part isn't necessary anymore due to the changes above...
**********************************************************************/
//while ( derives != myDerives->End() ) {
//derives.Value() = forces.Value();
//++derives;
//++forces;
//}
}
/* Generate Data */
template <typename TInputMesh, typename TOutputMesh>
void
DeformableMesh3DFilter<TInputMesh, TOutputMesh>
::GenerateData()
{
this->Initialize();
/***********************************************************************
Initialize the whole array so we can get pointers to correct data in
method SetMeshStiffness()...
I think its just a workaround an should be managed in another way,
if a solution is found.
************************************************************************/
for (int i=1; i< 10;i++)
{
SetStiffnessMatrix( m_StiffnessMatrix[0] ,i);
}
this->SetMeshStiffness();
....
}
--------------060502080703090308040000
Content-Type: text/plain;
name="mail_from2003_08_04_possible_error_in_DeformableMeshFilter.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="mail_from2003_08_04_possible_error_in_DeformableMeshFilter.txt"
Subject:
[Insight-users] Possible Error in itkDeformableMesh3DFilter
From:
Thomas Boettger <t . boettger at dkfz-heidelberg . de>
Date:
Mon, 04 Aug 2003 14:04:32 +0200
To:
"insight-users at public . kitware . com" <insight-users at public . kitware . com>
CC:
Ting Chen <chenting at graphics . cis . upenn . edu>
Hello,
I am currently writing a new filter based on the class itkDeformableMesh3DFilter. As I posted a few months ago I thought that the stiffness parameters are not taken into account. I downloaded the latest CVS version of the class and still experienced the same problem.
After debugging my program I found the following code in method
DeformableMesh3DFilter::SetMeshStiffness():
...
while (celldata != myCellData->End()){
x = celldata.Value();
m_K[j] = m_StiffnessMatrix+((int) x); <----- critical
++celldata;
j++;
}
...
I am not sure what the developer intended to do when adding the celldata value to the current pointer in m_StiffnessMatrix. But only the first pointer of m_StiffnessMatrix[10] is currently intialized at all. So by performing the pointer arithmetic the m_K[j] pointers all contain pointers to unintialized values. Strange that the filter produced such senseful results under Windows anyway. When running the code on linux it somtimes produced very strange results. Thats where I found the uninitialised values usage.
After a first modification in the itkDeformableMeshFilter where I initializes the whole m_StiffnessMatrix[] array. I saw the influence of the stiffness parameters to the segmentation result for the first time.
I am not sure how this bug should be handled, if it should just be fixed or if the algorithm should be checked for correctness.
My questions:
1. I would like to know on which paper this implementation is based.
2. What is the algorithmic sense of the applied pointer addition in SetMeshStiffness().
3. Or should the celldata value be added to the elements of the matrix?
4. What is the meaning of the celldata value anyway?
5. Why are the 10 stiffness matrices?
Thanks in advance,
Thomas Boettger
--
Dipl.-Inform. Thomas Boettger
Deutsches Krebsforschungszentrum (German Cancer Research Center)
Div. Medical and Biological Informatics H0100 Tel: (+49) 6221-42 2328
Im Neuenheimer Feld 280 Fax: (+49) 6221-42 2345
D-69120 Heidelberg e-mail: t . boettger at dkfz . de
Germany http://www . dkfz . de/mbi/people/thomasb . shtml
_______________________________________________
Insight-users mailing list
Insight-users at itk . org
http://www . itk . org/mailman/listinfo/insight-users
--------------060502080703090308040000--