[vtkusers] Unstructured grid to structured grid?

David Gobbi david.gobbi at gmail.com
Wed Sep 17 12:34:14 EDT 2014


Hi Chris,

The plane is already composed of cells, so you should not run it
through Delaunay.  If you want to convert its polygons (rectangles)
into triangles, you should do so with vtkTriangleFilter.

Delanay2D: triangulate unconnected points

TriangleFilter: triangulate cells (e.g. convert polys to triangles)

But there is no need to triangulate the plane.  It's original set of
rectangular polygons should be just fine for what you are doing.


As for the bad bounds you get from gridDelaunay, it looks like
you didn't call Update on the filter before getting the bounds from
its output.  The output of a filter is always invalid until you call
Update on the filter.

 - David


On Wed, Sep 17, 2014 at 10:20 AM, Chris Marsh <chris.marsh at usask.ca> wrote:
> Further to this:
>
> Tri bounds xmin:624714 xmax:629904 ymin:5.64502e+06 ymax:5.64877e+06
> gridpoints bounds xmin:624714 xmax:629904 ymin:5.64502e+06 ymax:5.64877e+06
> gridDelaunay bounds xmin:1 xmax:-1 ymin:1 ymax:-1
>
> It appears the final triangulation of the probed points is the problem
>
> vtkSmartPointer<vtkDelaunay2D> gridDelaunay =
> vtkSmartPointer<vtkDelaunay2D>::New();
>     gridDelaunay->SetInputConnection ( probeFilter->GetOutputPort() );
>
>     vtkSmartPointer<vtkXMLPolyDataWriter> gridWriter =
> vtkSmartPointer<vtkXMLPolyDataWriter>::New();
>     gridWriter->SetFileName ( file_name.c_str());
>     gridWriter->SetInputConnection ( gridDelaunay->GetOutputPort() );
>     gridWriter->Write();
>
> Is this final triangulation required? Is there a better way of doing this?
>
> Thanks, Chris
>
> On Tue, Sep 16, 2014 at 7:03 PM, Chris Marsh <chris.marsh at usask.ca> wrote:
>>
>> Ah I see resolution is just subdivisions.
>> http://www.vtk.org/doc/nightly/html/classvtkPlaneSource.html#acd65f04a453ba87b766048f0e7597710
>>
>> I have also confirmed bbox matches exactly the values that GetBounds()
>> reports.
>>
>> On Tue, Sep 16, 2014 at 6:16 PM, Chris Marsh <chris.marsh at usask.ca> wrote:
>>>
>>> Yes, I was trying to be optimistic :( Good call on the bounds:
>>> vtu export: xmin:624714 xmax:629904 ymin:5.64502e+06 ymax:5.64877e+06
>>> vtp export: xmin:1 xmax:-1 ymin:1 ymax:-1
>>>
>>> Presumably the problem. I'm unclear why it is doing this though?
>>>
>>> vtkSmartPointer<vtkPlaneSource> gridPoints =
>>> vtkSmartPointer<vtkPlaneSource>::New();
>>>
>>> int gridSize = 1000;
>>> gridPoints->SetResolution(gridSize, gridSize);
>>> gridPoints->SetOrigin(bbox.vertex(0).x(),  bbox.vertex(0).y(), 0);
>>> gridPoints->SetPoint1(bbox.vertex(2).x(), bbox.vertex(0).y(), 0);
>>> gridPoints->SetPoint2(bbox.vertex(0).x(), bbox.vertex(2).y(), 0);
>>> gridPoints->Update();
>>>
>>> where bbox is a bounding box (2D), index 0 is bottom left, indexing is
>>> increasing CCW. So 2 is upper right corner. As far as I can tell this is
>>> correct and corresponds to the vtu export numbers I had above.
>>>
>>> I feel like I perhaps don't understand resolution: I envision it as the
>>> size of each cell, however I don't think that is correct is it?
>>>
>>>
>>>
>>> On Tue, Sep 16, 2014 at 5:23 PM, David Gobbi <david.gobbi at gmail.com>
>>> wrote:
>>>>
>>>> On Tue, Sep 16, 2014 at 4:49 PM, Chris Marsh <chris.marsh at usask.ca>
>>>> wrote:
>>>>
>>>> > The interpolation doesn't look very good, even after increasing the
>>>> > resolution.  As well, it is flipped along the y-axis. Image is a
>>>> > mountain,
>>>> > variable is air temp (blue = cold, red = warm).
>>>> > See here: http://imgur.com/l1KpcS0
>>>>
>>>> Honestly, I can't see any correspondence at all between the mountain
>>>> and the interpolated data on that plane.  Are you sure the size and
>>>> position of the plane matches the size and position of your data?
>>>> An easy way to check is to call GetBounds() on both of them and
>>>> print out the bounds.  They should match.
>>>>
>>>> > Secondly: I need to iterate over the rows/cols of this interpolated
>>>> > grid in
>>>> > order to write it to an ArcGIS ascii file. What would be the best way
>>>> > to do
>>>> > this? The fact it is now in a delaunay triangulation somewhat doesn't
>>>> > help
>>>> > (now that I think about it!)
>>>>
>>>> The interpolated data can be retrieved via GetPointData()->GetScalars(),
>>>> but I'm not sure what ordering the vtkPlaneSource uses for its points.
>>>> You can print GetPoint(0) and GetPoint(1) to see whether the points
>>>> increase in the X direction or Y direction first.
>>>>
>>>> Also, watch for off-by-one mistakes.  SetResolution(1, 1) gives a plane
>>>> with 4 points (the corners).  So if you want a plane with 4x4 points,
>>>> you
>>>> would SetResolution(3, 3).
>>>>
>>>>  - David
>>>
>>>
>>
>


More information about the vtkusers mailing list