How to find in which triangle a point falls
Tim Hutton
T.Hutton at eastman.ucl.ac.uk
Thu Mar 2 09:34:36 EST 2000
Etienne wrote:
>Thanks for the reply, it looks just like what I need.
>
>You mention doing piecwise affine warps on images. I think what I am doing
>is basically the same thing. I have got a triangular mesh defining one
>projection of geographic data. I have got another corresponding mesh
>defining another projection. I want to convert a set of points from the one
>projection to the other using the two meshes. What it boils down to is
>finding in which triangle in the first mesh the point falls. Get the three
>node's dx and dy in relation to the second mesh. Now calculate the point's
>dx,dy in relation to the three node's dx and dy's.
>
>Now, my question is : if you have dx1,dy1 ; dx2,dy2 ; dx3,dy3 of the
>triangle's nodes, how do you calculate the dx,dy of the point that falls in
>the triangle. I'll probably figure out how to calculate the dx,dy with some
>convoluted trigonometry, but I doubt if that will be very efficient! Is
>there some "simple" formula (or even better - a method somewhere in VTK)
>that I can use to calculate that. I was thinking along the lines of using
>the inverses of the squares of the distances to the nodes from the point (or
>whatever :)
This is exactly what I'm doing. I've included the crucial bits of my code
below. Apologies for the state they're in, it's a bit of a mindless hack.
My initial triangulation comes from a Delaunay on the target points, the
source points correspond and so don't need to be triangulated. The
parameterized coordinates bit you talk about is in ComputeSourcePoint,
helpfully vtkCellLocator gives you the r and s values you need.
But wait! Before you dive into this code and try and use it, have you
considered using a different kind of warping? Transforming between
different geographical projections using piecewise affine warping (you're
right, it's the same) will of course introduce small errors as I'm sure
you're aware. (You will get away with it if the triangles are small
relative to the distortion.) One problem with piecewise affine is that it
can introduce kinks in previously straight lines, I'm thinking this might
look odd with your data. I think a better (and easier to code!) solution
for you might be to use thin-plate spline warping. Check out
vtkThinPlateSplineTransform in the nightlies. I'm not 100% exactly how to
use it on 2d data but I know it can be done, David Gobbi is the chap to ask.
Otherwise, is there absolutely no way to directly compute the transform
between the two projections, do you really need to use the meshes?
Cheers,
Tim.
----------------------------------------------------------------------------
--------
vtkDelaunay2D *mean_polys = vtkDelaunay2D::New();
mean_polys->SetInput(mean.GetPolyData(false));
mean_polys->SetOffset(2);
mean_polys->Update();
// sample across the mean shape, taking pixel values from the
corresponding place
// in the current shape
vtkCellLocator *locator = vtkCellLocator::New();
locator->SetDataSet(mean_polys->GetOutput());
locator->AutomaticOff();
locator->SetNumberOfCellsPerBucket( 5 );
locator->BuildLocator();
vtkPoints *current_points =
theDoc->m_current.GetPolyData(false)->GetPoints();
float t;
float pcoords[3]; // parameterised coordinates
int subid,cellid;
float a0[3],a1[3];
float input[3],output[3];
int x,y;
for(x=0;x<IMAGE_SIZE;x++)
{
for(y=0;y<IMAGE_SIZE;y++)
{
input[0]=x;
input[1]=y;
input[2]=0;
float p[3];
a0[0]=x;
a0[1]=y;
a0[2]=-0.1; // z coord below z=0 plane
a1[0]=x;
a1[1]=y;
a1[2]=0.1; // z coord above z=0 plane
int result =
locator->IntersectWithLine(a0,a1,0.1,t,p,pcoords,subid,cellid);
if(result==1) // is this point within any triangle?
{
ComputeSourcePoint(mean_polys->GetOutput(),current_points,cellid,pcoords,out
put);
// sample the source image
if(output[0]>=0 && output[0]<theDoc->m_image.GetX() &&
output[1]>=0 && theDoc->m_image.GetY()-output[1]<theDoc->m_image.GetY() )
{
//dc.SetPixel(output[0],output[1],RGB(0,255,0));
float intensity =
theDoc->m_image.getInterpolated(output[0],theDoc->m_image.GetY()-output[1]);
dc.SetPixel(x,y,RGB(intensity*255,intensity*255,intensity*255));
}
else dc.SetPixel(x,y,RGB(255,0,0)); // off the image
}
else
dc.SetPixel(x,y,RGB(0,0,0)); // outside the mean shape
}
}
// ...
}
void ComputeSourcePoint(vtkPolyData *mean,vtkPoints *source_points,int
cellid,float pcoords[3],float *output)
{
vtkIdList *ids = vtkIdList::New();
float v1[3],v2[3],v3[3];
mean->GetCellPoints(cellid,ids);
mean->GetPoints()->GetPoint(ids->GetId(0),v1);
mean->GetPoints()->GetPoint(ids->GetId(1),v2);
mean->GetPoints()->GetPoint(ids->GetId(2),v3);
C3dVector p0(v1[0],v1[1],v1[2]),p1(v2[0],v2[1],v2[2]),p2(v3[0],v3[1],v3[2]);
source_points->GetPoint(ids->GetId(0),v1);
source_points->GetPoint(ids->GetId(1),v2);
source_points->GetPoint(ids->GetId(2),v3);
C3dVector
p0d(v1[0],v1[1],v1[2]),p1d(v2[0],v2[1],v2[2]),p2d(v3[0],v3[1],v3[2]);
float r,s;
r=pcoords[0];
s=pcoords[1];
C3dVector out = p0d + (p1d-p0d)*r + (p2d-p0d)*s;
output[0]=out.x;
output[1]=out.y;
output[2]=0.0F;
ids->Delete();
}
>-----Original Message-----
>From: Tim Hutton [mailto:T.Hutton at eastman.ucl.ac.uk]
>Sent: Thursday, March 02, 2000 2:02 PM
>To: Etienne Labuschagne; vtkusers at public.kitware.com
>Subject: Re: How to find in which triangle a point falls
>
>
>At 11:58 02/03/00 +0200, you wrote:
>>Hi there
>>
>>I'm quite new to VTK and I thought the fastest way to decide if I must use
>>VTK or my own algorithm will be through the people who are already
>>experienced with VTK! Here is my question :
>>
>>I have a set of data with elements that basically look like this :
>>TriangleID, x1,y1,x2,y2,x3,y3
>>
>>The three coordinates are the three points of a triangle (if you couldn't
>>guess it from the description ;) and all the triangles make up an
>>interconnected mesh. (Much like a terrain mesh, but only 2D)
>>
>>If I have a point (x,y) I would like to know in which triangle it falls in
>>as short a time as possible. The total dataset I must eventually work
>>through (triangles and points falling in them) will be about 30GB so an
>>efficient algorithm is quite important for me.
>>
>>It seems as if VTK can help me with this, but it is such an all
>encompassing
>>library that I think it will be faster to implement my own algorithm that
>>learning to use the whole VTK! That is, unless you guys can help me by
>>suggesting which VTK classes will do the job best. For example, how do I
>>know if a point falls within a certain triangle - is there something in VTK
>>that can do that for me, or must I do it myself?
>>
>>I have played around a bit with VTK and have got a BASIC knowledge of it's
>>data structures, so I'm not completely in the dark (I have already used the
>>Delaunay2D class successfully to build a TIN).
>>
>>Any help would be greatly appreciated.
>>Etienne
>
>This is just one way to do it:
>
>vtkCellLocator will do exactly what you want, and I think quite efficiently
>using buckets. Put your triangles as VTK_TRIANGLE cells into a vtkPolyData,
>leaving all the z-coords as zero. Call vtkCellLocator::SetDataSet(my_polys)
>and then vtkCellLocator::BuildLocator(). Then for each sample point, call
>vtkCellLocator::IntersectWithLine(), passing it a line between +0.1 and
>-0.1 on the z-axis, say, to give it something to work with. One of the
>return parameters will be the id of the cell that was hit.
>
>It's a bit kludgy to get around the fact that vtkCellLocator is really for
>3d data but I use this method successfully. If anyone knows a 2d way then
>I'd be interested too.
>
>Actually worse than that, there used to be a bug (there still may be), that
>is all the z-coords are zero then vtkCellLocator doesn't work, so I perturb
>the first point by a tiny amount and then it works. This was still a
>problem in 2.4, may have been solved in 3.1
>
>I use this for doing piecewise affine warps of images using two
>corresponding triangulations.
>
>Hope this helps, vtk really can do what you want!
>
>Cheers,
>
>Tim.
>
>
>
>---------------------------------------------------------------------------
>Tim Hutton, Research Assistant Email: T.Hutton at eastman.ucl.ac.uk
>MINORI Project Eternal: T.Hutton at excite.co.uk
>Dental and Medical Informatics http://www.eastman.ucl.ac.uk/~dmi/MINORI
>
>Eastman Dental Institute, UCL Tel: [+44] (0207) 915 2344
>256 Gray's Inn Road, London WC1X 8LD, UK Fax: [+44] (0207) 915 2303
>---------------------------------------------------------------------------
>This email represents the views of the sender alone and must not be
>construed as representing the views of the Eastman Dental Institute. It may
>contain confidential information and may be protected by law as a legally
>privileged document and copyright work. Its content should not be disclosed
>and it should not be given or copied to anyone other than the person(s)
>named or referenced above. If you have received this email in error, please
>contact the sender.
>
>--------------------------------------------------------------------
>This is the private VTK discussion list. Please keep messages on-topic.
>Check the FAQ at: <http://public.kitware.com/cgi-bin/vtkfaq>
>To UNSUBSCRIBE, send message body containing "unsubscribe vtkusers" to
><majordomo at public.kitware.com>. For help, send message body containing
>"info vtkusers" to the same address.
>--------------------------------------------------------------------
>
>
--------------------------------------------------------------------
This is the private VTK discussion list. Please keep messages on-topic.
Check the FAQ at: <http://public.kitware.com/cgi-bin/vtkfaq>
To UNSUBSCRIBE, send message body containing "unsubscribe vtkusers" to
<majordomo at public.kitware.com>. For help, send message body containing
"info vtkusers" to the same address.
--------------------------------------------------------------------
More information about the vtkusers
mailing list