Commit ff4b2116 authored by Bill Lorensen's avatar Bill Lorensen
Browse files

ENH: Added gradient computation.

parent 54af22f4
......@@ -79,10 +79,9 @@ void vtkMarchingCubes::GenerateValues(int numContours, float r1, float r2)
this->GenerateValues(numContours,rng);
}
void ComputePointNormal(int i, int j, int k, short *s, int dims[3],
int sliceSize, float origin[3], float ar[3], float n[3])
void ComputePointGradient(int i, int j, int k, short *s, int dims[3],
int sliceSize, float origin[3], float aspectRatio[3], float n[3])
{
static vtkMath math;
float sp, sm;
// x-direction
......@@ -90,19 +89,19 @@ void ComputePointNormal(int i, int j, int k, short *s, int dims[3],
{
sp = s[i+1 + j*dims[0] + k*sliceSize];
sm = s[i + j*dims[0] + k*sliceSize];
n[0] = (sp - sm) / ar[0];
n[0] = (sp - sm) / aspectRatio[0];
}
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] = (sp - sm) / ar[0];
n[0] = (sp - sm) / aspectRatio[0];
}
else
{
sp = s[i+1 + j*dims[0] + k*sliceSize];
sm = s[i-1 + j*dims[0] + k*sliceSize];
n[0] = 0.5 * (sp - sm) / ar[0];
n[0] = 0.5 * (sp - sm) / aspectRatio[0];
}
// y-direction
......@@ -110,19 +109,19 @@ void ComputePointNormal(int i, int j, int k, short *s, int dims[3],
{
sp = s[i + (j+1)*dims[0] + k*sliceSize];
sm = s[i + j*dims[0] + k*sliceSize];
n[1] = (sp - sm) / ar[1];
n[1] = (sp - sm) / aspectRatio[1];
}
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] = (sp - sm) / ar[1];
n[1] = (sp - sm) / aspectRatio[1];
}
else
{
sp = s[i + (j+1)*dims[0] + k*sliceSize];
sm = s[i + (j-1)*dims[0] + k*sliceSize];
n[1] = 0.5 * (sp - sm) / ar[1];
n[1] = 0.5 * (sp - sm) / aspectRatio[1];
}
// z-direction
......@@ -130,22 +129,21 @@ void ComputePointNormal(int i, int j, int k, short *s, int dims[3],
{
sp = s[i + j*dims[0] + (k+1)*sliceSize];
sm = s[i + j*dims[0] + k*sliceSize];
n[2] = (sp - sm) / ar[2];
n[2] = (sp - sm) / aspectRatio[2];
}
else if ( k == (dims[1]-1) )
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] = (sp - sm) / ar[2];
n[2] = (sp - sm) / aspectRatio[2];
}
else
{
sp = s[i + j*dims[0] + (k+1)*sliceSize];
sm = s[i + j*dims[0] + (k-1)*sliceSize];
n[2] = 0.5 * (sp - sm) / ar[2];
n[2] = 0.5 * (sp - sm) / aspectRatio[2];
}
math.Normalize(n);
}
//
......@@ -153,16 +151,18 @@ void ComputePointNormal(int i, int j, int k, short *s, int dims[3],
//
void vtkMarchingCubes::Execute()
{
static vtkMath math;
vtkFloatPoints *newPts;
vtkCellArray *newPolys;
vtkShortScalars *newScalars;
vtkFloatNormals *newNormals;
vtkFloatVectors *newGradients;
vtkStructuredPoints *input=(vtkStructuredPoints *)this->Input;
vtkPointData *pd=input->GetPointData();
vtkScalars *inScalars=pd->GetScalars();
short *scalars, s[8], value;
int dims[3];
float ar[3], origin[3];
float aspectRatio[3], origin[3];
int i, j, k, sliceSize;
static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
TRIANGLE_CASES *triCase;
......@@ -170,7 +170,7 @@ void vtkMarchingCubes::Execute()
int contNum, jOffset, kOffset, idx, ii, jj, index, *vert;
int ptIds[3];
float t, *x1, *x2, x[3], *n1, *n2, n[3];
float pts[8][3], normals[8][3];
float pts[8][3], gradients[8][3];
static int edges[12][2] = { {0,1}, {1,2}, {2,3}, {3,0},
{4,5}, {5,6}, {6,7}, {7,4},
{0,4}, {1,5}, {3,7}, {2,6}};
......@@ -193,7 +193,7 @@ void vtkMarchingCubes::Execute()
}
input->GetDimensions(dims);
input->GetOrigin(origin);
input->GetAspectRatio(ar);
input->GetAspectRatio(aspectRatio);
if ( strcmp("short",inScalars->GetDataType()) )
{
......@@ -205,10 +205,11 @@ void vtkMarchingCubes::Execute()
newPts = new vtkFloatPoints(10000,50000);
newScalars = new vtkShortScalars(10000,50000);
newNormals = new vtkFloatNormals(10000,50000);
newGradients = new vtkFloatVectors(10000,50000);
newPolys = new vtkCellArray();
newPolys->Allocate(newPolys->EstimateSize(25000,3));
//
// Traverse all voxel cells, generating triangles and point normals
// Traverse all voxel cells, generating triangles and point gradients
// using marching cubes algorithm.
//
sliceSize = dims[0] * dims[1];
......@@ -218,11 +219,11 @@ void vtkMarchingCubes::Execute()
for ( k=0; k < (dims[2]-1); k++)
{
kOffset = k*sliceSize;
pts[0][2] = origin[2] + k*ar[2];
pts[0][2] = origin[2] + k*aspectRatio[2];
for ( j=0; j < (dims[1]-1); j++)
{
jOffset = j*dims[0];
pts[0][1] = origin[1] + j*ar[1];
pts[0][1] = origin[1] + j*aspectRatio[1];
for ( i=0; i < (dims[0]-1); i++)
{
//get scalar values
......@@ -244,45 +245,45 @@ void vtkMarchingCubes::Execute()
if ( index == 0 || index == 255 ) continue; //no surface
//create voxel points
pts[0][0] = origin[0] + i*ar[0];
pts[0][0] = origin[0] + i*aspectRatio[0];
pts[1][0] = pts[0][0] + ar[0];
pts[1][0] = pts[0][0] + aspectRatio[0];
pts[1][1] = pts[0][1];
pts[1][2] = pts[0][2];
pts[2][0] = pts[0][0] + ar[0];
pts[2][1] = pts[0][1] + ar[1];
pts[2][0] = pts[0][0] + aspectRatio[0];
pts[2][1] = pts[0][1] + aspectRatio[1];
pts[2][2] = pts[0][2];
pts[3][0] = pts[0][0];
pts[3][1] = pts[0][1] + ar[1];
pts[3][1] = pts[0][1] + aspectRatio[1];
pts[3][2] = pts[0][2];
pts[4][0] = pts[0][0];
pts[4][1] = pts[0][1];
pts[4][2] = pts[0][2] + ar[2];
pts[4][2] = pts[0][2] + aspectRatio[2];
pts[5][0] = pts[0][0] + ar[0];
pts[5][0] = pts[0][0] + aspectRatio[0];
pts[5][1] = pts[0][1];
pts[5][2] = pts[0][2] + ar[2];
pts[5][2] = pts[0][2] + aspectRatio[2];
pts[6][0] = pts[0][0] + ar[0];
pts[6][1] = pts[0][1] + ar[1];
pts[6][2] = pts[0][2] + ar[2];
pts[6][0] = pts[0][0] + aspectRatio[0];
pts[6][1] = pts[0][1] + aspectRatio[1];
pts[6][2] = pts[0][2] + aspectRatio[2];
pts[7][0] = pts[0][0];
pts[7][1] = pts[0][1] + ar[1];
pts[7][2] = pts[0][2] + ar[2];
//create normals
ComputePointNormal(i,j,k, scalars, dims, sliceSize, origin, ar, normals[0]);
ComputePointNormal(i+1,j,k, scalars, dims, sliceSize, origin, ar, normals[1]);
ComputePointNormal(i+1,j+1,k, scalars, dims, sliceSize, origin, ar, normals[2]);
ComputePointNormal(i,j+1,k, scalars, dims, sliceSize, origin, ar, normals[3]);
ComputePointNormal(i,j,k+1, scalars, dims, sliceSize, origin, ar, normals[4]);
ComputePointNormal(i+1,j,k+1, scalars, dims, sliceSize, origin, ar, normals[5]);
ComputePointNormal(i+1,j+1,k+1, scalars, dims, sliceSize, origin, ar, normals[6]);
ComputePointNormal(i,j+1,k+1, scalars, dims, sliceSize, origin, ar, normals[7]);
pts[7][1] = pts[0][1] + aspectRatio[1];
pts[7][2] = pts[0][2] + aspectRatio[2];
//create gradients
ComputePointGradient(i,j,k, scalars, dims, sliceSize, origin, aspectRatio, gradients[0]);
ComputePointGradient(i+1,j,k, scalars, dims, sliceSize, origin, aspectRatio, gradients[1]);
ComputePointGradient(i+1,j+1,k, scalars, dims, sliceSize, origin, aspectRatio, gradients[2]);
ComputePointGradient(i,j+1,k, scalars, dims, sliceSize, origin, aspectRatio, gradients[3]);
ComputePointGradient(i,j,k+1, scalars, dims, sliceSize, origin, aspectRatio, gradients[4]);
ComputePointGradient(i+1,j,k+1, scalars, dims, sliceSize, origin, aspectRatio, gradients[5]);
ComputePointGradient(i+1,j+1,k+1, scalars, dims, sliceSize, origin, aspectRatio, gradients[6]);
ComputePointGradient(i,j+1,k+1, scalars, dims, sliceSize, origin, aspectRatio, gradients[7]);
triCase = triCases + index;
edge = triCase->edges;
......@@ -295,8 +296,8 @@ void vtkMarchingCubes::Execute()
t = (float)(value - s[vert[0]]) / (s[vert[1]] - s[vert[0]]);
x1 = pts[vert[0]];
x2 = pts[vert[1]];
n1 = normals[vert[0]];
n2 = normals[vert[1]];
n1 = gradients[vert[0]];
n2 = gradients[vert[1]];
for (jj=0; jj<3; jj++)
{
x[jj] = x1[jj] + t * (x2[jj] - x1[jj]);
......@@ -304,6 +305,8 @@ void vtkMarchingCubes::Execute()
}
ptIds[ii] = newPts->InsertNextPoint(x);
newScalars->InsertScalar(ptIds[ii],value);
newGradients->InsertVector(ptIds[ii],n);
math.Normalize(n);
newNormals->InsertNormal(ptIds[ii],n);
}
newPolys->InsertNextCell(3,ptIds);
......@@ -323,6 +326,7 @@ void vtkMarchingCubes::Execute()
this->SetPoints(newPts);
this->SetPolys(newPolys);
this->PointData.SetScalars(newScalars);
this->PointData.SetVectors(newGradients);
this->PointData.SetNormals(newNormals);
this->Squeeze();
}
......
Supports Markdown
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