[vtkusers] Unstructured grid to structured grid?
Chris Marsh
chris.marsh at usask.ca
Tue Sep 16 18:49:56 EDT 2014
Aha! Excellent, thanks. This all seems to work, mostly.
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
The code that converts my triangulation into the vtk format is shared
between the mesh and plane code, so it should 'work' given the mesh is
right.
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!)
Thanks for all the help
On Tue, Sep 16, 2014 at 2:12 PM, David Gobbi <david.gobbi at gmail.com> wrote:
> You also have to call gridPoints->Update(), because you are getting
> its output directly, rather than connecting to its output port.
>
> On Tue, Sep 16, 2014 at 2:07 PM, Chris Marsh <chris.marsh at usask.ca> wrote:
> > No dice. With the examples code, I would get a blank plane uniform plane
> > with the data name 'elevation' as per the SetName() call. However, now, I
> > get no plane, and for coloring all I see (in paraview) is 'solid color'
> and
> > 'elevation' is no longer present. Which suggests to me something else is
> > perhaps wrong?
> >
> >
> >
> > On Tue, Sep 16, 2014 at 2:03 PM, David Gobbi <david.gobbi at gmail.com>
> wrote:
> >>
> >> Oops, the second "SetPoint1" should be "SetPoint2".
> >>
> >> On Tue, Sep 16, 2014 at 2:02 PM, David Gobbi <david.gobbi at gmail.com>
> >> wrote:
> >> > The default plane size is 1x1, so you might need this:
> >> >
> >> > gridPoints->SetResolution(gridSize, gridSize);
> >> > gridPoints->SetOrigin(0, 0, 0);
> >> > gridPoints->SetPoint1(gridSize, 0, 0);
> >> > gridPoints->SetPoint1(0, gridSize, 0);
> >> >
> >> > On Tue, Sep 16, 2014 at 1:49 PM, Chris Marsh <chris.marsh at usask.ca>
> >> > wrote:
> >> >> Hi David,
> >> >>
> >> >> To clarify: I want to take an unstructured grid (x,y) with face data
> >> >> (z) and
> >> >> make it a 2D structured grid. I have connectivity information
> >> >> (triangulation).
> >> >>
> >> >> Looking at vtkPlaneSource: looks like it should work/make my life
> >> >> easier
> >> >>
> >> >> However, even with the above linked demo with the slight
> modification,
> >> >> I
> >> >> cannot make this work. This is most likely a PEBKAC as a result of my
> >> >> lack
> >> >> of familiarity with VTK so I appreciate you taking the time to help.
> >> >> The
> >> >> minimum working example is below, based upon this code
> >> >>
> http://www.itk.org/Wiki/VTK/Examples/Cxx/PolyData/InterpolateMeshOnGrid
> >> >>
> >> >> vtkSmartPointer<vtkPoints> points =
> vtkSmartPointer<vtkPoints>::New();
> >> >> vtkSmartPointer<vtkFloatArray> data =
> >> >> vtkSmartPointer<vtkFloatArray>::New();
> >> >> data->SetName("elevation");
> >> >>
> >> >> unsigned int gridSize = 10;
> >> >> float maxHeight = 5;
> >> >> for ( unsigned int i = 0; i < 100; ++i )
> >> >> {
> >> >> double x = vtkMath::Random(0, gridSize);
> >> >> double y = vtkMath::Random(0, gridSize);
> >> >> double z = vtkMath::Random(0, maxHeight);
> >> >> points->InsertNextPoint ( x, y, 0);
> >> >> data->InsertNextValue(z);
> >> >> }
> >> >>
> >> >> // Add the grid points to a polydata object
> >> >> vtkSmartPointer<vtkPolyData> polydata =
> >> >> vtkSmartPointer<vtkPolyData>::New();
> >> >> polydata->SetPoints ( points );
> >> >> polydata->GetPointData()->SetScalars(data);
> >> >>
> >> >> // Triangulate the grid points. If you do not have a mesh (points
> >> >> // only), the output will not be interpolated!
> >> >> vtkSmartPointer<vtkDelaunay2D> tri =
> >> >> vtkSmartPointer<vtkDelaunay2D>::New();
> >> >> tri->SetInputData ( polydata );
> >> >> tri->Update();
> >> >>
> >> >> // Create a grid of points to interpolate over
> >> >> vtkSmartPointer<vtkPlaneSource> gridPoints =
> >> >> vtkSmartPointer<vtkPlaneSource>::New();
> >> >> gridPoints->SetResolution(gridSize, gridSize);
> >> >> gridPoints->SetOrigin(0, 0, 0);
> >> >>
> >> >> // Perform the interpolation
> >> >> vtkSmartPointer<vtkProbeFilter> probeFilter =
> >> >> vtkSmartPointer<vtkProbeFilter>::New();
> >> >> probeFilter->SetSourceConnection(tri->GetOutputPort());
> >> >> probeFilter->SetInputData(gridPoints->GetOutput()); //
> >> >> // Interpolate 'Source' at these points
> >> >> probeFilter->Update();
> >> >>
> >> >> vtkSmartPointer<vtkDelaunay2D> gridDelaunay =
> >> >> vtkSmartPointer<vtkDelaunay2D>::New();
> >> >> gridDelaunay->SetInputConnection ( probeFilter->GetOutputPort()
> );
> >> >>
> >> >> vtkSmartPointer<vtkXMLPolyDataWriter> gridWriter =
> >> >> vtkSmartPointer<vtkXMLPolyDataWriter>::New();
> >> >> gridWriter->SetFileName ( "gridSurface.vtp" );
> >> >> gridWriter->SetInputConnection ( gridDelaunay->GetOutputPort() );
> >> >> gridWriter->Write();
> >> >>
> >> >>
> >> >>
> >> >>
> >> >> On Tue, Sep 16, 2014 at 12:07 PM, David Gobbi <david.gobbi at gmail.com
> >
> >> >> wrote:
> >> >>>
> >> >>> If you want to create a polydata that is a regular grid, then
> >> >>> use vtkPlaneSource with SetResolution(gridSize, gridSize).
> >> >>>
> >> >>> If you want your regular grid to be a vtkImageData, then
> >> >>> use vtkImageGridSource with SetDataExtent(), SetDataSpacing(),
> >> >>> an SetDataOrigin() to create the needed dimensions.
> >> >>>
> >> >>> Note that if your point cloud is just points without any
> connectivity
> >> >>> between the points (i.e. without any cells), then you cannot use
> >> >>> vtkProbeFilter.
> >> >>>
> >> >>>
> >> >>>
> >> >>>
> >> >>> On Tue, Sep 16, 2014 at 11:52 AM, Chris Marsh <chris.marsh at usask.ca
> >
> >> >>> wrote:
> >> >>> > Mmh thanks. I think I'll stick with the probe filter.
> >> >>> >
> >> >>> > I've been able to get this example working
> >> >>> >
> >> >>> >
> http://www.itk.org/Wiki/VTK/Examples/Cxx/PolyData/InterpolateMeshOnGrid
> >> >>> > if I have a real-world pt cloud, I'm a bit unclear how I should
> make
> >> >>> > the
> >> >>> > grid.
> >> >>> >
> >> >>> > If xll and yll are the lower x and y coordinates of an axis
> aligned
> >> >>> > bounding
> >> >>> > box
> >> >>> >
> >> >>> > size_t gridSize = 100;
> >> >>> > for ( size_t x = 0; x < gridSize; x++ )
> >> >>> > {
> >> >>> > for ( size_t y = 0; y < gridSize; y++ )
> >> >>> > {
> >> >>> > gridPoints->InsertNextPoint ( x+xll, y+yll, 0);
> >> >>> > }
> >> >>> > }
> >> >>> >
> >> >>> > does not seem to produce any output.
> >> >>> >
> >> >>> > Any suggestions?
> >> >>> >
> >> >>> > On Mon, Sep 15, 2014 at 4:39 PM, David Gobbi <
> david.gobbi at gmail.com>
> >> >>> > wrote:
> >> >>> >>
> >> >>> >> Another filter that might work is vtkShepardMethod, which takes
> >> >>> >> unstructured points (e.g. an arbitrary set of samples) as input
> and
> >> >>> >> produces a vtkImageData as output.
> >> >>> >>
> >> >>> >> The difference between vtkProbeFilter and vtkShepardMethod is
> >> >>> >> that the former performs the interpolation based on the
> >> >>> >> connectivity
> >> >>> >> of the points (i.e. based on the cells), while the latter ignores
> >> >>> >> the
> >> >>> >> connectivity and does the interpolation purely based on distance.
> >> >>> >>
> >> >>> >> - David
> >> >>> >>
> >> >>> >> On Mon, Sep 15, 2014 at 4:23 PM, Chris Marsh <
> chris.marsh at usask.ca>
> >> >>> >> wrote:
> >> >>> >> > Thanks David, this looks perfect.
> >> >>> >> >
> >> >>> >> >
> >> >>> >> > On Mon, Sep 15, 2014 at 4:10 PM, David Gobbi
> >> >>> >> > <david.gobbi at gmail.com>
> >> >>> >> > wrote:
> >> >>> >> >>
> >> >>> >> >> Hi Chris,
> >> >>> >> >>
> >> >>> >> >> This sounds like a job for vtkProbeFilter. The probe filter
> >> >>> >> >> takes
> >> >>> >> >> two
> >> >>> >> >> inputs,
> >> >>> >> >> the "input" and the "source". The idea is that it resamples
> the
> >> >>> >> >> scalars
> >> >>> >> >> from
> >> >>> >> >> the "source" onto the geometry of the "input". So, in this
> >> >>> >> >> case,
> >> >>> >> >> the
> >> >>> >> >> "source"
> >> >>> >> >> is your unstructured grid, and the "input" is any vtkImageData
> >> >>> >> >> that
> >> >>> >> >> has
> >> >>> >> >> the
> >> >>> >> >> origin, spacing, and extent that you want to use for the
> >> >>> >> >> regridding.
> >> >>> >> >> The
> >> >>> >> >> output of vtkProbeFilter will be an image that has the same
> >> >>> >> >> geometry
> >> >>> >> >> as
> >> >>> >> >> the input, but that has pixel values that have been
> interpolated
> >> >>> >> >> from
> >> >>> >> >> your
> >> >>> >> >> unstructured grid. You should be able to find several
> examples
> >> >>> >> >> on
> >> >>> >> >> the
> >> >>> >> >> wiki.
> >> >>> >> >>
> >> >>> >> >> - David
> >> >>> >> >>
> >> >>> >> >>
> >> >>> >> >> On Mon, Sep 15, 2014 at 3:59 PM, chrism <chris.marsh at usask.ca
> >
> >> >>> >> >> wrote:
> >> >>> >> >> > Is it possible to (easily) convert an unstructured grid to a
> >> >>> >> >> > structured
> >> >>> >> >> > grid
> >> >>> >> >> > (Image data)?
> >> >>> >> >> >
> >> >>> >> >> > Is unstructured -> polydata -> image data the only way? If
> >> >>> >> >> > so,
> >> >>> >> >> > how
> >> >>> >> >> > example
> >> >>> >> >> >
> >> >>> >> >> >
> >> >>> >> >> >
> >> >>> >> >> > (
> http://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/DataSetSurfaceFilter
> >> >>> >> >> > ??)
> >> >>> >> >> >
> >> >>> >> >> > Ultimately, the unstructured grid represents a topography,
> and
> >> >>> >> >> > I
> >> >>> >> >> > am
> >> >>> >> >> > looking
> >> >>> >> >> > to 'rasterize' it for use in a GIS.
> >> >>> >> >
> >> >>> >> >
> >> >>> >
> >> >>> >
> >> >>
> >> >>
> >
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20140916/727f2859/attachment.html>
More information about the vtkusers
mailing list