Intersection of the vtkImageData with arbitrary box

Support
1

I have some 3D vtkImageData and want to find a “subvolume” (I’m not sure I use this term in the right way here) that is intersection of this vtkImageData and arbitrary box. This box can be translated and rotated. I wonder what data structures and filters should I use.

  1. I read it’s better to use vtkStructuredGrid instead of vtkImageData, right?
  2. Could I create some variable of type vtkStructuredGrid from vtkImageData and then translate and rotate it?
  3. The same for box. Should I create something like vtkCube and translate and rotate it?
  4. But then I want to find their intersection. I suspect it can be vtkStructuredGrid as well? Can I use some interpolation on the cells on the bounds?
  5. And can I then transform this vtkStructuredGrid into some minimal vtkImageData with zeros values in the extra cells?

Maybe I’m completely wrong and this problem has another/simpler solution. Could someone direct me, please?

2

If everything is axis aligned I believe it’s straightforward (famous last words :-)). Take the lower-left point of the box and determine the voxel cell/indices the point is contained in (via vtkImageData::ComputeStructuredCoordinates()), then take the upper-right corner and do the same thing. The resulting two sets of indices define a subvolume in the volume etc. If the volume and box are not aligned then you’ll have to use transformations to position them into the same space.

3

it’s impossible to tranform vtkImageData, I wrote about it

4

The box can be transformed. vtkImageData has direction matrices.

5

Unfortunately, your answer can’t help me. I don’t know which filter should I use to find intersection. Obviously the result can’t be vtkImageData. It should be vtkStructuredGrid?

6

You can use the vtkImageReslice filter to get an arbitrarily oriented subvolume of a vtkImageData. The output is a vtkImageData.

7

Thank you very much. But I can’t understand how works this filter even in the simplest case.

I have vtkImageData like this

00 01 02 03 04
05 06 07 08 09
10 11 12 13 14

But origin of this data is (2, 3). So, (1):

    vtkSmartPointer<vtkImageData> localImage = vtkSmartPointer<vtkImageData>::New();

    int minExtentX = 0;
    int maxExtentX = 4;
    int minExtentY = 0;
    int maxExtentY = 2;
    localImage->SetExtent(minExtentX, maxExtentX, minExtentY, maxExtentY, 0, 0);

    double originX = 2.0;
    double originY = 3.0;
    localImage->SetOrigin(originX, originY, 0.0);

    vtkFloatArray* floats = vtkFloatArray::New();
    int dataNumbers = (maxExtentY - minExtentY + 1) * (maxExtentX - minExtentX + 1);
    for (int i = 0; i < dataNumbers; ++i)
    {
        floats->InsertNextTuple1(i);
    }

    localImage->GetPointData()->SetScalars(floats);

Next I create a default vtkTransformation and apply it to source image (2):

    vtkSmartPointer<vtkTransform> localToGlobal = vtkSmartPointer<vtkTransform> ::New();

    vtkSmartPointer<vtkMatrix4x4> m = vtkSmartPointer<vtkMatrix4x4>::New();
    localToGlobal->GetMatrix(m);

    vtkSmartPointer<vtkImageReslice> reslicedImage = vtkSmartPointer<vtkImageReslice>::New();
    reslicedImage->SetInputData(localImage);
    reslicedImage->SetResliceTransform(localToGlobal);
    reslicedImage->GenerateStencilOutputOn();
    reslicedImage->Update();

Then I check the results (4):

double* reslicedOutputOrigion = reslicedImage->GetOutputOrigin();
int* reslicedOutputExtents = reslicedImage->GetOutputExtent();
double* reslicedOutputSpacing = reslicedImage->GetOutputSpacing();
vtkImageStencilData* stencil = reslicedImage->GetStencilOutput();

int stencilExtentMinX, stencilExtentMaxX, stencilExtentMinY, stencilExtentMaxY, stencilExtentMinZ, stencilExtentMaxZ;
stencil->GetExtent(stencilExtentMinX, stencilExtentMaxX, stencilExtentMinY, stencilExtentMaxY, stencilExtentMinZ, stencilExtentMaxZ);

for (int x = stencilExtentMinX; x <= stencilExtentMaxX; ++x){
    for (int y = stencilExtentMinY; y <= stencilExtentMaxY; ++y){
        for (int z = stencilExtentMinZ; z <= stencilExtentMaxZ; ++z){
            cout << "(" << x << ", " << y << ", " << z << ")" << (stencil->IsInside(x, y, z) ? " is inside" : " is outside") << endl;
            if (stencil->IsInside(x, y, z)){
                cout << *static_cast<float*>(reslicedImage->GetOutput()->GetScalarPointer(x, y, z)) << endl;
                }
            }
        }
    }

Stencil shows that data are the same but the origin of resliced image is (0, 0, 0).
But it is obviously that output origin has to be the same as before (i.e. (2, 3)).
All the extends are 0 as well. I do not understand why it is.
I tried to use reslicedImage->SetOutputExtends and reslicedImage->SetOutputOrigin but it confused me even more.

What I have missed?
How can I cut the part of transformed data and use them for new vtkImageData?

8

vtkImageReslice is probably the most complicated filter in VTK, because there are several distinct ways to use it and it has a complex internal state machine (so that it’s behavior depends on what order you call its methods). For example, you can enable automatic computation of output geometry in the filter, but if you then set anything related to the output geometry then the filter will ignore the automatic computation request.

I would recommend to start from a test or example that uses this filter; or step through all its setter methods and Update() method call in a debugger to see where things work differently compared to what you expect.

9

Yes, thank you very much.
But I spent few days in the debugger and don’t understand how you are able to force it to do what you need. I saw an example but don’t understand how it works still. Is there either “how to” or some “step by step” explanation? I saw few c++ and python examples and read vtk help but I am not able to write my own code properly

10

I agree, it is hard. I’ve been working with VTK full time for more than 10 years but still find usage of this filter very challenging.

However, all the tests and examples in VTK repository and the VTK Examples website work well, so as long as you choose one and follow it exactly (what methods of the filter are called, in what order) then your code is guaranteed to work well, too.

Spending several days with trying to understand how to use this filter may seem a lot, but this is a very versatile and quite nicely optimized filter that can take care of all image resampling needs (including various interpolations, resampling with non-linear transforms, etc), so it should worth the investment.