[vtkusers] Fwd: Click on surface and get image intensity value

[Renzo Ph]

Hello all,

I'm using vtkMarchingCubes to generate a surface of a niftii medical image.

    //Read image
    vtkSmartPointer<vtkNIFTIImageReader> reader = vtkSmartPointer<
vtkNIFTIImageReader>::New();
    reader->SetFileName (inputPath);
    reader->Update();

    //Trick to transform to vtkImageData
    vtkSmartPointer<vtkImageData> volume = vtkSmartPointer<vtkImageData>:
:New();
    volume->DeepCopy(reader->GetOutput());

    //Marching cubes
    vtkSmartPointer<vtkMarchingCubes> surface =  vtkSmartPointer<
vtkMarchingCubes>::New();
    surface->SetInputData(volume);
    surface->ComputeNormalsOn();
    surface->SetValue(0, 1);
    surface->Update();
    vtkPolyData* surfacePiece = surface->GetOutput();

Then, I try to get the original intensity values for every point of the
surface. This implies I have to map every surface point to its
corresponding voxel in the original medical image. As follows:

    for(int i = 0; i < surfacePiece->GetNumberOfPoints(); i++)
    {
       double p[3];
       p[0] = (p[0] - volume->GetOrigin()[0]) / volume->GetSpacing()[0];
       p[1] = (p[0] - volume->GetOrigin()[1]) / volume->GetSpacing()[1];
       p[2] = (p[0] - volume->GetOrigin()[2]) / volume->GetSpacing()[2];
       surfacePiece->GetPoint(i,p);

       double intensity = volume->GetScalarComponentAsDouble((int)(p[0]),
(int)(p[1]), (int)(p[2]), 0);
       cout << p[0] << " " << p[1] << " " << p[2] << " " << intensity <<
endl;
    }

However, I get an intensity value of 0 for most of the points of the
surface, which is not possible, as points with value 0 would not be part of
the surface. So, I think my mapping is not correct, and surface and
original image coordinates do not match.

I tried looking at close voxels, such as (int)(p[0]) + 1 or (int)(p[0]) -
1, but the result is the same. Also, I tried not to do the origin and
spacing conversion, but there were no improvements. I'm thinking about
orientation issues, but they do not seem to be the answer, and I cannot
find a way to check orientation by using VTK.

Has anybody had this problem before?

Sincerely,

Renzo Phellan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20160907/deee5ade/attachment.html>