Skip to content

GitLab

Commit 4e8b769f authored by Alexis Girault's avatar Alexis Girault
1 Browse files
Options

vtkMarchingCubes: Add support for image orientation

parent d17c2aab
Hide whitespace changes
Inline Side-by-side
Showing with 61 additions and 66 deletions
Filters/Core/vtkMarchingCubes.cxx
View file @ 4e8b769f
......@@ -18,6 +18,7 @@
#include "vtkCharArray.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkImageTransform.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
......@@ -83,7 +84,7 @@ vtkMTimeType vtkMarchingCubes::GetMTime()
// NOTE: We calculate the negative of the gradient for efficiency
template <class T>
void vtkMarchingCubesComputePointGradient(int i, int j, int k, T *s, int dims[3],
vtkIdType sliceSize, double spacing[3], double n[3])
vtkIdType sliceSize, double n[3])
{
double sp, sm;
......@@ -92,19 +93,19 @@ void vtkMarchingCubesComputePointGradient(int i, int j, int k, T *s, int dims[3]
{
sp = s[i+1 + j*dims[0] + k*sliceSize];
sm = s[i + j*dims[0] + k*sliceSize];
n[0] = (sm - sp) / spacing[0];
n[0] = sm - sp;
}
else if ( i == (dims[0]-1) )
{
sp = s[i + j*dims[0] + k*sliceSize];
sm = s[i-1 + j*dims[0] + k*sliceSize];
n[0] = (sm - sp) / spacing[0];
n[0] = sm - sp;
}
else
{
sp = s[i+1 + j*dims[0] + k*sliceSize];
sm = s[i-1 + j*dims[0] + k*sliceSize];
n[0] = 0.5 * (sm - sp) / spacing[0];
n[0] = 0.5 * (sm - sp);
}
// y-direction
......@@ -112,19 +113,19 @@ void vtkMarchingCubesComputePointGradient(int i, int j, int k, T *s, int dims[3]
{
sp = s[i + (j+1)*dims[0] + k*sliceSize];
sm = s[i + j*dims[0] + k*sliceSize];
n[1] = (sm - sp) / spacing[1];
n[1] = sm - sp;
}
else if ( j == (dims[1]-1) )
{
sp = s[i + j*dims[0] + k*sliceSize];
sm = s[i + (j-1)*dims[0] + k*sliceSize];
n[1] = (sm - sp) / spacing[1];
n[1] = sm - sp;
}
else
{
sp = s[i + (j+1)*dims[0] + k*sliceSize];
sm = s[i + (j-1)*dims[0] + k*sliceSize];
n[1] = 0.5 * (sm - sp) / spacing[1];
n[1] = 0.5 * (sm - sp);
}
// z-direction
......@@ -132,19 +133,19 @@ void vtkMarchingCubesComputePointGradient(int i, int j, int k, T *s, int dims[3]
{
sp = s[i + j*dims[0] + (k+1)*sliceSize];
sm = s[i + j*dims[0] + k*sliceSize];
n[2] = (sm - sp) / spacing[2];
n[2] = sm - sp;
}
else if ( k == (dims[2]-1) )
{
sp = s[i + j*dims[0] + k*sliceSize];
sm = s[i + j*dims[0] + (k-1)*sliceSize];
n[2] = (sm - sp) / spacing[2];
n[2] = sm - sp;
}
else
{
sp = s[i + j*dims[0] + (k+1)*sliceSize];
sm = s[i + j*dims[0] + (k-1)*sliceSize];
n[2] = 0.5 * (sm - sp) / spacing[2];
n[2] = 0.5 * (sm - sp);
}
}
......@@ -153,7 +154,6 @@ void vtkMarchingCubesComputePointGradient(int i, int j, int k, T *s, int dims[3]
//
template <class T>
void vtkMarchingCubesComputeGradient(vtkMarchingCubes *self,T *scalars, int dims[3],
double origin[3], double spacing[3],
vtkIncrementalPointLocator *locator,
vtkDataArray *newScalars,
vtkDataArray *newGradients,
......@@ -217,13 +217,13 @@ void vtkMarchingCubesComputeGradient(vtkMarchingCubes *self,T *scalars, int dims
break;
}
kOffset = k*sliceSize;
pts[0][2] = origin[2] + (k+extent[4]) * spacing[2];
zp = pts[0][2] + spacing[2];
pts[0][2] = k+extent[4];
zp = pts[0][2] + 1;
for ( j=0; j < (dims[1]-1); j++)
{
jOffset = j*dims[0];
pts[0][1] = origin[1] + (j+extent[2]) * spacing[1];
yp = pts[0][1] + spacing[1];
pts[0][1] = j+extent[2];
yp = pts[0][1] + 1;
for ( i=0; i < (dims[0]-1); i++)
{
//get scalar values
......@@ -246,8 +246,8 @@ void vtkMarchingCubesComputeGradient(vtkMarchingCubes *self,T *scalars, int dims
}
//create voxel points
pts[0][0] = origin[0] + (i+extent[0]) * spacing[0];
xp = pts[0][0] + spacing[0];
pts[0][0] = i+extent[0];
xp = pts[0][0] + 1;
pts[1][0] = xp;
pts[1][1] = pts[0][1];
......@@ -282,14 +282,14 @@ void vtkMarchingCubesComputeGradient(vtkMarchingCubes *self,T *scalars, int dims
//create gradients if needed
if (NeedGradients)
{
vtkMarchingCubesComputePointGradient(i,j,k, scalars, dims, sliceSize, spacing, gradients[0]);
vtkMarchingCubesComputePointGradient(i+1,j,k, scalars, dims, sliceSize, spacing, gradients[1]);
vtkMarchingCubesComputePointGradient(i+1,j+1,k, scalars, dims, sliceSize, spacing, gradients[2]);
vtkMarchingCubesComputePointGradient(i,j+1,k, scalars, dims, sliceSize, spacing, gradients[3]);
vtkMarchingCubesComputePointGradient(i,j,k+1, scalars, dims, sliceSize, spacing, gradients[4]);
vtkMarchingCubesComputePointGradient(i+1,j,k+1, scalars, dims, sliceSize, spacing, gradients[5]);
vtkMarchingCubesComputePointGradient(i+1,j+1,k+1, scalars, dims, sliceSize, spacing, gradients[6]);
vtkMarchingCubesComputePointGradient(i,j+1,k+1, scalars, dims, sliceSize, spacing, gradients[7]);
vtkMarchingCubesComputePointGradient(i,j,k, scalars, dims, sliceSize, gradients[0]);
vtkMarchingCubesComputePointGradient(i+1,j,k, scalars, dims, sliceSize, gradients[1]);
vtkMarchingCubesComputePointGradient(i+1,j+1,k, scalars, dims, sliceSize, gradients[2]);
vtkMarchingCubesComputePointGradient(i,j+1,k, scalars, dims, sliceSize, gradients[3]);
vtkMarchingCubesComputePointGradient(i,j,k+1, scalars, dims, sliceSize, gradients[4]);
vtkMarchingCubesComputePointGradient(i+1,j,k+1, scalars, dims, sliceSize, gradients[5]);
vtkMarchingCubesComputePointGradient(i+1,j+1,k+1, scalars, dims, sliceSize, gradients[6]);
vtkMarchingCubesComputePointGradient(i,j+1,k+1, scalars, dims, sliceSize, gradients[7]);
}
for (contNum=0; contNum < numValues; contNum++)
{
......@@ -389,7 +389,6 @@ int vtkMarchingCubes::RequestData(
vtkDataArray *inScalars;
int dims[3], extent[6];
vtkIdType estimatedSize;
double spacing[3], origin[3];
double bounds[6];
int numContours=this->ContourValues->GetNumberOfContours();
double *values=this->ContourValues->GetValues();
......@@ -418,8 +417,6 @@ int vtkMarchingCubes::RequestData(
return 1;
}
input->GetDimensions(dims);
input->GetOrigin(origin);
input->GetSpacing(spacing);
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
......@@ -436,8 +433,8 @@ int vtkMarchingCubes::RequestData(
// compute bounds for merging points
for ( int i=0; i<3; i++)
{
bounds[2*i] = origin[i] + extent[2*i] * spacing[i];
bounds[2*i+1] = origin[i] + extent[2*i+1] * spacing[i];
bounds[2*i] = extent[2*i];
bounds[2*i+1] = extent[2*i+1];
}
if ( this->Locator == nullptr )
{
......@@ -487,7 +484,7 @@ int vtkMarchingCubes::RequestData(
{
vtkTemplateMacro(
vtkMarchingCubesComputeGradient(this, static_cast<VTK_TT*>(scalars),
dims,origin,spacing,this->Locator,
dims,this->Locator,
newScalars,newGradients,
newNormals,newPolys,values,
numContours)
......@@ -504,7 +501,7 @@ int vtkMarchingCubes::RequestData(
inScalars->GetTuples(0,dataSize,image);
double *scalars = image->GetPointer(0);
vtkMarchingCubesComputeGradient(this,scalars,dims,origin,spacing,this->Locator,
vtkMarchingCubesComputeGradient(this,scalars,dims,this->Locator,
newScalars,newGradients,
newNormals,newPolys,values,numContours);
image->Delete();
......@@ -545,6 +542,8 @@ int vtkMarchingCubes::RequestData(
this->Locator->Initialize(); //free storage
}
vtkImageTransform::TransformPointSet(input, output);
return 1;
}
......
Filters/General/vtkDiscreteMarchingCubes.cxx
View file @ 4e8b769f
......@@ -19,6 +19,7 @@
#include "vtkCharArray.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkImageTransform.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
......@@ -60,7 +61,6 @@ vtkDiscreteMarchingCubes::~vtkDiscreteMarchingCubes() = default;
template <class T>
void vtkDiscreteMarchingCubesComputeGradient(
vtkDiscreteMarchingCubes *self,T *scalars, int dims[3],
double origin[3], double spacing[3],
vtkIncrementalPointLocator *locator,
vtkDataArray *newCellScalars,
vtkDataArray *newPointScalars,
......@@ -68,7 +68,7 @@ void vtkDiscreteMarchingCubesComputeGradient(
int numValues)
{
double s[8], value;
int i, j, k;
int i, j, k, pts[8][3], xp, yp, zp, *x1, *x2;
vtkIdType sliceSize, rowSize;
static const int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
vtkMarchingCubesTriangleCases *triCase, *triCases;
......@@ -79,8 +79,7 @@ void vtkDiscreteMarchingCubesComputeGradient(
int extent[6];
int ComputeScalars = newCellScalars != nullptr;
int ComputeAdjacentScalars = newPointScalars != nullptr;
double t, *x1, *x2, x[3], min, max;
double pts[8][3], xp, yp, zp;
double t, x[3], min, max;
static int edges[12][2] = { {0,1}, {1,2}, {3,2}, {0,3},
{4,5}, {5,6}, {7,6}, {4,7},
{0,4}, {1,5}, {3,7}, {2,6}};
......@@ -122,13 +121,13 @@ void vtkDiscreteMarchingCubesComputeGradient(
break;
}
kOffset = k*sliceSize;
pts[0][2] = origin[2] + (k+extent[4])*spacing[2];
zp = pts[0][2] + spacing[2];
pts[0][2] = k+extent[4];
zp = pts[0][2] + 1;
for ( j=0; j < (dims[1]-1); j++)
{
jOffset = j*rowSize;
pts[0][1] = origin[1] + (j+extent[2])*spacing[1];
yp = pts[0][1] + spacing[1];
pts[0][1] = j+extent[2];
yp = pts[0][1] + 1;
for ( i=0; i < (dims[0]-1); i++)
{
//get scalar values
......@@ -151,8 +150,8 @@ void vtkDiscreteMarchingCubesComputeGradient(
}
//create voxel points
pts[0][0] = origin[0] + (i+extent[0])*spacing[0];
xp = pts[0][0] + spacing[0];
pts[0][0] = i+extent[0];
xp = pts[0][0] + 1;
pts[1][0] = xp;
pts[1][1] = pts[0][1];
......@@ -277,7 +276,6 @@ int vtkDiscreteMarchingCubes::RequestData(
vtkDataArray *inScalars;
int dims[3], extent[6];
vtkIdType estimatedSize;
double spacing[3], origin[3];
double bounds[6];
vtkPolyData *output = vtkPolyData::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
......@@ -306,8 +304,6 @@ int vtkDiscreteMarchingCubes::RequestData(
return 1;
}
input->GetDimensions(dims);
input->GetOrigin(origin);
input->GetSpacing(spacing);
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);
......@@ -329,8 +325,8 @@ int vtkDiscreteMarchingCubes::RequestData(
// compute bounds for merging points
for ( int i=0; i<3; i++)
{
bounds[2*i] = origin[i] + extent[2*i] * spacing[i];
bounds[2*i+1] = origin[i] + extent[2*i+1] * spacing[i];
bounds[2*i] = extent[2*i];
bounds[2*i+1] = extent[2*i+1];
}
if ( this->Locator == nullptr )
{
......@@ -369,7 +365,7 @@ int vtkDiscreteMarchingCubes::RequestData(
vtkTemplateMacro(
vtkDiscreteMarchingCubesComputeGradient(this,
static_cast<VTK_TT*>(scalars),
dims, origin, spacing,
dims,
this->Locator, newCellScalars,
newPointScalars,
newPolys, values, numContours)
......@@ -391,8 +387,6 @@ int vtkDiscreteMarchingCubes::RequestData(
vtkDiscreteMarchingCubesComputeGradient(this,
scalars,
dims,
origin,
spacing,
this->Locator,
newCellScalars,
newPointScalars,
......@@ -432,5 +426,7 @@ int vtkDiscreteMarchingCubes::RequestData(
this->Locator->Initialize(); //free storage
}
vtkImageTransform::TransformPointSet(input, output);
return 1;
}
Filters/General/vtkImageMarchingCubes.cxx
View file @ 4e8b769f
......@@ -18,6 +18,7 @@
#include "vtkCommand.h"
#include "vtkFloatArray.h"
#include "vtkImageData.h"
#include "vtkImageTransform.h"
#include "vtkInformation.h"
#include "vtkInformationExecutivePortKey.h"
#include "vtkInformationVector.h"
......@@ -268,6 +269,8 @@ int vtkImageMarchingCubes::RequestData(
// Recover extra space.
output->Squeeze();
vtkImageTransform::TransformPointSet(inData,output);
// release the locators memory
this->DeleteLocator();
......@@ -336,7 +339,6 @@ int vtkImageMarchingCubesMakeNewPoint(vtkImageMarchingCubes *self,
int inc0, int inc1, int inc2,
T *ptr, int edge,
int *imageExtent,
double *spacing, double *origin,
double value)
{
int edgeAxis = 0;
......@@ -421,19 +423,19 @@ int vtkImageMarchingCubesMakeNewPoint(vtkImageMarchingCubes *self,
switch (edgeAxis)
{
case 0:
pt[0] = origin[0] + spacing[0] * ((double)idx0 + temp);
pt[1] = origin[1] + spacing[1] * ((double)idx1);
pt[2] = origin[2] + spacing[2] * ((double)idx2);
pt[0] = (double)idx0 + temp;
pt[1] = (double)idx1;
pt[2] = (double)idx2;
break;
case 1:
pt[0] = origin[0] + spacing[0] * ((double)idx0);
pt[1] = origin[1] + spacing[1] * ((double)idx1 + temp);
pt[2] = origin[2] + spacing[2] * ((double)idx2);
pt[0] = (double)idx0;
pt[1] = (double)idx1 + temp;
pt[2] = (double)idx2;
break;
case 2:
pt[0] = origin[0] + spacing[0] * ((double)idx0);
pt[1] = origin[1] + spacing[1] * ((double)idx1);
pt[2] = origin[2] + spacing[2] * ((double)idx2 + temp);
pt[0] = (double)idx0;
pt[1] = (double)idx1;
pt[2] = (double)idx2 + temp;
break;
}
......@@ -485,9 +487,9 @@ int vtkImageMarchingCubesMakeNewPoint(vtkImageMarchingCubes *self,
vtkImageMarchingCubesComputePointGradient(ptrB, gB, inc0, inc1, inc2,
b0, b1, b2);
// Interpolate Gradient
g[0] = (g[0] + temp * (gB[0] - g[0])) / spacing[0];
g[1] = (g[1] + temp * (gB[1] - g[1])) / spacing[1];
g[2] = (g[2] + temp * (gB[2] - g[2])) / spacing[2];
g[0] = g[0] + temp * (gB[0] - g[0]);
g[1] = g[1] + temp * (gB[1] - g[1]);
g[2] = g[2] + temp * (gB[2] - g[2]);
if (self->ComputeGradients)
{
self->Gradients->InsertNextTuple(g);
......@@ -578,8 +580,6 @@ void vtkImageMarchingCubesHandleCube(vtkImageMarchingCubes *self,
// If the point has not been created yet
if (pointIds[ii] == -1)
{
double *spacing = inData->GetSpacing();
double *origin = inData->GetOrigin();
int *extent =
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
......@@ -587,7 +587,7 @@ void vtkImageMarchingCubesHandleCube(vtkImageMarchingCubes *self,
cellX, cellY, cellZ,
inc0, inc1, inc2,
ptr, *edge, extent,
spacing, origin, value);
value);
self->AddLocatorPoint(cellX, cellY, *edge, pointIds[ii]);
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment