[vtkusers] probe volume with thick line ( intensity profile with line radius or thickness )

[David Gobbi]

Hi Greg,

I don't know of any diagrams for vtkImageReslice but I think that I can
answer your questions.  First, though, I see a couple simple fixes for the
code:

The length of a vector is the 2-norm of the vector (sqrt of sum of squares):
  np.linalg.norm(vv, 1)  should be
  np.linalg.norm(vv, 2) or simply np.linalg.norm(vv)

The call to vtkMath::Perpendiculars() from Python is done as follows, by
creating the output value arrays and letting the method fill in the values:

  u = np.zeros(3)
  w = np.zeros(3)
  vtk.vtkMath.Perpendiculars(nn, u, w, 0)

With VTK's wrapping, if the C++ method returns values via the parameters,
then the Python method always does the same.


With regards to SetOutputExtent(), the "Extent" in VTK is a set of six
values [Xmin, Xmax, Ymin, Ymax, Zmin, Zmax] that give the limits of a block
of voxels.  A line of N+1 voxels in the X direction that starts at zero is
(0, N, 0, 0, 0, 0).  The X, Y, Z for the extent are voxel indices, i.e.
they are integers.

So the Extent in VTK serves a similar purpose to the shape in numpy, with
three fundamental differences:

1) The Extent is always 3D, whereas numpy arrays are ND
2) Extent has six values, since it gives the lower and upper value in each
dimension rather than just giving the size in each dimension
3) The numpy shape is ordered from the slowest-varying to the
fastest-varying index, whereas Extent is the opposite

An Extent of (0, Xsize-1, 0, Ysize-1, 0, Zsize-1) is equivalent to numpy
shape of (Zsize, Ysize, Xsize) if the numpy array is C-contiguous.

Hopefully that didn't muddy up the issue too much.


About the Origin, its purpose here is as follows.  Let's say that you want
to use three lines instead of one, to give the line some thickness in one
direction:

  reslice.SetOutputExtent(0, N, 0, 2, 0, 0)  # gives a (N+1)*3*1 block of
voxels

This block needs to be shifted so the middle line is the one at the desired
position, and this shift is done by setting the Origin:

  reslice.SetOutputOrigin(0.0, -1.0*spacing[1], 0.0)

Setting the OutputOrigin like this causes the block to be shifted down by
one sample spacing in the y direction.  More lines would require a larger
shift.  Remember that the idea is that you the ResliceAxes transformation
is going to put one end of this block at the desired starting point for
your line.

And note that this block has to define a line (or lines) that go along the
X-axis!  That's the whole point, since the transformation matrix will
rotate this X-axis to the desired line orientation for your input image.
So get rid of the "if" statements that set up the extent in different
directions, since the line orientation is already dealt with by the
rotation in the transformation.

 - David



On Sat, Feb 3, 2018 at 10:16 PM, gregthom992 <gregthom992 at gmail.com> wrote:

> Hi David,
> So I am almost there but there is some magic you have to explain to me. Is
> there some graphics that explans how imagereslice works ?
>
> For instance:
> How do you know to do
> reslice.SetOutputExtent(0, N, 0, 0, 0, 0)  # this gives one line of output
> and not
> reslice.SetOutputExtent(0, 0, 0, 0, 0, N)  # this gives one line of output
>
> Both give one line of output but I am lost here.
>
> Also:
>
> Please I know I can get multiple lines for instance 2 lines by:
> reslice.SetOutputExtent(0, N, 0, 1, 0, 0)  # this gives two lines of output
>
> What were you saying about fiddling with origin ? Below is my complete code
> which doesn't give me the expected results
>
>
>
> ===================
>
>             vtkimage = SITKtoVTK(itkimage)
>             vtkimage.SetOrigin(itkimage.GetOrigin())
>             print itkimage
>             print vtkimage
>
>             p0 = np.array(prf_start)
>             p1 = np.array(prf_stop)
>             vv = p1 - p0
>             nn = vv/np.linalg.norm(vv, 1)
>
>             u, w = vtkMathPerpendiculars(nn, 3*[0], 3*[0], 0)
>
>             #print vtkimage, itkimage
>             print vv
>             print nn
>             print u, ' ',  w
>             #
>             oM2 = vtk.vtkMatrix4x4()
>             oM2.DeepCopy((
>                 nn[0], u[0], w[0], p0[0],
>                 nn[1], u[1], w[1], p0[1],
>                 nn[2], u[2], w[2], p0[2],
>                 0, 0, 0, 1 ))
>             print oM2
>             #
>             reslice = vtk.vtkImageReslice()
>             reslice.SetInput(vtkimage)
>             #reslice.SetOutputOrigin(itkimage.GetOrigin()[0],
> itkimage.GetOrigin()[1], itkimage.GetOrigin()[2])
>             reslice.SetOutputOrigin(0, 0, 0)
>             reslice.SetOutputSpacing(itkimage.GetSpacing()[0],
> itkimage.GetSpacing()[1], itkimage.GetSpacing()[2])
>             if prf_params[0] in ['depth']:
>                 reslice.SetOutputExtent(0, itkimage.GetDepth()-1, 0, 0, 0,
> 0)  # this gives one line of output
>             if prf_params[0] in ['inline']:
>                 reslice.SetOutputExtent(0, 0, 0, itkimage.GetHeight()-1, 0,
> 0)  # this gives one line of output
>             if prf_params[0] in ['crossline']:
>                 reslice.SetOutputExtent(0, 0, 0, 0, 0,
> itkimage.GetWidth()-1)  # this gives one line of output
>             reslice.SetInterpolationModeToLinear()
>             reslice.SetResliceAxes(oM2)
>             reslice.Update()
>             vtkimage2 = reslice.GetOutput()
>
>             print vtkimage2
>
>             numpy_data =
> numpy_support.vtk_to_numpy(vtkimage2.GetPointData().GetScalars())
>
>             print numpy_data
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://vtk.org/pipermail/vtkusers/attachments/20180204/d45f0323/attachment.html>