[vtkusers] vtkImageData point extraction
David Gobbi
david.gobbi at gmail.com
Sat Jan 1 18:02:39 EST 2011
Hi Kouroush,
It sounds to me like you want to flatten a structured grid, i.e. throw
away the geometry and just keep the scalars. I don't think there is
any filter for doing this, but you can try something like this:
1) Get the scalars by calling GetPointData()->GetScalars(), and
DeepCopy the scalars into a new array.
2) Generate a new vtkImageData from scratch using the scalars.
But for doing curved reformats like this, I would probably use a
different technique altogether. I'd write a new vtkWarpTransform
subclass that does the plane-to-curve transform that I need, and then
use vtkImageReslice to directly perform the curve-to-plane
interpolation. Either that, or I would create a vtkGridTransform that
maps from the planes to the desired curves.
David
On Sat, Jan 1, 2011 at 10:27 PM, Jafari, Kourosh <kjafari at rad.hfh.edu> wrote:
> Thanks David!
>
> Let me give you an example that I am currently working on:
>
> Generating panoramic dental images from 3D CT images and displaying them. Here are the main steps, in simple words, to create a panoramic image:
>
> 1. Cut the 3D CT image through parallel and equi-spaced planes (non-flat planes, e.g., semi-cylindrical) that pass through all teeth. The sampled points generate a volume (of course all the planes have similar and corresponding number of points).
>
> 2. Calculate the projection of the above generated volume in the third direction (the direction of increasing plane number). This is simply the average of all the above planes). The outcome is a set of points on a plane that are the intensities of the panoramic image. (here I just applied simple averaging, but I could apply a different formula instead of just averaging)
>
> 3. Display the panoramic image.
>
> I think now you know why I try to use vtkStructuredGrid to probe a vtkImageData and what I'd like to display.
>
> Thanks a lot!
> Kourosh
>
> ________________________________________
> From: David Gobbi [david.gobbi at gmail.com]
> Sent: Friday, December 31, 2010 1:37 PM
> To: Jafari, Kourosh
> Cc: vtkusers at vtk.org
> Subject: Re: [vtkusers] vtkImageData point extraction
>
> Hi Kourosh,
>
> The GetImageDataOutput() method does not do conversion. It will
> return NULL unless the data is already image data.
>
> I'm not sure what you mean when you say "display the output as an
> image". The output is a set of points. Do you want draw small circle
> in the image at each point location? Or draw a polygon in the image
> for each polygon in the data set?
>
> Polydata-to-image conversion is non-specific. There are several ways
> to do it, depending on what a person wants to achieve.
>
> David
>
>
> On Fri, Dec 31, 2010 at 6:04 PM, Jafari, Kourosh <kjafari at rad.hfh.edu> wrote:
>> Hi David,
>>
>> You were right in your original e-mail about Source and Input for vtkProbeFilter. The error in my code was coming from displaying the output. To summarize, as you said, this is how we should use vtkProbeFilter:
>>
>> probeFilter->SetSourceConnection(<data-to-be-interpolated goes here>)
>> probeFilter->SetInputConnection(<point input goes here>)
>>
>> In my code, the source is an image with vtkImageData type and the points are in a vtkStructuredGrid. In the following test code, I have manually generated an image since for some reason the code did not work with my loaded image, which still I have to figure out why. I have now problem with displaying the output of the probe. Ideally, I want the output to be in vtkImageData type in order to easily display it. But it is fine if I display it using any other method. Unfortunately, probe->GetImageDataOutput() whose output is supposed to be of type vtkImageData did not work. I have tried different methods to display the output but they all failed. I can, however, print the interpolated values. The probe->GetOutput() is of type vtkDataSet. How do I display it as an image?
>>
>> Thanks a lot!
>> Kourosh
>>
>>
>> P.S.
>>
>> The test code:
>>
>> // Manually create an image of type vtkImageData
>> vtkImageData *image2 = vtkImageData::New();
>> image2->SetScalarTypeToDouble();
>> image2->SetNumberOfScalarComponents(1);
>> image2->SetExtent(0, 9, 0, 9, 0, 9);
>> unsigned int GridSize = 10;
>> for ( unsigned int x = 0; x < GridSize; x++ )
>> {
>> for ( unsigned int y = 0; y < GridSize; y++ )
>> {
>> for ( unsigned int z = 0; z < GridSize; z++ )
>> {
>> double val = vtkMath::Random(0.0, 1.0)*100;
>> image2->SetScalarComponentFromDouble(x,y,z,0,val);
>> }
>> }
>> }
>>
>> // Create a structured grid and insert the points at which I would like to interpolate the image.
>> int i, j, k, kOffset, jOffset, offset;
>> float x[3], v[3], rMin=0, rMax=3.0, deltaRad, deltaZ;
>> float radius, theta;
>> static int dims[3]={10,10,1};
>> vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
>> sgrid->SetDimensions(dims);
>> vtkPoints *points = vtkPoints::New();
>> points->Allocate(dims[0]*dims[1]*dims[2]);
>> deltaZ = 2.0 / (dims[2]-1);
>> deltaRad = (rMax-rMin) / (dims[1]-1);
>> v[2]=0.0;
>> for ( k=0; k<dims[2]; k++)
>> {
>> x[2] = (-1.0 + k*deltaZ)*10;
>> kOffset = k * dims[0] * dims[1];
>> for (j=0; j<dims[1]; j++)
>> {
>> radius = rMin + j*deltaRad;
>> jOffset = j * dims[0];
>> for (i=0; i<dims[0]; i++)
>> {
>> x[0]=i/2;x[1]=j/2;x[2]=k/2;
>> offset = i + jOffset + kOffset;
>> points->InsertPoint(offset,x);
>> }
>> }
>> }
>> sgrid->SetPoints(points);
>> points->Delete();
>>
>> // Interpolate the image at points specified by the structured grid using vtkProbeFilter
>> vtkProbeFilter* probe = vtkProbeFilter::New();
>> probe->SetInput(sgrid);
>> probe->SetSource(image2);
>> probe->SpatialMatchOn();
>> probe->Update();
>>
>> // I can print the interpolated values:
>> vtkDataArray* data = probe->GetOutput()->GetPointData()->GetScalars();
>> vtkDoubleArray* doubleData = vtkDoubleArray::SafeDownCast (data);
>> cout << doubleData->GetNumberOfTuples()<< endl;
>> for(int i = 0; i < doubleData->GetNumberOfTuples(); i++)
>> {
>> double val = doubleData->GetValue(i);
>> cout << "Interpolation using ProbeFilter ";
>> cout << "doubleData->GetValue(" << i << "): " << val << endl;
>> }
>>
>>
>> ________________________________________
>> From: David Gobbi [david.gobbi at gmail.com]
>> Sent: Tuesday, December 28, 2010 12:53 PM
>> To: Jafari, Kourosh
>> Cc: vtkusers at vtk.org
>> Subject: Re: [vtkusers] vtkImageData point extraction
>>
>> Hi Karoush,
>>
>> You, it looks like I did have the Source and Input switched in my last email.
>> The "outside the whole extent" error is a bit surprising to me, because the
>> vtkProbeFilter is supposed to automatically set the UpdateExtent to be
>> exactly the same as the WholeExtent. So this error might be coming from
>> some earlier filter in the pipeline.
>>
>> One thing you can check is that GetBounds() on your image data returns a
>> bounds that is at least as large as GetBounds() on your structured grid.
>> I'm just making a guess here, but the Spacing and Origin of your image
>> data might not be correct. You should print them out (after calling Update,
>> of course) just to be sure.
>>
>> David
>>
>>
>> On Tue, Dec 28, 2010 at 9:31 AM, Jafari, Kourosh <kjafari at rad.hfh.edu> wrote:
>>> Thanks David,
>>>
>>> I had actually set up the points of the structured grid. I had replaced it with "..." in the e-mail to make it brief. I applied your suggestions but still it crashes. Following is the actual code (the structured grid here is just an example to test the code). Are you sure SetInput is for points in which we want to resample and SetSource is the volume? Actually, when I exchange their inputs the program does not crash but give the error message:
>>>
>>> ERROR: In ..\..\vtk-5.6.0\Filtering\vtkStreamingDemandDrivenPipeline.cxx, line 879
>>> vtkStreamingDemandDrivenPipeline (03CF52C8): The update extent specified in the information for output port 0 on algorithm vtkTrivialProducer(0259FD08) is 0 249 0 249 0 249, which is outside the whole extent 0 99 0 99 0 99.
>>>
>>> Thanks a lot,
>>> Kourosh
>>>
>>> ----------------------------------------------------------
>>> The actual code that causes the program to crash:
>>>
>>> vtkProbeFilter* probe = vtkProbeFilter::New();
>>> // Create the structured grid
>>> int i, j, k, kOffset, jOffset, offset;
>>> float x[3], v[3], rMin=0, rMax=3.0, deltaRad, deltaZ;
>>> float radius, theta;
>>> static int dims[3]={100,100,100};
>>> vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
>>> sgrid->SetDimensions(dims);
>>> vtkFloatArray *vectors = vtkFloatArray::New();
>>> vectors->SetNumberOfComponents(3);
>>> vectors->SetNumberOfTuples(dims[0]*dims[1]*dims[2]);
>>> vtkPoints *points = vtkPoints::New();
>>> points->Allocate(dims[0]*dims[1]*dims[2]);
>>> deltaZ = 2.0 / (dims[2]-1);
>>> deltaRad = (rMax-rMin) / (dims[1]-1);
>>> v[2]=0.0;
>>> for ( k=0; k<dims[2]; k++)
>>> {
>>> x[2] = (-1.0 + k*deltaZ)*10;
>>> kOffset = k * dims[0] * dims[1];
>>> for (j=0; j<dims[1]; j++)
>>> {
>>> radius = rMin + j*deltaRad;
>>> jOffset = j * dims[0];
>>> for (i=0; i<dims[0]; i++)
>>> {
>>> v[0] = -x[1];
>>> v[1] = x[0];
>>> v[2] = 2;
>>> x[0]=i;x[1]=j;x[2]=k;
>>> offset = i + jOffset + kOffset;
>>> points->InsertPoint(offset,x);
>>> vectors->InsertTuple(offset,v);
>>> }
>>> }
>>> }
>>> sgrid->SetPoints(points);
>>> points->Delete();
>>> sgrid->GetPointData()->SetVectors(vectors);
>>> vectors->Delete();
>>>
>>> probe->SetInput(sgrid);
>>> probe->SetSource(connector->GetOutput());
>>> probe->SpatialMatchOff();
>>> probe->Update();
>>> vtkImageData *imageP = vtkImageData::New();
>>> imageP->DeepCopy(probe->GetImageDataOutput());
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> ________________________________________
>>> From: David Gobbi [david.gobbi at gmail.com]
>>> Sent: Tuesday, December 28, 2010 2:45 AM
>>> To: Jafari, Kourosh
>>> Cc: vtkusers at vtk.org
>>> Subject: Re: [vtkusers] vtkImageData point extraction
>>>
>>> Hi Kourosh,
>>>
>>> Don't use SpatialMatchOn() until you first have it working properly
>>> with SpatialMatchOff(). Also, don't forget to do probe->Update()
>>> before you use the probe filter's output.
>>>
>>> If sgrid->SetDimensions(dims) the only thing that you did to set up
>>> the structured grid, then that isn't enough. The structured grid
>>> needs you to set the points (I've never used vtkStructuredGrid,
>>> myself, so I can't provide an example, but I'm sure that it takes a
>>> few more steps to create the data).
>>>
>>> David
>>>
>>>
>>> On Sat, Dec 25, 2010 at 8:38 PM, Jafari, Kourosh <kjafari at rad.hfh.edu> wrote:
>>>> Hi David,
>>>>
>>>> I have similar a question. I have a 3D image volume and would like to sample it on a structured grid. Here is how I use it:
>>>>
>>>> static int dims[3]={100,100,100};
>>>> vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
>>>> sgrid->SetDimensions(dims);
>>>>
>>>> ....
>>>>
>>>> probe->SetInput(sgrid);
>>>> probe->SetSource(connector->GetOutput());
>>>> probe->SpatialMatchOn();
>>>>
>>>> vtkImageData *imageP = vtkImageData::New();
>>>> imageP->DeepCopy(probe->GetImageDataOutput());
>>>>
>>>>
>>>> The type of connector->GetOutput() is vtkImageData which is the volume I'd like to sample. I later display imageP. There is no compilation error but it crashes when I run it. If I exchange the inputs of SetInput and SetSource, it does not crash but it does not do anything either. What am I missing?
>>>>
>>>> Thanks!!
>>>> Kourosh
>>>>
>>>>
>>>> ________________________________________
>>>> From: vtkusers-bounces at vtk.org [vtkusers-bounces at vtk.org] On Behalf Of David Gobbi [david.gobbi at gmail.com]
>>>> Sent: Friday, December 24, 2010 4:51 PM
>>>> To: Maxwell, Thomas P. (GSFC-606.2)[SCIENCE APPLICATIONS INTL CORP]
>>>> Cc: vtkusers at vtk.org
>>>> Subject: Re: [vtkusers] vtkImageData point extraction
>>>>
>>>> Hi Tom,
>>>>
>>>> I think you want vtkProbeFilter. It takes two inputs, one input for
>>>> volumetric data and another input for a set of points whose scalars
>>>> will be interpolated from the first input.
>>>>
>>>> probeFilter->SetSourceConnection(<data-to-be-interpolated goes here>)
>>>> probeFilter->SetInputConnection(<point input goes here>)
>>>>
>>>> Is this what you were looking for?
>>>>
>>>> David
>>>>
>>>>
>>>> On Thu, Dec 23, 2010 at 9:43 AM, Maxwell, Thomas P.
>>>> (GSFC-606.2)[SCIENCE APPLICATIONS INTL CORP] <thomas.maxwell at nasa.gov>
>>>> wrote:
>>>>>
>>>>> Is there a VTK filter that will take as input a vtkImageData and a list of
>>>>> (unstructured) points and then compute (as output) the image values at the
>>>>> specified point locations?
>>>>>
>>>>> Thanks,
>>>>> Tom
>>>> _______________________________________________
>>>> 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://www.vtk.org/mailman/listinfo/vtkusers
More information about the vtkusers
mailing list