[vtkusers] vtkProbeFilter not probing volume

[Seth Gilchrist]

I played with this for quite a while, but never could get it to work. In
the end I made an application specific probe. The code, for others who have
the same issue is:

#include <vtkUnstructuredGridReader.h>
#include <vtkPolyData.h>
#include <vtkUnstructuredGrid.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkSmartPointer.h>
#include <vtkPointData.h>
#include <vtkDoubleArray.h>
#include <vtkCellLocator.h>
#include <vtkWedge.h>
#include <vtkCell.h>


int main(int argc, char **argv)
{
    // Read the volume
    vtkSmartPointer<vtkUnstructuredGridReader> volumeReader =
        vtkSmartPointer<vtkUnstructuredGridReader>::New();
    volumeReader->SetFileName("/home/seth/Desktop/volume.vtk");
    volumeReader->Update();

    // Read the surface
    vtkSmartPointer<vtkXMLPolyDataReader> polyReader =
        vtkSmartPointer<vtkXMLPolyDataReader>::New();
    polyReader->SetFileName("/home/seth/Desktop/aligned.vtp");
    polyReader->Update();


    // remove pre-existing array data
    vtkSmartPointer<vtkPolyData> inputSurface = polyReader->GetOutput();
    unsigned int numberOfArrays =
inputSurface->GetPointData()->GetNumberOfArrays();
    for (unsigned int i = 0; i < numberOfArrays; ++i)
    {
        inputSurface->GetPointData()->RemoveArray(i);
    }

    std::cout<<"Images Read"<<std::endl;

    // create a new array to hold the data
    vtkSmartPointer<vtkDoubleArray> newArray =
vtkSmartPointer<vtkDoubleArray>::New();
    newArray->SetNumberOfTuples(inputSurface->GetNumberOfPoints());
    newArray->SetNumberOfComponents(1);
    newArray->SetName("Extracted MinPStrain");
    inputSurface->GetPointData()->AddArray(newArray);

    // create a cell locator to aid in finding the cells that points belong
to
    vtkSmartPointer<vtkCellLocator> cellLocator =
        vtkSmartPointer<vtkCellLocator>::New();
    cellLocator->SetDataSet(volumeReader->GetOutput());
    cellLocator->BuildLocator();

    // these will be used in the loop to hold data
    double blankTuple = 0.0;
    vtkSmartPointer<vtkWedge> cCell;
    double cPoint[3];
    double closestPoint[3];
    int subId;
    double pcoords[3];
    double dist2;
    double weights[6];
    double cellData[6];
    double cValue;
    vtkSmartPointer<vtkIdList> cellPoints;

    // loop through each point in the input surface
    for (unsigned int i = 0; i < inputSurface->GetNumberOfPoints(); ++i)
    {
        // get the point location
        inputSurface->GetPoint(i,cPoint);
        // find the cell that contains the point
        vtkIdType cCellNo = cellLocator->FindCell(cPoint);
        if (cCellNo == -1)  // if the point is outside of cells, enter a
value of zero in the data array
        {

inputSurface->GetPointData()->GetArray(0)->SetTuple(i,&blankTuple);
            continue;
        }
        // get the point IDs that define the containing cell
        cellPoints =
volumeReader->GetOutput()->GetCell(cCellNo)->GetPointIds();
        // get the strain values for each point in the containing cell
        cellData[0] =
*volumeReader->GetOutput()->GetPointData()->GetArray(0)->GetTuple(cellPoints->GetId(0));
        cellData[1] =
*volumeReader->GetOutput()->GetPointData()->GetArray(0)->GetTuple(cellPoints->GetId(1));
        cellData[2] =
*volumeReader->GetOutput()->GetPointData()->GetArray(0)->GetTuple(cellPoints->GetId(2));
        cellData[3] =
*volumeReader->GetOutput()->GetPointData()->GetArray(0)->GetTuple(cellPoints->GetId(3));
        cellData[4] =
*volumeReader->GetOutput()->GetPointData()->GetArray(0)->GetTuple(cellPoints->GetId(4));
        cellData[5] =
*volumeReader->GetOutput()->GetPointData()->GetArray(0)->GetTuple(cellPoints->GetId(5));
        // use the EvaluatePosition method of the cell to get the
interpolation function weight values. The rest of the data is not used

volumeReader->GetOutput()->GetCell(cCellNo)->EvaluatePosition(cPoint,closestPoint,subId,pcoords,dist2,weights);
        // calculate the value at the point by multiplying each weight by
the data
        cValue = cellData[0]*weights[0] + cellData[1]*weights[1] +
cellData[2]*weights[2] + cellData[3]*weights[3] + cellData[4]*weights[4] +
cellData[5]*weights[5];
        // set the data.
        inputSurface->GetPointData()->GetArray(0)->SetTuple(i,&cValue);
    }

    // write the probed surface out
    vtkSmartPointer<vtkXMLPolyDataWriter> polyWirter =
        vtkSmartPointer<vtkXMLPolyDataWriter>::New();
    polyWirter->SetInput(inputSurface);
    polyWirter->SetFileName("/home/seth/Desktop/testProbe.vtp");
    polyWirter->Write();
}

*Seth Gilchrist, MASc, PhD Candidate
+---------------------------------------+
Orthopaedic and Injury Biomechanics Group
UBC Department of Mechanical Engineering
6th Floor-2635 Laurel Street
Vancouver, BC V5Z-1M9
+---------------------------------------+
seth at mech.ubc.ca
Fax 604-675-2576*


On Thu, May 23, 2013 at 12:08 PM, Seth Gilchrist <seth at mech.ubc.ca> wrote:

> Hello all,
> I'm having a problem with the vtkProbeFilter. I have a vtkUnstructuredGrid
> volume with one data field called MinPStrain and I am trying to probe it
> with a vtkPolyData surface. The output indicates that the two interact in
> only a very small portion of the surface, when in reality they almost
> completely overlap.
>
> There are some screen shots at:
> https://plus.google.com/photos/101839812291734892344/albums/5881253191058882961?authkey=CIOdsYjQi8zKjQE
>
> On the left is the surface and volume in their actual positions. On the
> right is the output of the probe, translated so that it can be viewed.
>
> Most recently I have tried using the vtkAppendFilter to convert the
> vtkPolyData to a vtkUnstructuredGrid, but to no avail. The latest code is:
>
> #include <vtkUnstructuredGridReader.h>
> #include <vtkXMLUnstructuredGridWriter.h>
> #include <vtkPolyData.h>
> #include <vtkUnstructuredGrid.h>
> #include <vtkXMLPolyDataReader.h>
> #include <vtkXMLPolyDataWriter.h>
> #include <vtkProbeFilter.h>
> #include <vtkSmartPointer.h>
> #include <vtkAppendFilter.h>
> #include <vtkPointData.h>
>
>
> int main(int argc, char **argv)
> {
>     // Read the volume
>     vtkSmartPointer<vtkUnstructuredGridReader> volumeReader =
>         vtkSmartPointer<vtkUnstructuredGridReader>::New();
>     volumeReader->SetFileName("/home/seth/Desktop/volume.vtk");
>     volumeReader->Update();
>
>     // Read the surface
>     vtkSmartPointer<vtkXMLPolyDataReader> polyReader =
>         vtkSmartPointer<vtkXMLPolyDataReader>::New();
>     polyReader->SetFileName("/home/seth/Desktop/aligned.vtp");
>     polyReader->Update();
>
>     // remove pre-existing array data
>     vtkSmartPointer<vtkPolyData> inputSurface = polyReader->GetOutput();
>     unsigned int numberOfArrays =
> inputSurface->GetPointData()->GetNumberOfArrays();
>     for (unsigned int i = 0; i < numberOfArrays; ++i)
>     {
>         inputSurface->GetPointData()->RemoveArray(i);
>     }
>
>     std::cout<<"Images Read"<<std::endl;
>
>     // convert the surface from vtkPolyData to vtkUnstructuredGrid
>     vtkSmartPointer<vtkAppendFilter> appender =
>         vtkSmartPointer<vtkAppendFilter>::New();
>     appender->SetInput(inputSurface);
>     appender->Update();
>
>     // probe the volume with the surface
>     vtkSmartPointer<vtkProbeFilter> probe =
>         vtkSmartPointer<vtkProbeFilter>::New();
>     probe->SetInput(appender->GetOutput());
>     probe->SetSource(volumeReader->GetOutput());
>
>     // write the probed surface out
>     vtkSmartPointer<vtkXMLUnstructuredGridWriter> polyWirter =
>         vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
>     polyWirter->SetInputConnection(probe->GetOutputPort());
>     polyWirter->SetFileName("/home/seth/Desktop/testProbe.vtu");
>     polyWirter->Write();
> }
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20130523/80bd21a1/attachment.htm>