[vtkusers] Unstructured grid to structured grid?

Bill Lorensen bill.lorensen at gmail.com
Wed Sep 17 12:59:53 EDT 2014


Not to confuse things, but this example uses two different methods to
interpolate:
http://vtk.org/Wiki/VTK/Examples/Cxx/PolyData/InterpolateTerrain

Bill


On Wed, Sep 17, 2014 at 12:45 PM, Chris Marsh <chris.marsh at usask.ca> wrote:
> Hi David,
>
> Thanks for the reply, this is slowly starting to make sense. That call did
> in fact fix the bounds issues. However, the output still looks like it did
> in the above screenshot.
>
>>The plane is already composed of cells, so you should not run it
> through Delaunay.
>
> This is what perplexed me about the above linked example.
>
> What is the best way to write the probe results to a file then?
>
>
> Saving the after-probe'd grid does not produce the expected output.
> gridWriter->SetInputConnection ( gridPoints->GetOutputPort());
> gridWriter->Write();
>
> I feel like I am missing something simple. Thanks for the help.
>
> Chris Marsh
> PhD Student
> chrismarsh.ca
>
> 13 Kirk Hall
> University of Saskatchewan
>
> On Wed, Sep 17, 2014 at 10:34 AM, David Gobbi <david.gobbi at gmail.com> wrote:
>>
>> 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
>> >>>
>> >>>
>> >>
>> >
>
>
>
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/vtkusers
>



-- 
Unpaid intern in BillsBasement at noware dot com


More information about the vtkusers mailing list