[vtkusers] 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.

[Andaharoo]

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: 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

[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 <vtkusers-***@public.kitware.com> On Behalf Of Andaharoo
Sent: Sunday, December 2, 2018 11:27 PM
To: ***@vtk.org
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

[Radosław Furmaniak]

Thanks a lot!
The solution from Andaharoo was what I was looking for.

Post by 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-----
Sent: Sunday, December 2, 2018 11:27 PM
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
// 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
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;
}
--
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
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
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
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
http://www.vtk.org/Wiki/VTK_FAQ
Search the list archives at: http://markmail.org/search/?q=vtkusers
https://public.kitware.com/mailman/listinfo/vtkusers