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