Commit 0a9fb2bd authored by Alexis Girault's avatar Alexis Girault Committed by Kitware Robot
Browse files

Merge topic 'orientation-image-isocontour'

1aa07e3e vtkFlyingEdges2D: Add support for image orientation
b5dc4423 vtkFlyingEdges3D: Add support for image orientation
691e1528 vtkMarchingSquares: Add support for image orientation
9414d605 vtkMarchingCubes: Update tests to account for orientation
4e8b769f vtkMarchingCubes: Add support for image orientation
d17c2aab vtkMarchingCubes: Add python test
c392f54f New helper class to transform out of image data
9e179c2a

 Add vtkMatrix4x4::IsIdentity() method
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Rejected-by: Alexis Girault's avatarAlexis Girault <alexis.girault@kitware.com>
Acked-by: Ken Martin's avatarKen Martin <ken.martin@kitware.com>
Merge-request: !5590
parents 7b3b4837 1aa07e3e
Pipeline #138491 failed with stage
in 0 seconds
......@@ -91,6 +91,7 @@ set(classes
vtkHyperTreeGridOrientedGeometryCursor
vtkImageData
vtkImageIterator
vtkImageTransform
vtkImplicitBoolean
vtkImplicitDataSet
vtkImplicitFunction
......
/*=========================================================================
Program: Visualization Toolkit
Module: vtkImageTransform.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkImageTransform.h"
#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkImageData.h"
#include "vtkMatrix3x3.h"
#include "vtkMatrix4x4.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPointSet.h"
#include "vtkSMPTools.h"
vtkStandardNewMacro(vtkImageTransform);
//============================================================================
//----------------------------------------------------------------------------
// Functors to support threaded execution
namespace { //anonymous
template <typename T>
struct InPlaceTranslatePoints
{
T *Points;
double *Translation;
InPlaceTranslatePoints(double t[3], T *pts) : Points(pts), Translation(t) {}
void operator()(vtkIdType ptId, vtkIdType endPtId)
{
T *pIn = this->Points + 3*ptId;
T *pOut = pIn;
for ( ; ptId < endPtId; ++ptId )
{
*pIn++ = *pOut++ + this->Translation[0];
*pIn++ = *pOut++ + this->Translation[1];
*pIn++ = *pOut++ + this->Translation[2];
}
}
// Interface to vtkSMPTools
static void Execute(double t[3], vtkIdType num, T *pts)
{
InPlaceTranslatePoints<T> translate(t,pts);
vtkSMPTools::For(0,num,translate);
}
};//InPlaceTransformPoints
template <typename T>
struct InPlaceTransformPoints
{
T *Points;
vtkMatrix4x4 *M4;
InPlaceTransformPoints(vtkMatrix4x4 *m4, T *pts) : Points(pts), M4(m4) {}
void operator()(vtkIdType ptId, vtkIdType endPtId)
{
T *pIn = this->Points + 3*ptId;
T tmp[3] = {0, 0, 0};
for ( ; ptId < endPtId; ++ptId )
{
tmp[0] = M4->GetElement(0, 0) * pIn[0] +
M4->GetElement(0, 1) * pIn[1] +
M4->GetElement(0, 2) * pIn[2] +
M4->GetElement(0, 3);
tmp[1] = M4->GetElement(1, 0) * pIn[0] +
M4->GetElement(1, 1) * pIn[1] +
M4->GetElement(1, 2) * pIn[2] +
M4->GetElement(1, 3);
tmp[2] = M4->GetElement(2, 0) * pIn[0] +
M4->GetElement(2, 1) * pIn[1] +
M4->GetElement(2, 2) * pIn[2] +
M4->GetElement(2, 3);
*pIn++ = tmp[0];
*pIn++ = tmp[1];
*pIn++ = tmp[2];
}
}
// Interface to vtkSMPTools
static void Execute(vtkMatrix4x4 *m4, vtkIdType num, T *pts)
{
InPlaceTransformPoints<T> transform(m4,pts);
vtkSMPTools::For(0,num,transform);
}
};//InPlaceTransformPoints
template <typename T>
struct InPlaceTransformNormals
{
T *Normals;
vtkMatrix3x3 *M3;
double Determinant;
double *Spacing;
InPlaceTransformNormals(vtkMatrix3x3 *m3, double* spacing, T *n) :
Normals(n),
M3(m3),
Determinant(m3->Determinant()),
Spacing(spacing)
{}
void operator()(vtkIdType ptId, vtkIdType endPtId)
{
T *nIn = this->Normals + 3*ptId;
T tmp[3] = {0, 0, 0};
T toUnit = 0;
for ( ; ptId < endPtId; ++ptId )
{
nIn[0] = nIn[0] / this->Spacing[0];
nIn[1] = nIn[1] / this->Spacing[1];
nIn[2] = nIn[2] / this->Spacing[2];
tmp[0] = M3->GetElement(0, 0) * nIn[0] +
M3->GetElement(0, 1) * nIn[1] +
M3->GetElement(0, 2) * nIn[2];
tmp[1] = M3->GetElement(1, 0) * nIn[0] +
M3->GetElement(1, 1) * nIn[1] +
M3->GetElement(1, 2) * nIn[2];
tmp[2] = M3->GetElement(2, 0) * nIn[0] +
M3->GetElement(2, 1) * nIn[1] +
M3->GetElement(2, 2) * nIn[2];
tmp[0] *= this->Determinant;
tmp[1] *= this->Determinant;
tmp[2] *= this->Determinant;
toUnit = 1 / sqrt(tmp[0]*tmp[0] + tmp[1]*tmp[1] + tmp[2]*tmp[2]);
*nIn++ = tmp[0] * toUnit;
*nIn++ = tmp[1] * toUnit;
*nIn++ = tmp[2] * toUnit;
}
}
// Interface to vtkSMPTools
static void Execute(vtkMatrix3x3 *m3, double* spacing, vtkIdType num, T *n)
{
InPlaceTransformNormals<T> transform(m3,spacing,n);
vtkSMPTools::For(0,num,transform);
}
};//InPlaceTransformNormals
template <typename T>
struct InPlaceTransformVectors
{
T *Vectors;
vtkMatrix3x3 *M3;
double *Spacing;
InPlaceTransformVectors(vtkMatrix3x3 *m3, double* spacing, T *v) :
Vectors(v),
M3(m3),
Spacing(spacing)
{}
void operator()(vtkIdType ptId, vtkIdType endPtId)
{
T *nIn = this->Vectors + 3*ptId;
T tmp[3] = {0, 0, 0};
for ( ; ptId < endPtId; ++ptId )
{
nIn[0] = nIn[0] / this->Spacing[0];
nIn[1] = nIn[1] / this->Spacing[1];
nIn[2] = nIn[2] / this->Spacing[2];
tmp[0] = M3->GetElement(0, 0) * nIn[0] +
M3->GetElement(0, 1) * nIn[1] +
M3->GetElement(0, 2) * nIn[2];
tmp[1] = M3->GetElement(1, 0) * nIn[0] +
M3->GetElement(1, 1) * nIn[1] +
M3->GetElement(1, 2) * nIn[2];
tmp[2] = M3->GetElement(2, 0) * nIn[0] +
M3->GetElement(2, 1) * nIn[1] +
M3->GetElement(2, 2) * nIn[2];
*nIn++ = tmp[0];
*nIn++ = tmp[1];
*nIn++ = tmp[2];
}
}
// Interface to vtkSMPTools
static void Execute(vtkMatrix3x3 *m3, double* spacing, vtkIdType num, T *v)
{
InPlaceTransformVectors<T> transform(m3,spacing,v);
vtkSMPTools::For(0,num,transform);
}
};//InPlaceTransformVectors
}//anonymous namespace
//============================================================================
//----------------------------------------------------------------------------
// Here is the VTK class proper.
vtkImageTransform::vtkImageTransform()
{
}
//----------------------------------------------------------------------------
void vtkImageTransform::TransformPointSet(vtkImageData *im, vtkPointSet *ps)
{
// Check input
if ( im == nullptr || ps == nullptr )
{
return;
}
// Nothing to do if the direction matrix is the identity
vtkMatrix4x4 *m4 = im->GetIndexToPhysicalMatrix();
if ( m4->IsIdentity() )
{
return;
}
// Make sure points are available
vtkIdType numPts = ps->GetNumberOfPoints();
if ( numPts < 1 )
{
return;
}
// Grab the points-related-data and process as appropriate
vtkDataArray *pts = ps->GetPoints()->GetData();
vtkMatrix3x3 *m3 = im->GetDirectionMatrix();
double *ar = im->GetSpacing();
// If there is no rotation or spacing, only translate
if ( m3->IsIdentity() && ar[0] == 1 && ar[1] == 1 && ar[2] == 1 )
{
vtkImageTransform::TranslatePoints(im->GetOrigin(),pts);
return;
}
vtkImageTransform::TransformPoints(m4,pts);
vtkDataArray *normals = ps->GetPointData()->GetNormals();
if ( normals != nullptr )
{
vtkImageTransform::TransformNormals(m3,ar,normals);
}
vtkDataArray *vectors = ps->GetPointData()->GetVectors();
if ( vectors != nullptr )
{
vtkImageTransform::TransformVectors(m3,ar,vectors);
}
// Grab the cells-related-data and process as appropriate
normals = ps->GetCellData()->GetNormals();
if ( normals != nullptr )
{
vtkImageTransform::TransformNormals(m3,ar,normals);
}
vectors = ps->GetCellData()->GetVectors();
if ( vectors != nullptr )
{
vtkImageTransform::TransformVectors(m3,ar,vectors);
}
}
//----------------------------------------------------------------------------
void vtkImageTransform::TranslatePoints(double *t, vtkDataArray *da)
{
void *pts = da->GetVoidPointer(0);
vtkIdType num = da->GetNumberOfTuples();
switch ( da->GetDataType() )
{
vtkTemplateMacro(InPlaceTranslatePoints<VTK_TT>::
Execute(t,num,(VTK_TT*)pts));
}
}
//----------------------------------------------------------------------------
void vtkImageTransform::TransformPoints(vtkMatrix4x4 *m4,vtkDataArray *da)
{
void *pts = da->GetVoidPointer(0);
vtkIdType num = da->GetNumberOfTuples();
switch ( da->GetDataType() )
{
vtkTemplateMacro(InPlaceTransformPoints<VTK_TT>::
Execute(m4,num,(VTK_TT*)pts));
}
}
//----------------------------------------------------------------------------
void vtkImageTransform::TransformNormals(vtkMatrix3x3 *m3,
double spacing[3],
vtkDataArray *da)
{
void *n = da->GetVoidPointer(0);
vtkIdType num = da->GetNumberOfTuples();
switch ( da->GetDataType() )
{
vtkTemplateMacro(InPlaceTransformNormals<VTK_TT>::
Execute(m3,spacing,num,(VTK_TT*)n));
}
}
//----------------------------------------------------------------------------
void vtkImageTransform::TransformVectors(vtkMatrix3x3 *m3,
double spacing[3],
vtkDataArray *da)
{
void *v = da->GetVoidPointer(0);
vtkIdType num = da->GetNumberOfTuples();
switch ( da->GetDataType() )
{
vtkTemplateMacro(InPlaceTransformVectors<VTK_TT>::
Execute(m3,spacing,num,(VTK_TT*)v));
}
}
//----------------------------------------------------------------------------
void vtkImageTransform::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkImageTransform.h
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
/**
* @class vtkImageTransform
* @brief helper class to transform output of non-axis-aligned images
*
* vtkImageTransform is a helper class to transform the output of
* image filters (i.e., filter that input vtkImageData) by applying the
* Index to Physical transformation frmo the input image, which can
* include origin, spacing, direction. The transformation process is
* threaded with vtkSMPTools for performance.
*
* Typically in application the single method TransformPointSet() is
* invoked to transform the output of an image algorithm (assuming
* that the image's direction/orientation matrix is
* non-identity). Note that vtkPointSets encompass vtkPolyData as well
* as vtkUnstructuredGrids. In the future other output types may be
* added. Note that specific methods for transforming points, normals,
* and vectors is also provided by this class in case additional
* output data arrays need to be transformed (since
* TransformPointSet() only processes data arrays labeled as points,
* normals, and vectors).
*
* @warning
* This class has been threaded with vtkSMPTools. Using TBB or other
* non-sequential type (set in the CMake variable
* VTK_SMP_IMPLEMENTATION_TYPE) may improve performance significantly.
*
*/
#ifndef vtkImageTransform_h
#define vtkImageTransform_h
#include "vtkCommonDataModelModule.h" // For export macro
#include "vtkObject.h"
class vtkDataArray;
class vtkImageData;
class vtkMatrix3x3;
class vtkMatrix4x4;
class vtkPointSet;
class VTKCOMMONDATAMODEL_EXPORT vtkImageTransform : public vtkObject
{
public:
//@{
/**
* Standard methods for construction, type information, printing.
*/
static vtkImageTransform *New();
vtkTypeMacro(vtkImageTransform,vtkObject);
void PrintSelf(ostream& os, vtkIndent indent) override;
//@}
/**
* Given a vtkImageData (and hence its associated orientation
* matrix), and an instance of vtkPointSet, transform its points, as
* well as any normals and vectors, associated with the
* vtkPointSet. This is a convenience function, internally it calls
* TranslatePoints(), TransformPoints(), TransformNormals(), and/or
* TransformVectors() as appropriate. Note that both the normals and
* vectors associated with the point and cell data are transformed.
*/
static void TransformPointSet(vtkImageData *im, vtkPointSet *ps);
/**
* Given x-y-z points represented by a vtkDataArray,
* translate the points using the image origin. This
* method is useful if there is no orientation or
* spacing to apply.
*/
static void TranslatePoints(double *t, vtkDataArray *da);
/**
* Given x-y-z points represented by a vtkDataArray,
* transform the points using the matrix provided.
*/
static void TransformPoints(vtkMatrix4x4 *m4, vtkDataArray *da);
/**
* Given three-component normals represented by a vtkDataArray,
* transform the normals using the matrix provided.
*/
static void TransformNormals(vtkMatrix3x3 *m3,
double spacing[3],
vtkDataArray *da);
/**
* Given three-component vectors represented by a vtkDataArray,
* transform the vectors using the matrix provided.
*/
static void TransformVectors(vtkMatrix3x3 *m3,
double spacing[3],
vtkDataArray *da);
protected:
vtkImageTransform();
~vtkImageTransform() {}
private:
vtkImageTransform(const vtkImageTransform&) = delete;
void operator=(const vtkImageTransform&) = delete;
};
#endif
......@@ -87,6 +87,11 @@ public:
{ vtkMatrix4x4::Identity(*this->Element); this->Modified();}
static void Identity(double elements[16]);
/**
* Returns true if this matrix is equal to the identity matrix.
*/
bool IsIdentity();
/**
* Matrix Inversion (adapted from Richard Carling in "Graphics Gems,"
* Academic Press, 1990).
......@@ -228,5 +233,16 @@ inline void vtkMatrix4x4::SetElement(int i, int j, double value)
}
}
//----------------------------------------------------------------------------
inline bool vtkMatrix4x4::IsIdentity()
{
double *M = *this->Element;
return
M[0] == 1.0 && M[1] == 0.0 && M[2] == 0.0 && M[3] == 0.0 &&
M[4] == 0.0 && M[5] == 1.0 && M[6] == 0.0 && M[7] == 0.0 &&
M[8] == 0.0 && M[9] == 0.0 && M[10] == 1.0 && M[11] == 0.0 &&
M[12] == 0.0 && M[13] == 0.0 && M[14] == 0.0 && M[15] == 1.0;
}
#endif
......@@ -16,6 +16,7 @@
#include "vtkImageData.h"
#include "vtkCellArray.h"
#include "vtkImageTransform.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
......@@ -81,9 +82,7 @@ public:
// Internal variables used by the various algorithm methods. Interfaces VTK
// image data in a form more convenient to the algorithm.
vtkIdType Dims[2];
double Origin[3];
double Spacing[3];
double Z;
int K;
int Axis0;
int Min0;
int Max0;
......@@ -103,14 +102,6 @@ public:
// Instantiate and initialize key data members.
vtkFlyingEdges2DAlgorithm();
// Adjust the origin to the lower-left corner of the volume (if necessary)
void AdjustOrigin(int updateExt[6])
{
this->Origin[0] = this->Origin[0] + this->Spacing[0]*updateExt[0];
this->Origin[1] = this->Origin[1] + this->Spacing[1]*updateExt[2];
this->Origin[2] = this->Origin[2] + this->Spacing[2]*updateExt[4];
}
// The three passes of the algorithm.
void ProcessXEdge(double value, T* inPtr, vtkIdType row); //PASS 1
void ProcessYEdges(vtkIdType row); //PASS 2
......@@ -165,24 +156,24 @@ public:
}
// Interpolate along a pixel axes edge.
void InterpolateAxesEdge(double value, T *s0, float x0[3], T* s1,
float x1[3], vtkIdType vId)
void InterpolateAxesEdge(double value, T *s0, int ijk0[3], T* s1,
int ijk1[3], vtkIdType vId)
{
double t = (value - *s0) / (*s1 - *s0);
float *x = this->NewPoints + 3*vId;
x[0] = x0[0] + t*(x1[0]-x0[0]);
x[1] = x0[1] + t*(x1[1]-x0[1]);
x[2] = this->Z;
x[0] = ijk0[0] + t*(ijk1[0]-ijk0[0]) + this->Min0;
x[1] = ijk0[1] + t*(ijk1[1]-ijk0[1]) + this->Min1;
x[2] = this->K;
}
// Interpolate along an arbitrary edge, typically one that may be on the
// volume boundary. This means careful computation of stuff requiring
// neighborhood information (e.g., gradients).
void InterpolateEdge(double value, T *s, float x[3], unsigned char edgeNum,
void InterpolateEdge(double value, T *s, int ijk[3], unsigned char edgeNum,
unsigned char edgeUses[4], vtkIdType *eIds);
// Produce the output points on the pixel axes for this pixel cell.
void GeneratePoints(double value, unsigned char loc, T *sPtr, float x[3],
void GeneratePoints(double value, unsigned char loc, T *sPtr, int ijk[3],
unsigned char *edgeUses, vtkIdType *eIds);
// Helper function to set up the point ids on pixel edges.
......@@ -321,7 +312,7 @@ vtkFlyingEdges2DAlgorithm():XCases(nullptr),EdgeMetaData(nullptr),Scalars(nullpt
// Interpolate a new point along a boundary edge. Make sure to consider
// proximity to boundary when computing gradients, etc.
template <class T> void vtkFlyingEdges2DAlgorithm<T>::
InterpolateEdge(double value, T *s, float x[3],unsigned char edgeNum,
InterpolateEdge(double value, T *s, int ijk[3], unsigned char edgeNum,
unsigned char edgeUses[12], vtkIdType *eIds)
{
// if this edge is not used then get out
......@@ -338,42 +329,42 @@ InterpolateEdge(double value, T *s, float x[3],unsigned char edgeNum,
const unsigned char *offsets = this->VertOffsets[vertMap[0]];
s0 = s + offsets[0]*this->Inc0 + offsets[1]*this->Inc1;
x0[0] = x[0] + offsets[0]*this->Spacing[this->Axis0];
x0[1] = x[1] + offsets[1]*this->Spacing[this->Axis1];
x0[0] = ijk[0] + offsets[0];
x0[1] = ijk[1] + offsets[1];
offsets = this->VertOffsets[vertMap[1]];
s1 = s + offsets[0]*this->Inc0 + offsets[1]*this->Inc1;
x1[0] = x[0] + offsets[0]*this->Spacing[this->Axis0];
x1[1] = x[1] + offsets[1]*this->Spacing[this->Axis1];
x1[0] = ijk[0] + offsets[0];
x1[1] = ijk[1] + offsets[1];
// Okay interpolate
double t = (value - *s0) / (*s1 - *s0);
float *xPtr = this->NewPoints + 3*vId;
xPtr[0] = x0[0] + t*(x1[0]-x0[0]);
xPtr[1] = x0[1] + t*(x1[1]-x0[1]);