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