Color PolyData based on VolumeData

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Color PolyData based on VolumeData

Radosław Furmaniak
Hi All,
I am trying to prepare a PolyData object from a series od DCM images and color it using some transfer function. I managed to create PolyData using vtkContourFilter.
I also have vtkVolume object with mapped colors, so I know that I have correct color transfer function and I am able to map intensities onto desired colors.
 What I am trying to do is achieve the same results colorwise with the PolyData from ContourFilter. However, I do not know how to tackle this problem. I have tried using ProbeFilter with PolyData as input and RGB ImageData as the source, but the output is the same as input.
Any suggestions how to approach this problem or if I am going the right way? I am fairly new to the VTK, so I do not know many of its possibilites yet. 
I'm working on Python 3.6 with VTK 8.1.1.
Any help is appreciated.

_______________________________________________
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

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
https://public.kitware.com/mailman/listinfo/vtkusers
Reply | Threaded
Open this post in threaded view
|

Re: Color PolyData based on VolumeData

Andaharoo
This post was updated on .
So, if I understand right, you want to color the scalar/whatever field onto
the poly? Just provide vtkPolyData with scalars and set the polyData's
mapper scalar visibility to on.

If you give it 1 component scalars you can give the mapper a color function
to map them. If you give it 3 component uchar it should just display them as
color (as long as scalar visibility is on). Like so:


// Create the scalars array
vtkDoubleArray* scalars = vtkDoubleArray::New();
scalars->SetName("Scalars");
// For every point in the poly
for (vtkIdType i = 0; i < outputPoly->GetNumberOfPoints(); i++)
{
        double pt[3];
        outputPoly->GetPoint(i, pt);

        scalars->InsertTuple1(i, trilinearSamplePoint(inputImage, pt[0], pt[1],
pt[2], 1, 0));
}
// Add the scalar data to the poly data
outputPoly->GetPointData()->AddArray(scalars);
outputPoly->GetPointData()->SetActiveScalars("Scalars");


This just goes through all the points in the polygon and gets the color from
the image. Getting the color at a point in an image can be done by simply
using the nearest color from the image (nearest neighbor) or better would be
to use trilinear interpolation of the 8 nearest colors which I used in my
example. Also note my code expects float image input:
static const int calcIndex(int x, int y, int z, int width, int height) {
return x + width * (y + height * z); }

static float trilinearSamplePoint(vtkImageData* imageData, float x, float y,
float z, int numComps, int comp)
{
        float* imagePtr = static_cast<float*>(imageData->GetScalarPointer());
        double* spacing = imageData->GetSpacing();
        double* origin = imageData->GetOrigin();
        int* dim = imageData->GetDimensions();

        // We assume point x, y, z is in the image
        int xi = (x - origin[0]) / spacing[0];
        int yi = (y - origin[1]) / spacing[1];
        int zi = (z - origin[2]) / spacing[2];

        // Get the intensities at the 8 nearest neighbors
        float i000 = imagePtr[calcIndex(xi, yi, zi, dim[0], dim[1]) * numComps +
comp];
        float i100 = imagePtr[calcIndex(xi + 1, yi, zi, dim[0], dim[1]) * numComps
+ comp];
        float i110 = imagePtr[calcIndex(xi + 1, yi + 1, zi, dim[0], dim[1] *
numComps + comp)];
        float i010 = imagePtr[calcIndex(xi, yi + 1, zi, dim[0], dim[1]) * numComps
+ comp];

        float i001 = imagePtr[calcIndex(xi, yi, zi + 1, dim[0], dim[1]) * numComps
+ comp];
        float i101 = imagePtr[calcIndex(xi + 1, yi, zi + 1, dim[0], dim[1]) *
numComps + comp];
        float i111 = imagePtr[calcIndex(xi + 1, yi + 1, zi + 1, dim[0], dim[1]) *
numComps + comp];
        float i011 = imagePtr[calcIndex(xi, yi + 1, zi + 1, dim[0], dim[1]) *
numComps + comp];

        // Get the fractional/unit distance from nearest neighbor 000
        float rx = xi * spacing[0] + origin[0]; // Position of node
        rx = (x - rx) / spacing[0]; // (Node - actual point position) / voxel width
        float ry = yi * spacing[1] + origin[1];
        ry = (y - ry) / spacing[1];
        float rz = zi * spacing[2] + origin[2];
        rz = (z - rz) / spacing[2];

        // Now we do the trilinear interpolation
        float ax = i000 + (i100 - i000) * rx;
        float bx = i010 + (i110 - i010) * rx;
        float cy = ax + (bx - ax) * ry;

        float dx = i001 + (i101 - i001) * rx;
        float ex = i011 + (i111 - i011) * rx;
        float fy = dx + (ex - dx) * ry;

        float gz = cy + (fy - cy) * rz;
        return gz;
}

I believe if you get the vtkCell from the data you can get VTK to do the interpolation for you as well but I've never tried it.



--
Sent from: http://vtk.1045678.n5.nabble.com/VTK-Users-f1224199.html
_______________________________________________
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

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
https://public.kitware.com/mailman/listinfo/vtkusers
Reply | Threaded
Open this post in threaded view
|

Re: Color PolyData based on VolumeData

Andras Lasso
You can copy voxel values from a volume to a polydata using vtkProbeFilter.

However, this won't help you if you obtained your polydata using a vtkContourFilter, because the contour filter generates isosurface, which means all extracted points have the same intensity value in the input volume. No matter what transfer function you apply, you'll always end up with a solid colored surface.

It is a common desire to replace volume rendering by displaying a colored surface, but this is not feasible (https://discourse.slicer.org/t/save-volume-rendering-as-stl-file/524/20).

If your goal is 3D printing then you may use voxel printing technique, which allows printing directly from volume rendering, without segmenting any surfaces (https://discourse.slicer.org/t/printing-volume-renderings-in-plastic/3017).

If your goal is volume-rendering-like display (for example, you want to show "volume rendering" in Unity) then you need to use an actual volume renderer.

Andras


-----Original Message-----
From: vtkusers <[hidden email]> On Behalf Of Andaharoo
Sent: Sunday, December 2, 2018 11:27 PM
To: [hidden email]
Subject: Re: [vtkusers] Color PolyData based on VolumeData

So, if I understand right, you want to color the scalar/whatever field onto the poly? Just provide vtkPolyData with scalars and set the polyData's mapper scalar visibility to on.

If you give it 1 component scalars you can give the mapper a color function to map them. If you give it 3 component uchar it should just display them as color (as long as scalar visibility is on). Like so:


// Create the scalars array
vtkDoubleArray* scalars = vtkDoubleArray::New();
scalars->SetName("Scalars");
// For every point in the poly
for (vtkIdType i = 0; i < outputPoly->GetNumberOfPoints(); i++) {
        double pt[3];
        outputPoly->GetPoint(i, pt);

        scalars->InsertTuple1(i, trilinearSamplePoint(inputImage, pt[0], pt[1], pt[2], 1, 0)); } // Add the scalar data to the poly data
outputPoly->GetPointData()->AddArray(scalars);
outputPoly->GetPointData()->SetActiveScalars("Scalars");


This just goes through all the points in the polygon and gets the color from the image. Getting the color at a point in an image can be done by simply using the nearest color from the image (nearest neighbor) or better would be to use trilinear interpolation of the 8 nearest colors which I used in my example. Also note my code expects float image input:
static const int calcIndex(int x, int y, int z, int width, int height) { return x + width * (y + height * z); }

static float trilinearSamplePoint(vtkImageData* imageData, float x, float y, float z, int numComps, int comp) {
        float* imagePtr = static_cast<float*>(imageData->GetScalarPointer());
        double* spacing = imageData->GetSpacing();
        double* origin = imageData->GetOrigin();
        int* dim = imageData->GetDimensions();

        // We assume point x, y, z is in the image
        int xi = (x - origin[0]) / spacing[0];
        int yi = (y - origin[1]) / spacing[1];
        int zi = (z - origin[2]) / spacing[2];

        // Get the intensities at the 8 nearest neighbors
        float i000 = imagePtr[calcIndex(xi, yi, zi, dim[0], dim[1]) * numComps + comp];
        float i100 = imagePtr[calcIndex(xi + 1, yi, zi, dim[0], dim[1]) * numComps
+ comp];
        float i110 = imagePtr[calcIndex(xi + 1, yi + 1, zi, dim[0], dim[1] * numComps + comp)];
        float i010 = imagePtr[calcIndex(xi, yi + 1, zi, dim[0], dim[1]) * numComps
+ comp];

        float i001 = imagePtr[calcIndex(xi, yi, zi + 1, dim[0], dim[1]) * numComps
+ comp];
        float i101 = imagePtr[calcIndex(xi + 1, yi, zi + 1, dim[0], dim[1]) * numComps + comp];
        float i111 = imagePtr[calcIndex(xi + 1, yi + 1, zi + 1, dim[0], dim[1]) * numComps + comp];
        float i011 = imagePtr[calcIndex(xi, yi + 1, zi + 1, dim[0], dim[1]) * numComps + comp];

        // Get the fractional/unit distance from nearest neighbor 000
        float rx = xi * spacing[0] + origin[0]; // Position of node
        rx = (x - rx) / spacing[0]; // (Node - actual point position) / voxel width
        float ry = yi * spacing[1] + origin[1];
        ry = (y - ry) / spacing[1];
        float rz = zi * spacing[2] + origin[2];
        rz = (z - rz) / spacing[2];

        // Now we do the trilinear interpolation
        float ax = i000 + (i100 - i000) * rx;
        float bx = i010 + (i110 - i010) * rx;
        float cy = ax + (bx - ax) * ry;

        float dx = i001 + (i101 - i001) * rx;
        float ex = i011 + (i111 - i011) * rx;
        float fy = dx + (ex - dx) * ry;

        float gz = cy + (fy - cy) * rz;
        return gz;
}




--
Sent from: https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fvtk.1045678.n5.nabble.com%2FVTK-Users-f1224199.html&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091371821&amp;sdata=iziee8ltwEAWna3RvM3dma7A3p94NWMHCfIwclRpR0I%3D&amp;reserved=0
_______________________________________________
Powered by https://na01.safelinks.protection.outlook.com/?url=www.kitware.com&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091371821&amp;sdata=6Nv%2FWdIIIaOq%2B3peZlXxuGdAwkOwzORn7AfWthXw5Xo%3D&amp;reserved=0

Visit other Kitware open-source projects at https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.kitware.com%2Fopensource%2Fopensource.html&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=oEIpua2iKsLpUdJaY%2BmQdOhCHfusGYUZDRd3Ru67Qoo%3D&amp;reserved=0

Please keep messages on-topic and check the VTK FAQ at: https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.vtk.org%2FWiki%2FVTK_FAQ&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=7h18pSmCa3LapcNQJsqqxU7kyWVyaQzKbPYZPMCbDiw%3D&amp;reserved=0

Search the list archives at: https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fmarkmail.org%2Fsearch%2F%3Fq%3Dvtkusers&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=luU4EgKZZ9nD1OOxN5EGfBpA5HnK%2FyuVoN6C7PE%2Bzhc%3D&amp;reserved=0

Follow this link to subscribe/unsubscribe:
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpublic.kitware.com%2Fmailman%2Flistinfo%2Fvtkusers&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=URr5vqvvkUuWmwaPIIGXXJzsmW2LFrYs7Z7jDcvt4ek%3D&amp;reserved=0
_______________________________________________
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

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
https://public.kitware.com/mailman/listinfo/vtkusers
Reply | Threaded
Open this post in threaded view
|

Re: Color PolyData based on VolumeData

Radosław Furmaniak
Thanks a lot! 
The solution from Andaharoo was what I was looking for. 

pon., 3 gru 2018 o 07:50 Andras Lasso <[hidden email]> napisał(a):
You can copy voxel values from a volume to a polydata using vtkProbeFilter.

However, this won't help you if you obtained your polydata using a vtkContourFilter, because the contour filter generates isosurface, which means all extracted points have the same intensity value in the input volume. No matter what transfer function you apply, you'll always end up with a solid colored surface.

It is a common desire to replace volume rendering by displaying a colored surface, but this is not feasible (https://discourse.slicer.org/t/save-volume-rendering-as-stl-file/524/20).

If your goal is 3D printing then you may use voxel printing technique, which allows printing directly from volume rendering, without segmenting any surfaces (https://discourse.slicer.org/t/printing-volume-renderings-in-plastic/3017).

If your goal is volume-rendering-like display (for example, you want to show "volume rendering" in Unity) then you need to use an actual volume renderer.

Andras


-----Original Message-----
From: vtkusers <[hidden email]> On Behalf Of Andaharoo
Sent: Sunday, December 2, 2018 11:27 PM
To: [hidden email]
Subject: Re: [vtkusers] Color PolyData based on VolumeData

So, if I understand right, you want to color the scalar/whatever field onto the poly? Just provide vtkPolyData with scalars and set the polyData's mapper scalar visibility to on.

If you give it 1 component scalars you can give the mapper a color function to map them. If you give it 3 component uchar it should just display them as color (as long as scalar visibility is on). Like so:


// Create the scalars array
vtkDoubleArray* scalars = vtkDoubleArray::New();
scalars->SetName("Scalars");
// For every point in the poly
for (vtkIdType i = 0; i < outputPoly->GetNumberOfPoints(); i++) {
        double pt[3];
        outputPoly->GetPoint(i, pt);

        scalars->InsertTuple1(i, trilinearSamplePoint(inputImage, pt[0], pt[1], pt[2], 1, 0)); } // Add the scalar data to the poly data
outputPoly->GetPointData()->AddArray(scalars);
outputPoly->GetPointData()->SetActiveScalars("Scalars");


This just goes through all the points in the polygon and gets the color from the image. Getting the color at a point in an image can be done by simply using the nearest color from the image (nearest neighbor) or better would be to use trilinear interpolation of the 8 nearest colors which I used in my example. Also note my code expects float image input:
static const int calcIndex(int x, int y, int z, int width, int height) { return x + width * (y + height * z); }

static float trilinearSamplePoint(vtkImageData* imageData, float x, float y, float z, int numComps, int comp) {
        float* imagePtr = static_cast<float*>(imageData->GetScalarPointer());
        double* spacing = imageData->GetSpacing();
        double* origin = imageData->GetOrigin();
        int* dim = imageData->GetDimensions();

        // We assume point x, y, z is in the image
        int xi = (x - origin[0]) / spacing[0];
        int yi = (y - origin[1]) / spacing[1];
        int zi = (z - origin[2]) / spacing[2];

        // Get the intensities at the 8 nearest neighbors
        float i000 = imagePtr[calcIndex(xi, yi, zi, dim[0], dim[1]) * numComps + comp];
        float i100 = imagePtr[calcIndex(xi + 1, yi, zi, dim[0], dim[1]) * numComps
+ comp];
        float i110 = imagePtr[calcIndex(xi + 1, yi + 1, zi, dim[0], dim[1] * numComps + comp)];
        float i010 = imagePtr[calcIndex(xi, yi + 1, zi, dim[0], dim[1]) * numComps
+ comp];

        float i001 = imagePtr[calcIndex(xi, yi, zi + 1, dim[0], dim[1]) * numComps
+ comp];
        float i101 = imagePtr[calcIndex(xi + 1, yi, zi + 1, dim[0], dim[1]) * numComps + comp];
        float i111 = imagePtr[calcIndex(xi + 1, yi + 1, zi + 1, dim[0], dim[1]) * numComps + comp];
        float i011 = imagePtr[calcIndex(xi, yi + 1, zi + 1, dim[0], dim[1]) * numComps + comp];

        // Get the fractional/unit distance from nearest neighbor 000
        float rx = xi * spacing[0] + origin[0]; // Position of node
        rx = (x - rx) / spacing[0]; // (Node - actual point position) / voxel width
        float ry = yi * spacing[1] + origin[1];
        ry = (y - ry) / spacing[1];
        float rz = zi * spacing[2] + origin[2];
        rz = (z - rz) / spacing[2];

        // Now we do the trilinear interpolation
        float ax = i000 + (i100 - i000) * rx;
        float bx = i010 + (i110 - i010) * rx;
        float cy = ax + (bx - ax) * ry;

        float dx = i001 + (i101 - i001) * rx;
        float ex = i011 + (i111 - i011) * rx;
        float fy = dx + (ex - dx) * ry;

        float gz = cy + (fy - cy) * rz;
        return gz;
}




--
Sent from: https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fvtk.1045678.n5.nabble.com%2FVTK-Users-f1224199.html&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091371821&amp;sdata=iziee8ltwEAWna3RvM3dma7A3p94NWMHCfIwclRpR0I%3D&amp;reserved=0
_______________________________________________
Powered by https://na01.safelinks.protection.outlook.com/?url=www.kitware.com&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091371821&amp;sdata=6Nv%2FWdIIIaOq%2B3peZlXxuGdAwkOwzORn7AfWthXw5Xo%3D&amp;reserved=0

Visit other Kitware open-source projects at https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.kitware.com%2Fopensource%2Fopensource.html&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=oEIpua2iKsLpUdJaY%2BmQdOhCHfusGYUZDRd3Ru67Qoo%3D&amp;reserved=0

Please keep messages on-topic and check the VTK FAQ at: https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.vtk.org%2FWiki%2FVTK_FAQ&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=7h18pSmCa3LapcNQJsqqxU7kyWVyaQzKbPYZPMCbDiw%3D&amp;reserved=0

Search the list archives at: https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fmarkmail.org%2Fsearch%2F%3Fq%3Dvtkusers&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=luU4EgKZZ9nD1OOxN5EGfBpA5HnK%2FyuVoN6C7PE%2Bzhc%3D&amp;reserved=0

Follow this link to subscribe/unsubscribe:
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpublic.kitware.com%2Fmailman%2Flistinfo%2Fvtkusers&amp;data=02%7C01%7Classo%40queensu.ca%7C34224811bdf846b9c4a808d658d7898f%7Cd61ecb3b38b142d582c4efb2838b925c%7C1%7C0%7C636794080091381830&amp;sdata=URr5vqvvkUuWmwaPIIGXXJzsmW2LFrYs7Z7jDcvt4ek%3D&amp;reserved=0
_______________________________________________
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

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
https://public.kitware.com/mailman/listinfo/vtkusers

_______________________________________________
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

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
https://public.kitware.com/mailman/listinfo/vtkusers