Hi
Agatte,<div><br></div><div> Let's me try to answer each question as following. And let's assume that you have n landmarks in 3D.</div><div><br></div><div>1.) You want the TPS parameters</div><div> The TPS parameters in ITK implementation include </div>
<div> -an nx3 D matrix represents the TPS warp coefficients </div><div> -a 3x3 A matrix represents the affine coefficients (rotation, scale, etc)</div><div> -a 3x1 B matrix represents the translation</div><div> The matrix dimensions might be transposed. You cannot save these parameters using the current implementation in ITK but you can make a class deriving</div>
<div>from itk::<span style="background-color:rgb(255,255,255);color:rgb(34,34,34);font-family:arial,sans-serif;font-size:13px">ThinPlateSplineKernelTransform and add a function to save these parameters</span></div><div><span style="background-color:rgb(255,255,255);color:rgb(34,34,34);font-family:arial,sans-serif;font-size:13px"><br>
</span></div><div><span style="background-color:rgb(255,255,255);color:rgb(34,34,34);font-family:arial,sans-serif;font-size:13px">2) To compute points after transformation.</span></div><div><font color="#222222" face="arial, sans-serif"> You can set a for loop iterating through all of your points and call the function TransformPoints() to transform these points like you did previously.</font></div>
<div><font color="#222222" face="arial, sans-serif"><br></font></div><div><font color="#222222" face="arial, sans-serif"> </font><font face="arial, sans-serif" color="#000099">// to save points after transform</font></div>
<div><div><font color="#000099"> PointSetType::Pointer transformed_points = PointSetType::New();</font></div><div><font color="#000099"> PointSetType::PointsContainer::Pointer transformed_points_container = transformed_points->GetPoints(); </font></div>
<div><font color="#000099"> TPSTransformType::OutputPointType trackerPointNewPosition;</font></div><div><font color="#000099"><br></font></div><div><font color="#000099"> for(unsigned int i = 0; i < sourceLandMarks->GetNumberOfPoints(); i++){</font></div>
<div><font color="#000099"> PointType src_pnt = sourceLandMarks->GetPoint(i);</font></div><div><font color="#000099"> trackerPointNewPosition = tps->TransformPoint(trackerPoint);</font></div><div><font color="#000099"> transformed_points_container->InsertElement(i,trackerPointNewPosition); </font></div>
<div><font color="#000099"> }</font></div></div><div><font color="#222222" face="arial, sans-serif"><br></font></div><div><font color="#222222" face="arial, sans-serif"> </font></div><div><font color="#222222" face="arial, sans-serif">3.) To compute the registration error (RE), you know that the TPS transform in ITK is mapping from source->target.</font><font color="#222222" face="arial, sans-serif"> I</font>f your tracker points are the points defined in the SOURCE coordinate system, <span style="color:rgb(34,34,34);font-family:arial,sans-serif">you can </span><span style="color:rgb(34,34,34);font-family:arial,sans-serif">only </span><span style="color:rgb(34,34,34);font-family:arial,sans-serif">compute </span>|T(Point _tracker) - Point_image |. We cannot compute the inverse of the TPS transform. However, if you converse the TPS transform to a displacement field, there are classes in ITK to inverse the displacement field and you can compute |invT(Point_image) - Point_tracker|. (</div>
<div>theoretically these should give you approximately the same result. I suggest that you should just compute <font color="#222222" face="arial, sans-serif"> </font>|T(Point _tracker) - Point_image | )</div><div><br></div>
<div> You can add a few more lines in the previous for loop to compute the RE. Maybe something like this:</div><div><div> </div><div> <font color="#000099"> PointSetType::Pointer transformed_points = PointSetType::New();</font></div>
<div><font color="#000099"> PointSetType::PointsContainer::Pointer transformed_points_container = transformed_points->GetPoints();</font></div><div><font color="#009900"> std::vector<double> regis_errors;</font></div>
<div><font color="#000099"><br></font></div><div><font color="#000099"> for(unsigned int i = 0; i < sourceLandMarks->GetNumberOfPoints(); i++){</font></div><div><font color="#000099"> PointType src_pnt = sourceLandMarks->GetPoint(i);</font></div>
<div><font color="#000099"> TPSTransformType::OutputPointType trackerPointNewPosition;</font></div><div><font color="#000099"> trackerPointNewPosition = tps->TransformPoint(trackerPoint);</font></div><div><font color="#000099"> transformed_points_container->InsertElement(i,trackerPointNewPosition); </font></div>
<div><font color="#000099"><br></font></div><div><font color="#009900"> PointType target_pnt = targetLandMarks->GetPoint(i);</font></div><div><font color="#009900"> regis_errors.push_back(trackerPointNewPosition.EuclideanDistanceTo(target_pnt));</font></div>
<div><font color="#000099"> }</font></div></div><div><br></div><div>Hope this help.</div><div><br></div><div>Regards,</div><div><br></div><div>Ja</div><div><br></div><div><div class="gmail_quote">On Sat, Aug 11, 2012 at 11:57 AM, Agata Krasoń <span dir="ltr"><<a href="mailto:agatakrason@gmail.com" target="_blank">agatakrason@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi
Sureerat ,<div><br></div><div>Thanks for Your reply.</div><div>My problem is :</div><div>I have set of 5 sources points (from tracker - Polaris) and set of 5 targets points from (image - dicom).</div><div>I put in all in tab[30] :</div>
<div><br></div><div><div>double tab[30] = </div><div><span style="white-space:pre-wrap"> </span>{</div><div><span style="white-space:pre-wrap"> </span> // POLARIS</div><div class="im"><div><span style="white-space:pre-wrap"> </span>-81.29,-31.07,<span style="white-space:pre-wrap"> </span>-770.58, </div>
<div> -83.11,<span style="white-space:pre-wrap"> </span>-21.26,<span style="white-space:pre-wrap"> </span>-822.64,</div><div> -93.45,<span style="white-space:pre-wrap"> </span>-32.44,<span style="white-space:pre-wrap"> </span>-858.72,</div>
<div> -68.08,<span style="white-space:pre-wrap"> </span>-126.89,-813.07,</div><div> -61.04,<span style="white-space:pre-wrap"> </span>75.74,<span style="white-space:pre-wrap"> </span>-808.36,</div>
<div><br></div></div><div><span style="white-space:pre-wrap"> </span> // DICOM</div><div class="im"><div> 140.6,<span style="white-space:pre-wrap"> </span>230.7,<span style="white-space:pre-wrap"> </span>-30.5,</div>
<div> 140.2,<span style="white-space:pre-wrap"> </span>231.7,<span style="white-space:pre-wrap"> </span>-71.1,</div><div> 144.8,<span style="white-space:pre-wrap"> </span>235.9,<span style="white-space:pre-wrap"> </span>-116.1,</div>
<div> 45.8,<span style="white-space:pre-wrap"> </span>220.2,<span style="white-space:pre-wrap"> </span>-66.7,</div><div> 231.6,<span style="white-space:pre-wrap"> </span>211.3,<span style="white-space:pre-wrap"> </span>-66.1</div>
<div><span style="white-space:pre-wrap"> </span></div><div><span style="white-space:pre-wrap"> </span>};</div><div><span style="white-space:pre-wrap"> </span></div></div></div>
<div>I want to receive transform (matrix of TPS with parameters) or directly computed points after transformation. </div><div>After TPS registration I want to check errors - difference between |T(Point _tracker) - Point_image | & |T(Point_image) - Point_tracker|.</div>
<div>Difference between measured point and computed point. </div><div><br></div><div>Yes, I am trying to compute these points here :</div><div class="im"><div><br></div><div><div>TPSTransformType::OutputPointType trackerPointNewPosition;</div>
<div>trackerPointNewPosition = tps->TransformPoint(trackerPoint);</div><div>std::cout<<"trackerPointNewPosition: " <<trackerPointNewPosition<<std::endl;</div></div><div><br></div></div><div>
But I compute here only one point. When I choose tracker point I received only one point on image.</div>
<div>I think I receive the same values of computed image point as measured point (231.6,<span style="white-space:pre-wrap"> </span>211.3,<span style="white-space:pre-wrap"> </span>-66.1).</div>
<div>I attach a file with all my code.</div><div><br></div><div>I would appreciate for any help please. </div><div><br></div><div><br></div><div>Best,</div><div>Agatte</div><div class="HOEnZb"><div class="h5"><div><br></div>
<div><br></div><div><br></div><div>
<br><div class="gmail_quote">2012/8/11 Sureerat Reaungamornrat <span dir="ltr"><<a href="mailto:sureerat.r@gmail.com" target="_blank">sureerat.r@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi,<div><br></div><div> I do not quite understand your question. I think you want to know how to compute the tracker points after registration </div><div>using TPS. If this is your question, I think you have already solved it with this line of your code.</div>
<div><div><br><br><div> TPSTransformType::OutputPointType trackerPointNewPosition;</div><div> trackerPointNewPosition = tps->TransformPoint(trackerPoint);</div><div> std::cout<<"trackerPointNewPosition: " <<trackerPointNewPosition<<std::endl;</div>
<div><br></div></div><div> You can create a new PointSetType object to keep your tracker points after TPS transform.</div><div>If this is not your question, could you make your question a little clearer and maybe more specific?</div>
<div><br></div><div><br></div><div>Regards,</div><div><br></div><div>Ja </div><div><div><br><br><div class="gmail_quote">On Sat, Aug 11, 2012 at 7:05 AM, agatte <span dir="ltr"><<a href="mailto:wiatrak11@poczta.onet.pl" target="_blank">wiatrak11@poczta.onet.pl</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
Hi,<br>
<br>
I have a question about TPS transformation in ITK.<br>
I have set of 5 sources and target points(landmarks).<br>
I want to compute TPS transformation.<br>
I used ThinPlateSpline example for it.<br>
I want to compute new positions of point after transformation.<br>
What I can compute new position of each landmark after transformation ?<br>
I would be appreciate for any help please.<br>
<br>
Here I attach code :<br>
<br>
int main()<br>
{<br>
<br>
double tab[30] =<br>
{<br>
// sources landmarks from tracker<br>
-81.29,-31.07, -770.58,<br>
-83.11, -21.26, -822.64,<br>
-93.45, -32.44, -858.72,<br>
-68.08, -126.89,-813.07,<br>
-61.04, 75.74, -808.36,<br>
<br>
// target landmarks from image - dicom<br>
140.6, 230.7, -30.5,<br>
140.2, 231.7, -71.1,<br>
144.8, 235.9, -116.1,<br>
45.8, 220.2, -66.7,<br>
231.6, 211.3, -66.1<br>
<br>
};<br>
<br>
<br>
const unsigned int Dimension = 3;<br>
typedef unsigned char PixelType;<br>
typedef double CoordinateRepType;<br>
typedef itk::ThinPlateSplineKernelTransform<<br>
CoordinateRepType,Dimension> TPSTransformType;<br>
typedef itk::Point< CoordinateRepType,Dimension > PointType;<br>
typedef std::vector< PointType > PointArrayType;<br>
typedef TPSTransformType::PointSetType PointSetType;<br>
typedef PointSetType::Pointer PointSetPointer;<br>
typedef PointSetType::PointIdentifier PointIdType;<br>
<br>
<br>
<br>
<br>
// Landmarks correspondances may be associated with the<br>
SplineKernelTransforms<br>
// via Point Set containers. Let us define containers for the landmarks.<br>
PointSetType::Pointer sourceLandMarks = PointSetType::New();<br>
PointSetType::Pointer targetLandMarks = PointSetType::New();<br>
PointType trackerPoint; PointType imagePoint;<br>
PointSetType::PointsContainer::Pointer sourceLandMarkContainer =<br>
sourceLandMarks->GetPoints();<br>
PointSetType::PointsContainer::Pointer targetLandMarkContainer =<br>
targetLandMarks->GetPoints();<br>
<br>
<br>
<br>
<br>
// 1 Landmark<br>
trackerPoint[0] = tab[0];<br>
trackerPoint[1] = tab[1];<br>
trackerPoint[2] = tab[2];<br>
imagePoint[0] = tab[15];<br>
imagePoint[1] = tab[16];<br>
imagePoint[2] = tab[17];<br>
sourceLandMarkContainer->InsertElement( 0,trackerPoint);<br>
targetLandMarkContainer->InsertElement(0,imagePoint);<br>
<br>
// 2 Landmark<br>
trackerPoint[0] = tab[3];<br>
trackerPoint[1] = tab[4];<br>
trackerPoint[2] = tab[5];<br>
imagePoint[0] = tab[18];<br>
imagePoint[1] = tab[19];<br>
imagePoint[2] = tab[20];<br>
sourceLandMarkContainer->InsertElement(1,trackerPoint);<br>
targetLandMarkContainer->InsertElement(1,imagePoint);<br>
<br>
// 3 Landmark<br>
trackerPoint[0] = tab[6];<br>
trackerPoint[1] = tab[7];<br>
trackerPoint[2] = tab[8];<br>
imagePoint[0] = tab[21];<br>
imagePoint[1] = tab[22];<br>
imagePoint[2] = tab[23];<br>
sourceLandMarkContainer->InsertElement( 2,trackerPoint);<br>
targetLandMarkContainer->InsertElement(2,imagePoint);<br>
<br>
// 4 Landmark<br>
trackerPoint[0] = tab[9];<br>
trackerPoint[1] = tab[10];<br>
trackerPoint[2] = tab[11];<br>
imagePoint[0] = tab[24];<br>
imagePoint[1] = tab[25];<br>
imagePoint[2] = tab[26];<br>
sourceLandMarkContainer->InsertElement( 3,trackerPoint);<br>
targetLandMarkContainer->InsertElement(3,imagePoint);<br>
<br>
// 5 Landmark<br>
trackerPoint[0] = tab[12];<br>
trackerPoint[1] = tab[13];<br>
trackerPoint[2] = tab[14];<br>
imagePoint[0] = tab[27];<br>
imagePoint[1] = tab[28];<br>
imagePoint[2] = tab[29];<br>
sourceLandMarkContainer->InsertElement(4,trackerPoint);<br>
targetLandMarkContainer->InsertElement(4,imagePoint);<br>
<br>
TPSTransformType::Pointer tps = TPSTransformType::New();<br>
tps->SetSourceLandmarks(sourceLandMarks);<br>
tps->SetTargetLandmarks(targetLandMarks);<br>
// ComputeMAtrix<br>
tps->ComputeWMatrix();<br>
<br>
TPSTransformType::OutputPointType trackerPointNewPosition;<br>
trackerPointNewPosition = tps->TransformPoint(trackerPoint);<br>
std::cout<<"trackerPointNewPosition:<br>
"<<trackerPointNewPosition<<std::endl;<br>
<br>
// TPSTransformType::OutputVectorType pointsAfterTr;<br>
// pointsAfterTr = tps->TransformVector(sourceLandMarks);<br>
<br>
<br>
// save transformation<br>
typedef itk::TransformFileWriter TransformWriterType;<br>
TransformWriterType::Pointer transformWriter = TransformWriterType::New();<br>
transformWriter->AddTransform( tps);<br>
transformWriter->SetInput(tps);<br>
transformWriter->SetFileName("tps.txt");<br>
transformWriter->Update();<br>
<br>
<br>
<br>
return EXIT_SUCCESS;<br>
<br>
<br>
}<br>
<br>
<br>
<br>
<br>
--<br>
View this message in context: <a href="http://old.nabble.com/TPS--compute-new-position-of-points--tp34285325p34285325.html" target="_blank">http://old.nabble.com/TPS--compute-new-position-of-points--tp34285325p34285325.html</a><br>
Sent from the ITK - Users mailing list archive at Nabble.com.<br>
<br>
_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
</blockquote></div><br></div></div></div>
<br>_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
<br></blockquote></div><br></div>
</div></div></blockquote></div><br></div>