You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
5327 lines
172 KiB
5327 lines
172 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkVolumeShearWarpMapper.cxx,v $
|
|
|
|
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 "vtkVolumeShearWarpMapper.h"
|
|
#include "vtkVolumeShearWarpDataStructure.h"
|
|
|
|
#include "vtkVolume.h"
|
|
#include "vtkRenderer.h"
|
|
#include "vtkFiniteDifferenceGradientEstimator.h"
|
|
#include "vtkRenderWindow.h"
|
|
#include "vtkGraphicsFactory.h"
|
|
#include "vtkMatrix4x4.h"
|
|
#include "vtkTransform.h"
|
|
#include "vtkTimerLog.h"
|
|
#include "vtkPlaneCollection.h"
|
|
|
|
#include "vtkEncodedGradientShader.h"
|
|
#include "vtkEncodedGradientEstimator.h"
|
|
#include "vtkStructuredPoints.h"
|
|
#include "vtkCamera.h"
|
|
#include "vtkImageData.h"
|
|
#include "vtkVolumeProperty.h"
|
|
#include "vtkPiecewiseFunction.h"
|
|
|
|
#include "vtkOpenGLVolumeShearWarpMapper.h"
|
|
|
|
#define vtkVSWMultiplyPointMacro( A, B, M ) \
|
|
B[0] = A[0]*M[0] + A[1]*M[1] + A[2]*M[2] + M[3]; \
|
|
B[1] = A[0]*M[4] + A[1]*M[5] + A[2]*M[6] + M[7]; \
|
|
B[2] = A[0]*M[8] + A[1]*M[9] + A[2]*M[10] + M[11]; \
|
|
B[3] = A[0]*M[12] + A[1]*M[13] + A[2]*M[14] + M[15]; \
|
|
if ( B[3] != 1.0 ) { B[0] /= B[3]; B[1] /= B[3]; B[2] /= B[3]; }
|
|
|
|
#define vtkVSWMultiplyNormalMacro( A, B, M ) \
|
|
B[0] = A[0]*M[0] + A[1]*M[4] + A[2]*M[8]; \
|
|
B[1] = A[0]*M[1] + A[1]*M[5] + A[2]*M[9]; \
|
|
B[2] = A[0]*M[2] + A[1]*M[6] + A[2]*M[10]
|
|
|
|
vtkCxxRevisionMacro(vtkVolumeShearWarpMapper, "$Revision: 1.1 $");
|
|
vtkCxxRevisionMacro(vtkShearWarpBase,"$Revision: 1.1 $");
|
|
vtkCxxRevisionMacro(vtkShearWarpOctree<unsigned char>,"$Revision: 1.1 $");
|
|
vtkCxxRevisionMacro(vtkShearWarpOctree<unsigned short>,"$Revision: 1.1 $");
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Needed when we don't use the vtkStandardNewMacro.
|
|
vtkInstantiatorNewMacro(vtkVolumeShearWarpMapper);
|
|
//----------------------------------------------------------------------------
|
|
|
|
// Factor the view matrix into shear and warp
|
|
void vtkVolumeShearWarpMapper::FactorViewMatrix()
|
|
{
|
|
|
|
ComputeViewportMatrix();
|
|
ComputeViewMatrix();
|
|
|
|
if (this->ParallelProjection)
|
|
ComputePrincipalAxisParallel();
|
|
else if (this->MyPerspectiveProjection)
|
|
ComputePrincipalAxisParallel();
|
|
else
|
|
// ComputePrincipalAxisParallel();
|
|
ComputePrincipalAxisPerspective();
|
|
|
|
|
|
ComputePermutationMatrix();
|
|
|
|
if (this->ParallelProjection)
|
|
ComputeShearMatrixParallel();
|
|
else if (this->MyPerspectiveProjection)
|
|
ComputeShearMatrixParallel();
|
|
else
|
|
// ComputeShearMatrixParallel();
|
|
ComputeShearMatrixPerspective();
|
|
|
|
ComputeWarpMatrix();
|
|
}
|
|
|
|
|
|
// Compute the view matrix for parallel projection
|
|
void vtkVolumeShearWarpMapper::ComputeViewMatrix()
|
|
{
|
|
vtkTransform *view = vtkTransform::New();
|
|
view->SetMatrix(VoxelsToViewMatrix);
|
|
view->Inverse();
|
|
|
|
WorldViewingDirection[0] = 0.0;
|
|
WorldViewingDirection[1] = 0.0;
|
|
WorldViewingDirection[2] = 1.0;
|
|
WorldViewingDirection[3] = 0.0;
|
|
|
|
WorldEyePosition[0] = 0.0;
|
|
WorldEyePosition[1] = 0.0;
|
|
WorldEyePosition[2] = -1.0;
|
|
WorldEyePosition[3] = 0.0;
|
|
|
|
// Compute viewing direction in object space (for parallel projection)
|
|
view->MultiplyPoint(WorldViewingDirection,ObjectViewingDirection);
|
|
|
|
// Compute eye position in object space (for perspective projection)
|
|
view->MultiplyPoint(WorldEyePosition,ObjectEyePosition);
|
|
|
|
}
|
|
|
|
// Compute the viewport matrix
|
|
void vtkVolumeShearWarpMapper::ComputeViewportMatrix()
|
|
{
|
|
ViewportMatrix->Identity();
|
|
ViewportMatrix->Element[0][0] = 0.5 * (double) ImageViewportSize[0];
|
|
ViewportMatrix->Element[0][3] = 0.5 * (double) ImageViewportSize[0];
|
|
ViewportMatrix->Element[1][1] = 0.5 * (double) ImageViewportSize[1];
|
|
ViewportMatrix->Element[1][3] = 0.5 * (double) ImageViewportSize[1];
|
|
}
|
|
|
|
|
|
// Compute the principal viewing axis for parallel projection
|
|
void vtkVolumeShearWarpMapper::ComputePrincipalAxisParallel()
|
|
{
|
|
double x = fabs(ObjectViewingDirection[0]);
|
|
double y = fabs(ObjectViewingDirection[1]);
|
|
double z = fabs(ObjectViewingDirection[2]);
|
|
|
|
if (x >= y)
|
|
{
|
|
if (x >= z)
|
|
this->MajorAxis = VTK_X_AXIS;
|
|
else
|
|
this->MajorAxis = VTK_Z_AXIS;
|
|
}
|
|
else
|
|
{
|
|
if (y >= z)
|
|
this->MajorAxis = VTK_Y_AXIS;
|
|
else
|
|
this->MajorAxis = VTK_Z_AXIS;
|
|
|
|
}
|
|
|
|
if (ObjectViewingDirection[this->MajorAxis] > 0)
|
|
this->ReverseOrder = 0;
|
|
else
|
|
this->ReverseOrder = 1;
|
|
}
|
|
|
|
// Compute the principal viewing axis for perspective projection
|
|
void vtkVolumeShearWarpMapper::ComputePrincipalAxisPerspective()
|
|
{
|
|
double vertex[3];
|
|
double distance[3];
|
|
double eye[3];
|
|
double ax, ay, az;
|
|
double maximumDistance;
|
|
int maximumCount;
|
|
int count[3];
|
|
int order[3];
|
|
int axis[8];
|
|
int i;
|
|
|
|
eye[0] = ObjectEyePosition[0] / ObjectEyePosition[3];
|
|
eye[1] = ObjectEyePosition[1] / ObjectEyePosition[3];
|
|
eye[2] = ObjectEyePosition[2] / ObjectEyePosition[3];
|
|
|
|
// Find principal axes:
|
|
for (i=0; i<8; i++)
|
|
{
|
|
// Generate volume corners:
|
|
vertex[0] = -0.5 + (double)(i % 2);
|
|
vertex[1] = -0.5 + (double)((i/2) % 2);
|
|
vertex[2] = -0.5 + (double)((i/4) % 2);
|
|
|
|
vertex[0] *= this->GetInput()->GetDimensions()[0];
|
|
vertex[1] *= this->GetInput()->GetDimensions()[1];
|
|
vertex[2] *= this->GetInput()->GetDimensions()[2];
|
|
|
|
distance[0] = vertex[0] - eye[0];
|
|
distance[1] = vertex[1] - eye[1];
|
|
distance[2] = vertex[2] - eye[2];
|
|
|
|
// Determine the principal viewing axis and the stacking order:
|
|
ax = fabs(distance[0]);
|
|
|
|
ay = fabs(distance[1]);
|
|
az = fabs(distance[2]);
|
|
|
|
maximumDistance = ax;
|
|
|
|
if (ay > maximumDistance)
|
|
maximumDistance = ay;
|
|
|
|
if (az > maximumDistance)
|
|
maximumDistance = az;
|
|
|
|
if (maximumDistance == ax)
|
|
{
|
|
axis[i] = VTK_X_AXIS;
|
|
order[0] = (distance[0] < 0.0) ? 1 : 0;
|
|
}
|
|
else if (maximumDistance == ay)
|
|
{
|
|
axis[i] = VTK_Y_AXIS;
|
|
order[1] = (distance[1] < 0.0) ? 1 : 0;
|
|
}
|
|
else
|
|
{
|
|
axis[i] = VTK_Z_AXIS;
|
|
order[2] = (distance[2] < 0.0) ? 1 : 0;
|
|
}
|
|
}
|
|
|
|
// Find the dominating principal axis:
|
|
for (i=0; i<3; i++)
|
|
count[i] = 0;
|
|
|
|
for (i=0; i<8; i++)
|
|
{
|
|
switch (axis[i])
|
|
{
|
|
case VTK_X_AXIS:
|
|
count[0]++;
|
|
break;
|
|
|
|
case VTK_Y_AXIS:
|
|
count[1]++;
|
|
break;
|
|
|
|
case VTK_Z_AXIS:
|
|
count[2]++;
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Assign the dominant axis for the principal axis (favor the Z axis for ties):
|
|
maximumCount = count[0];
|
|
|
|
if (count[1] > maximumCount)
|
|
maximumCount = count[1];
|
|
|
|
if (count[2] > maximumCount)
|
|
maximumCount = count[2];
|
|
|
|
if (maximumCount == count[2])
|
|
{
|
|
this->MajorAxis = VTK_Z_AXIS;
|
|
this->ReverseOrder = order[2];
|
|
}
|
|
else if (maximumCount==count[1])
|
|
{
|
|
this->MajorAxis = VTK_Y_AXIS;
|
|
this->ReverseOrder = order[1];
|
|
}
|
|
else
|
|
{
|
|
this->MajorAxis = VTK_X_AXIS;
|
|
this->ReverseOrder = order[0];
|
|
}
|
|
}
|
|
|
|
|
|
// Compute the permutation matrix (transformation from object space to standard object space)
|
|
void vtkVolumeShearWarpMapper::ComputePermutationMatrix()
|
|
{
|
|
PermutationMatrix->Zero();
|
|
|
|
int size[3];
|
|
GetInput()->GetDimensions(size);
|
|
|
|
switch (this->MajorAxis)
|
|
|
|
{
|
|
case VTK_X_AXIS:
|
|
PermutationMatrix->Element[0][1] = 1.0f;
|
|
PermutationMatrix->Element[1][2] = 1.0f;
|
|
PermutationMatrix->Element[2][0] = 1.0f;
|
|
PermutationMatrix->Element[3][3] = 1.0f;
|
|
|
|
this->CountI = int(float(size[1]) / this->ImageSampleDistance);
|
|
this->CountJ = int(float(size[2]) / this->ImageSampleDistance);
|
|
this->CountK = int(float(size[0]) / this->ImageSampleDistance);
|
|
break;
|
|
|
|
case VTK_Y_AXIS:
|
|
PermutationMatrix->Element[0][2] = 1.0f;
|
|
|
|
PermutationMatrix->Element[1][0] = 1.0f;
|
|
PermutationMatrix->Element[2][1] = 1.0f;
|
|
PermutationMatrix->Element[3][3] = 1.0f;
|
|
|
|
this->CountI = int(float(size[2]) / this->ImageSampleDistance);
|
|
this->CountJ = int(float(size[0]) / this->ImageSampleDistance);
|
|
this->CountK = int(float(size[1]) / this->ImageSampleDistance);
|
|
break;
|
|
|
|
case VTK_Z_AXIS:
|
|
default:
|
|
PermutationMatrix->Element[0][0] = 1.0f;
|
|
PermutationMatrix->Element[1][1] = 1.0f;
|
|
PermutationMatrix->Element[2][2] = 1.0f;
|
|
PermutationMatrix->Element[3][3] = 1.0f;
|
|
|
|
this->CountI = int(float(size[0]) / this->ImageSampleDistance);
|
|
this->CountJ = int(float(size[1]) / this->ImageSampleDistance);
|
|
this->CountK = int(float(size[2]) / this->ImageSampleDistance);
|
|
break;
|
|
}
|
|
|
|
this->MaximumIntermediateDimension = size[0];
|
|
|
|
if (size[1] > this->MaximumIntermediateDimension)
|
|
this->MaximumIntermediateDimension = size[1];
|
|
|
|
if (size[2] > this->MaximumIntermediateDimension)
|
|
this->MaximumIntermediateDimension = size[2];
|
|
|
|
this->MaximumIntermediateDimension *= 2;
|
|
|
|
// Compute the viewing direction in standard object space (for parallel projection)
|
|
this->PermutationMatrix->MultiplyPoint(ObjectViewingDirection,StandardViewingDirection);
|
|
|
|
// Compute the eye position in standard object space (for perspective projection)
|
|
this->PermutationMatrix->MultiplyPoint(ObjectEyePosition,StandardEyePosition);
|
|
|
|
// Compute the permuted view to voxel matrix
|
|
vtkMatrix4x4::Multiply4x4(this->PermutationMatrix,this->ViewToVoxelsMatrix,this->PermutedViewToVoxelsMatrix);
|
|
|
|
// Compute the permuted voxel to view matrix
|
|
vtkMatrix4x4::Multiply4x4(this->PermutationMatrix,this->VoxelsToViewMatrix,this->PermutedVoxelsToViewMatrix);
|
|
|
|
|
|
// Get depth cueing factors
|
|
/*
|
|
this->DepthI = this->PermutedVoxelsToViewMatrix->Element[2][0];
|
|
this->DepthJ = this->PermutedVoxelsToViewMatrix->Element[2][1];
|
|
this->DepthK = this->PermutedVoxelsToViewMatrix->Element[2][2];
|
|
this->Depth0 = this->PermutedVoxelsToViewMatrix->Element[2][3];
|
|
*/
|
|
|
|
}
|
|
|
|
// Compute the shear matrix (transformation from object to intermediate image space)
|
|
void vtkVolumeShearWarpMapper::ComputeShearMatrixParallel()
|
|
{
|
|
vtkMatrix4x4 *conv = vtkMatrix4x4::New(); // conversion to intermediate image coordinate system
|
|
vtkMatrix4x4 *shear = vtkMatrix4x4::New(); // shear standard object space to intermediate image space
|
|
|
|
// Compute shear factors:
|
|
this->ShearI = - StandardViewingDirection[0] / StandardViewingDirection[2];
|
|
this->ShearJ = - StandardViewingDirection[1] / StandardViewingDirection[2];
|
|
this->Scale = 1.0f;
|
|
|
|
/* compute the intermediate image size */
|
|
this->IntermediateWidth = this->CountI + 1 + (int) ceil((this->CountK-1)*fabs(this->ShearI));
|
|
this->IntermediateHeight = this->CountJ + 1 + (int) ceil((this->CountK-1)*fabs(this->ShearJ));
|
|
|
|
/* compute the translation coefficients */
|
|
if (this->ShearI >= 0.0)
|
|
this->TranslationI = 1.0;
|
|
else
|
|
|
|
this->TranslationI = 1.0 - this->ShearI * (double) (this->CountK - 1);
|
|
|
|
if (this->ShearJ >= 0.0)
|
|
this->TranslationJ = 1.0;
|
|
else
|
|
this->TranslationJ = 1.0 - this->ShearJ * (double) (this->CountK - 1);
|
|
|
|
// Assemble standard object space shear matrix from shear factors
|
|
shear->Identity();
|
|
shear->Element[0][2] = this->ShearI;
|
|
shear->Element[1][2] = this->ShearJ;
|
|
|
|
// Add scale factor depending on object size:
|
|
// shear->Scale((double) size[0] / (double) this->CountI, (double) size[1] / (double) this->CountJ, (double) size[2] / (double) this->CountK);
|
|
|
|
// Create conversion matrix for intermediate image coordinates
|
|
conv->Identity();
|
|
conv->Element[0][3] = 0.5 * (double) this->IntermediateWidth;
|
|
conv->Element[1][3] = 0.5 * (double) this->IntermediateHeight;
|
|
|
|
vtkTransform *shearTransform = vtkTransform::New();
|
|
shearTransform->SetMatrix(this->PermutationMatrix);
|
|
shearTransform->PostMultiply();
|
|
shearTransform->Concatenate(shear);
|
|
shearTransform->Concatenate(conv);
|
|
this->ShearMatrix->DeepCopy(shearTransform->GetMatrix());
|
|
shearTransform->Delete();
|
|
|
|
shear->Delete();
|
|
conv->Delete();
|
|
|
|
}
|
|
|
|
// Compute the shear matrix (transformation from object to intermediate image space)
|
|
void vtkVolumeShearWarpMapper::ComputeShearMatrixPerspective()
|
|
{
|
|
|
|
vtkMatrix4x4 *conv = vtkMatrix4x4::New(); // conversion to intermediate image coordinate system
|
|
vtkMatrix4x4 *shear = vtkMatrix4x4::New(); // shear standard object space to intermediate image space
|
|
vtkMatrix4x4 *scale = vtkMatrix4x4::New();
|
|
|
|
// Compute shear factors: fabs
|
|
this->ShearI = - StandardEyePosition[0] / StandardEyePosition[2];
|
|
this->ShearJ = - StandardEyePosition[1] / StandardEyePosition[2];
|
|
this->Scale = - StandardEyePosition[3] / StandardEyePosition[2];
|
|
|
|
float scaleFactor = 1.0;
|
|
|
|
// if(this->Scale < 0.0)
|
|
// this->Scale = - this->Scale;
|
|
|
|
|
|
if (this->ReverseOrder==1)
|
|
scaleFactor = 1.0f / (1.0f - (this->CountK-1) * this->Scale);
|
|
|
|
/* compute the intermediate image size */
|
|
this->IntermediateWidth = this->CountI + 1 + (int) ceil((this->CountK-1)*fabs(this->ShearI));// - (int)(this->CountI - 1 - scaleFactor*(this->CountI-1));
|
|
this->IntermediateHeight = this->CountJ + 1 + (int) ceil((this->CountK-1)*fabs(this->ShearJ));// - (int)(this->CountJ - 1 - scaleFactor*(this->CountJ-1));
|
|
|
|
/* compute the translation coefficients */
|
|
|
|
if (this->ShearI >= 0.0)
|
|
this->TranslationI = 1.0;
|
|
else
|
|
this->TranslationI = 1.0 - this->ShearI * (double) (this->CountK - 1);
|
|
|
|
if (this->ShearJ >= 0.0)
|
|
this->TranslationJ = 1.0;
|
|
else
|
|
this->TranslationJ = 1.0 - this->ShearJ * (double) (this->CountK - 1);
|
|
|
|
// Assemble standard object space shear matrix from shear factors
|
|
shear->Identity();
|
|
shear->Element[0][2] = this->ShearI;
|
|
shear->Element[1][2] = this->ShearJ;
|
|
shear->Element[3][2] = this->Scale;
|
|
|
|
// Add scale factor depending on object size:
|
|
// shear->Scale((double) size[0] / (double) this->CountI, (double) size[1] / (double) this->CountJ, (double) size[2] / (double) this->CountK);
|
|
|
|
double sf;
|
|
|
|
if (this->ReverseOrder)
|
|
sf = 1.0 + this->Scale * (double)(this->CountK - 1);// - shear->Element[3][2];
|
|
else
|
|
sf = 1.0;//shear->Element[3][2] + 1;
|
|
|
|
/*
|
|
double sf;
|
|
|
|
if (this->ReverseOrder==1)
|
|
sf = 1.0 - shear->Element[3][2];
|
|
else
|
|
sf = shear->Element[3][2] + 1.0;
|
|
|
|
sf = 1.0 / sf; // invert scale factor
|
|
*/
|
|
scale->Identity();
|
|
scale->Element[0][0] = sf;
|
|
scale->Element[1][1] = sf;
|
|
|
|
// Create conversion matrix for intermediate image coordinates
|
|
conv->Identity();
|
|
conv->Element[0][3] = 0.5 * (double) this->IntermediateWidth;
|
|
conv->Element[1][3] = 0.5 * (double) this->IntermediateHeight;
|
|
|
|
vtkTransform *shearTransform = vtkTransform::New();
|
|
shearTransform->SetMatrix(this->PermutationMatrix);
|
|
shearTransform->PostMultiply();
|
|
// shearTransform->Concatenate(scale);
|
|
shearTransform->Concatenate(shear);
|
|
shearTransform->Concatenate(conv);
|
|
this->ShearMatrix->DeepCopy(shearTransform->GetMatrix());
|
|
shearTransform->Delete();
|
|
|
|
shear->Delete();
|
|
scale->Delete();
|
|
conv->Delete();
|
|
|
|
|
|
|
|
|
|
/*
|
|
double lb[4] = { 0, 0, this->CountK, 1.0 };
|
|
double rb[4] = { this->CountI, this->CountJ, this->CountK, 1.0 };
|
|
double tlt[4],tlb[4],trb[4];
|
|
shearTransform->MultiplyPoint(lb,tlb);
|
|
shearTransform->MultiplyPoint(rb,trb);
|
|
tlb[0] /= tlb[3];
|
|
tlb[1] /= tlb[3];
|
|
trb[0] /= trb[3];
|
|
trb[1] /= trb[3];
|
|
this->IntermediateWidth = trb[0] - tlb[0];
|
|
this->IntermediateHeight = trb[1] - tlb[1];
|
|
conv->Identity();
|
|
conv->Element[0][3] = 0.5 * (double) this->IntermediateWidth;
|
|
conv->Element[1][3] = 0.5 * (double) this->IntermediateHeight;
|
|
*/
|
|
|
|
|
|
}
|
|
|
|
// Compute a the two-dimensional warp matrix
|
|
void vtkVolumeShearWarpMapper::ComputeWarpMatrix()
|
|
{
|
|
|
|
vtkTransform *warp = vtkTransform::New();
|
|
/*
|
|
vtkTransform *inverse = vtkTransform::New();
|
|
inverse->SetMatrix(this->ShearMatrix);
|
|
inverse->Inverse();
|
|
|
|
// Invert shear matrix
|
|
warp->SetMatrix(VoxelsToViewMatrix);
|
|
|
|
// Compute warp matrix
|
|
warp->PreMultiply();
|
|
warp->Concatenate(inverse);
|
|
|
|
|
|
warp->PreMultiply();
|
|
warp->Concatenate(ViewportMatrix);
|
|
*/
|
|
|
|
// Compute inverse of shear matrix:
|
|
warp->SetMatrix(ShearMatrix);
|
|
warp->Inverse();
|
|
|
|
// Compute warp matrices:
|
|
warp->PostMultiply();
|
|
warp->Concatenate(VoxelsToViewMatrix);
|
|
warp->Concatenate(ViewportMatrix);
|
|
|
|
WarpMatrix->DeepCopy(warp->GetMatrix());
|
|
warp->Delete();
|
|
}
|
|
|
|
/*
|
|
// Compute a lookup table used for depth cueing
|
|
void vtkVolumeShearWarpMapper::ComputeDepthTable(int first, int last)
|
|
{
|
|
for (int i = first; i <= last; i++)
|
|
this->DepthTable[i] = this->FrontDepth * exp(-this->DepthDensity*(1.0 - i*this->DeltaDepth));
|
|
|
|
this->DepthTableSize = last - first;
|
|
}
|
|
|
|
|
|
// Perfrom second step for fast depth cueing (multiply intermediate image by second depth factor)
|
|
void vtkVolumeShearWarpMapper::DepthCueImage (vtkShearWarpPixelData *im, int slice)
|
|
{
|
|
|
|
double depthQuant = 1.0 / this->DeltaDepth;
|
|
|
|
float pixelDepth;
|
|
int pixelDepthInteger;
|
|
|
|
float uSlice = this->ShearI * slice + this->TranslationI;
|
|
float vSlice = this->ShearJ * slice + this->TranslationJ;
|
|
float leftDepth = this->Depth0 + this->DepthK*slice - uSlice*this->DepthI - vSlice*this->DepthJ;
|
|
|
|
if (this->DepthI > 0)
|
|
{
|
|
if (this->DepthJ > 0)
|
|
max_depth = left_depth + this->DepthI * width + this->DepthJ * height;
|
|
else
|
|
max_depth = left_depth + this->DepthI * width;
|
|
}
|
|
else
|
|
{
|
|
if (depth_dj > 0)
|
|
max_depth = left_depth + this->DepthJ * height;
|
|
else
|
|
max_depth = left_depth;
|
|
}
|
|
|
|
max_depth_int = max_depth * depth_quant;
|
|
|
|
if (max_depth_int >= vpc->dc_table_len)
|
|
{
|
|
// Resize table
|
|
}
|
|
|
|
float di = this->DepthI * depthQuant;
|
|
float dj = this->DepthJ * depthQuant;
|
|
|
|
leftDepth *= depthQuant;
|
|
|
|
for (int j = this->IntermediateHeight; j > 0; j--)
|
|
{
|
|
pixelDepth = leftDepth;
|
|
leftDepth += dj;
|
|
|
|
for (int i = this->IntermediateWidth; i > 0; i--)
|
|
{
|
|
pixelDepthInteger = pixelDepth;
|
|
pixelDepth += di;
|
|
|
|
if (pixelDepthInteger < 0)
|
|
pixelDepthInteger = 0;
|
|
|
|
if (pixelDepthInteger >= this->DepthTableSize)
|
|
{
|
|
// This shouldn't happen
|
|
pixelDepthInteger = this->DepthTableSize - 1;
|
|
}
|
|
|
|
im->Red *= this->DepthTable[pixelDepthInteger];
|
|
im->Green *= this->DepthTable[pixelDepthInteger];
|
|
im->Blue *= this->DepthTable[pixelDepthInteger];
|
|
im++;
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
|
|
// Simple parallel projection shear-warp without runlength encoded volume using bilinear interpolation
|
|
template <class T>
|
|
|
|
void CompositeIntermediateNearestSimple(vtkShearWarpRLEImage *image, vtkVolumeShearWarpMapper *myThis)
|
|
{
|
|
vtkShearWarpPixelData *pixels;
|
|
|
|
float sampledRed,sampledGreen,sampledBlue,sampledOpacity;
|
|
float sampledValue;
|
|
float sampledGradientMagnitude;
|
|
float oldRed,oldGreen,oldBlue,oldOpacity;
|
|
float newRed,newGreen,newBlue,newOpacity;
|
|
float redDiffuse,greenDiffuse,blueDiffuse;
|
|
float redSpecular,greenSpecular,blueSpecular;
|
|
float gradientOpacity;
|
|
float isoRed;
|
|
float isoGreen;
|
|
float isoBlue;
|
|
|
|
int skipped;
|
|
int i,j,k;
|
|
int vi,vj,vk;
|
|
|
|
float *SOTF;
|
|
float *CTF;
|
|
float *GTF;
|
|
float *GOTF;
|
|
float gradientOpacityConstant;
|
|
int gradientOpacityIsConstant;
|
|
unsigned short encodedNormal;
|
|
unsigned char gradientMagnitude;
|
|
|
|
float uSlice;
|
|
float vSlice;
|
|
|
|
int uSliceInteger;
|
|
int vSliceInteger;
|
|
|
|
T value;
|
|
|
|
int kStart;
|
|
int kEnd;
|
|
int kIncrement;
|
|
int vkStart;
|
|
int viIncrement;
|
|
int vjIncrement;
|
|
int vkIncrement;
|
|
|
|
T *dptr = (T*) myThis->GetInput()->GetScalarPointer();
|
|
|
|
unsigned short *nptr = myThis->GradientEstimator->GetEncodedNormals();
|
|
unsigned char *gptr = myThis->GradientEstimator->GetGradientMagnitudes();
|
|
|
|
int *dimensions = myThis->GetInput()->GetDimensions();
|
|
int location;
|
|
int plane = dimensions[0]*dimensions[1];
|
|
int halfDistance = myThis->ImageSampleDistance / 2;
|
|
|
|
//float depthCueFactor;
|
|
//float depthCueRatio;
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
kStart = myThis->CountK - 1;
|
|
kEnd = -1 + halfDistance;
|
|
kIncrement = -1;
|
|
}
|
|
else
|
|
{
|
|
kStart = 0;
|
|
kEnd = myThis->CountK - halfDistance;
|
|
kIncrement = 1;
|
|
}
|
|
|
|
SOTF = myThis->Volume->GetCorrectedScalarOpacityArray();
|
|
CTF = myThis->Volume->GetRGBArray();
|
|
GTF = myThis->Volume->GetGrayArray();
|
|
GOTF = myThis->Volume->GetGradientOpacityArray();
|
|
gradientOpacityConstant = myThis->Volume->GetGradientOpacityConstant();
|
|
|
|
if (gradientOpacityConstant > 0.0f)
|
|
gradientOpacityIsConstant = 1;
|
|
else
|
|
gradientOpacityIsConstant = 0;
|
|
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
isoRed = CTF[int(myThis->IsoValue)*3 + 0];
|
|
isoGreen = CTF[int(myThis->IsoValue)*3 + 1];
|
|
isoBlue = CTF[int(myThis->IsoValue)*3 + 2];
|
|
}
|
|
|
|
switch (myThis->MajorAxis)
|
|
{
|
|
case VTK_X_AXIS:
|
|
viIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vjIncrement = plane * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * myThis->ImageSampleDistance;
|
|
vkStart = kStart * myThis->ImageSampleDistance;
|
|
break;
|
|
|
|
case VTK_Y_AXIS:
|
|
viIncrement = plane * myThis->ImageSampleDistance;
|
|
vjIncrement = myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * dimensions[0] * myThis->ImageSampleDistance;
|
|
vkStart = kStart * dimensions[0] * myThis->ImageSampleDistance;
|
|
|
|
break;
|
|
|
|
|
|
|
|
case VTK_Z_AXIS:
|
|
default:
|
|
viIncrement = myThis->ImageSampleDistance;
|
|
vjIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * plane * myThis->ImageSampleDistance;
|
|
vkStart = kStart * plane * myThis->ImageSampleDistance;
|
|
|
|
break;
|
|
}
|
|
|
|
/*
|
|
float delta = myThis->DepthK - myThis->DepthI * myThis->ShearI - myThis->DepthJ * myThis->ShearJ;
|
|
|
|
cout << "\n\nDELTA: " << delta << "\n\n";
|
|
|
|
if (myThis->ReverseOrder == 1)
|
|
delta = -delta;
|
|
|
|
depthCueRatio = exp (myThis->DepthDensity * delta);
|
|
depthCueFactor = 1.0f;
|
|
*/
|
|
for (k = kStart, vk = vkStart; k != kEnd; k += kIncrement, vk += vkIncrement)
|
|
{
|
|
uSlice = k*myThis->ShearI + myThis->TranslationI;
|
|
vSlice = k*myThis->ShearJ + myThis->TranslationJ;
|
|
|
|
uSliceInteger = (int) ceil(uSlice) - 1;
|
|
vSliceInteger = (int) ceil(vSlice) - 1;
|
|
|
|
// Composite one slice into the intermediate image
|
|
for (j=0, vj = halfDistance; j < myThis->CountJ-halfDistance; j++, vj += vjIncrement)
|
|
{
|
|
image->Position(pixels,uSliceInteger + (vSliceInteger+j)*myThis->IntermediateWidth);
|
|
|
|
for (i=0, vi = halfDistance; i < myThis->CountI-halfDistance; )
|
|
{
|
|
// Skip opaque pixels in intermediate image
|
|
skipped = image->Skip(pixels);
|
|
|
|
// Update both runs if to be aligned with intermediate pixels
|
|
if (skipped > 0)
|
|
{
|
|
i += skipped;
|
|
vi += viIncrement * skipped;
|
|
}
|
|
else
|
|
{
|
|
if (myThis->IntermixIntersectingGeometry)
|
|
{
|
|
float depth = myThis->IntermediateZBuffer[myThis->ImageSampleDistance * (uSliceInteger + i) + myThis->ImageSampleDistance * (vSliceInteger + j) * myThis->IntermediateWidth * myThis->ImageSampleDistance];
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
if (k*myThis->ImageSampleDistance <= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
else
|
|
{
|
|
if (k*myThis->ImageSampleDistance >= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
|
|
}
|
|
|
|
// Only process non-opaque pixels
|
|
if (pixels->Offset == 0)
|
|
{
|
|
if (myThis->IsVoxelClipped(i*myThis->ImageSampleDistance,j*myThis->ImageSampleDistance,k*myThis->ImageSampleDistance) == 1)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
continue;
|
|
}
|
|
|
|
|
|
oldOpacity = pixels->Opacity;
|
|
oldRed = pixels->Red;
|
|
oldGreen = pixels->Green;
|
|
oldBlue = pixels->Blue;
|
|
|
|
location = vi + vj + vk;
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
sampledOpacity = 0.0f;
|
|
sampledRed = 0.0f;
|
|
sampledGreen = 0.0f;
|
|
sampledBlue = 0.0f;
|
|
|
|
value = dptr[location];
|
|
sampledOpacity += SOTF[value];
|
|
sampledRed += CTF[value*3 + 0];
|
|
sampledGreen += CTF[value*3 + 1];
|
|
sampledBlue += CTF[value*3 + 2];
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
sampledGradientMagnitude = 0.0f;
|
|
gradientOpacity = gradientOpacityConstant;
|
|
|
|
encodedNormal = nptr[location];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal];
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal];
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal];
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal];
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal];
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal];
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[location];
|
|
sampledGradientMagnitude += float(gradientMagnitude);
|
|
|
|
if (sampledGradientMagnitude > 255.0f)
|
|
gradientOpacity = GOTF[255];
|
|
else if (sampledGradientMagnitude < 0.0f)
|
|
gradientOpacity = GOTF[0];
|
|
else
|
|
gradientOpacity = GOTF[(unsigned char) sampledGradientMagnitude];
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
sampledOpacity *= gradientOpacity;
|
|
}
|
|
|
|
/*
|
|
if (1)
|
|
{
|
|
sampledRed *= depthCueFactor;
|
|
sampledGreen *= depthCueFactor;
|
|
sampledBlue *= depthCueFactor;
|
|
}
|
|
*/
|
|
|
|
// Alpha Compositing
|
|
newRed = oldRed + sampledOpacity * sampledRed * (1.0f - oldOpacity);
|
|
newGreen = oldGreen + sampledOpacity * sampledGreen * (1.0f - oldOpacity);
|
|
newBlue = oldBlue + sampledOpacity * sampledBlue * (1.0f - oldOpacity);
|
|
newOpacity = oldOpacity + sampledOpacity * (1.0f - oldOpacity);
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[location];
|
|
sampledValue += float(value);
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue > pixels->Value)
|
|
{
|
|
|
|
newRed = CTF[int(sampledValue)*3+0];
|
|
newGreen = CTF[int(sampledValue)*3+1];
|
|
newBlue = CTF[int(sampledValue)*3+2];
|
|
newOpacity = SOTF[int(sampledValue)];
|
|
pixels->Value = sampledValue;
|
|
}
|
|
else
|
|
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[location];
|
|
sampledValue += float(value);
|
|
|
|
if (sampledValue >= myThis->IsoValue)
|
|
{
|
|
sampledRed = isoRed;
|
|
sampledGreen = isoGreen;
|
|
|
|
sampledBlue = isoBlue;
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
|
|
encodedNormal = nptr[location];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal];
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal];
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal];
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal];
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal];
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal];
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
}
|
|
|
|
newRed = sampledRed;
|
|
newGreen = sampledGreen;
|
|
newBlue = sampledBlue;
|
|
newOpacity = 1.0f;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
|
|
pixels->Red = newRed;
|
|
pixels->Green = newGreen;
|
|
pixels->Blue = newBlue;
|
|
pixels->Opacity = newOpacity;
|
|
|
|
if (newOpacity >= 1.0f)
|
|
{
|
|
// The current intermediate pixel is opaque, so exit
|
|
// loop and skip opaque pixels
|
|
pixels->Offset = 1;
|
|
}
|
|
else
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
depthCueFactor *= (pow(depthCueRatio,myThis->ImageSampleDistance));
|
|
cout << depthCueFactor << ", ";
|
|
*/
|
|
}
|
|
|
|
//myThis->DepthCueImage(image->GetPixelData(),(kEnd-kIncrement) * myThis->ImageSampleDistance);
|
|
}
|
|
|
|
|
|
// Simple parallel projection shear-warp without runlength encoded volume using bilinear interpolation
|
|
template <class T>
|
|
void CompositeIntermediateLinearSimple(vtkShearWarpRLEImage *image, vtkVolumeShearWarpMapper *myThis)
|
|
{
|
|
vtkShearWarpPixelData *pixels = 0;
|
|
|
|
float sampledRed,sampledGreen,sampledBlue,sampledOpacity;
|
|
float sampledValue;
|
|
float sampledGradientMagnitude;
|
|
float oldRed,oldGreen,oldBlue,oldOpacity;
|
|
float newRed,newGreen,newBlue,newOpacity;
|
|
float redDiffuse,greenDiffuse,blueDiffuse;
|
|
float redSpecular,greenSpecular,blueSpecular;
|
|
float gradientOpacity;
|
|
|
|
|
|
float isoRed;
|
|
float isoGreen;
|
|
float isoBlue;
|
|
|
|
int skipped;
|
|
|
|
int i,j,k;
|
|
int vi,vj,vk;
|
|
|
|
float *SOTF;
|
|
float *CTF;
|
|
float *GTF;
|
|
float *GOTF;
|
|
float gradientOpacityConstant;
|
|
int gradientOpacityIsConstant;
|
|
unsigned short encodedNormal;
|
|
unsigned char gradientMagnitude;
|
|
|
|
float uSlice;
|
|
float vSlice;
|
|
|
|
int uSliceInteger;
|
|
int vSliceInteger;
|
|
|
|
float uSliceFractional;
|
|
float vSliceFractional;
|
|
|
|
|
|
float weightTopLeft;
|
|
float weightBottomLeft;
|
|
float weightTopRight;
|
|
float weightBottomRight;
|
|
|
|
float adjustedTL;
|
|
float adjustedTR;
|
|
float adjustedBL;
|
|
float adjustedBR;
|
|
|
|
T value;
|
|
|
|
int kStart;
|
|
int vkStart;
|
|
int kEnd;
|
|
int kIncrement;
|
|
int viIncrement;
|
|
int vjIncrement;
|
|
int vkIncrement;
|
|
|
|
T *dptr = (T*) myThis->GetInput()->GetScalarPointer();
|
|
|
|
unsigned short *nptr = myThis->GradientEstimator->GetEncodedNormals();
|
|
unsigned char *gptr = myThis->GradientEstimator->GetGradientMagnitudes();
|
|
|
|
int *dimensions = myThis->GetInput()->GetDimensions();
|
|
int locationTL,locationTR,locationBL,locationBR;
|
|
int plane = dimensions[0]*dimensions[1];
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
kStart = myThis->CountK - 1;
|
|
kEnd = -1;
|
|
kIncrement = -1;
|
|
}
|
|
else
|
|
{
|
|
kStart = 0;
|
|
kEnd = myThis->CountK;
|
|
kIncrement = 1;
|
|
}
|
|
|
|
SOTF = myThis->Volume->GetCorrectedScalarOpacityArray();
|
|
CTF = myThis->Volume->GetRGBArray();
|
|
GTF = myThis->Volume->GetGrayArray();
|
|
GOTF = myThis->Volume->GetGradientOpacityArray();
|
|
gradientOpacityConstant = myThis->Volume->GetGradientOpacityConstant();
|
|
|
|
if (gradientOpacityConstant > 0.0f)
|
|
|
|
gradientOpacityIsConstant = 1;
|
|
else
|
|
gradientOpacityIsConstant = 0;
|
|
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
isoRed = CTF[int(myThis->IsoValue)*3 + 0];
|
|
isoGreen = CTF[int(myThis->IsoValue)*3 + 1];
|
|
isoBlue = CTF[int(myThis->IsoValue)*3 + 2];
|
|
}
|
|
|
|
switch (myThis->MajorAxis)
|
|
{
|
|
case VTK_X_AXIS:
|
|
viIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vjIncrement = plane * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * myThis->ImageSampleDistance;
|
|
vkStart = kStart * myThis->ImageSampleDistance;
|
|
break;
|
|
|
|
case VTK_Y_AXIS:
|
|
viIncrement = plane * myThis->ImageSampleDistance;
|
|
vjIncrement = myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * dimensions[0] * myThis->ImageSampleDistance;
|
|
vkStart = kStart * dimensions[0] * myThis->ImageSampleDistance;
|
|
break;
|
|
|
|
|
|
case VTK_Z_AXIS:
|
|
default:
|
|
viIncrement = myThis->ImageSampleDistance;
|
|
vjIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * plane * myThis->ImageSampleDistance;
|
|
vkStart = kStart * plane * myThis->ImageSampleDistance;
|
|
break;
|
|
}
|
|
|
|
|
|
|
|
for (k = kStart, vk = vkStart; k != kEnd; k += kIncrement, vk += vkIncrement)
|
|
{
|
|
uSlice = k*myThis->ShearI + myThis->TranslationI;
|
|
vSlice = k*myThis->ShearJ + myThis->TranslationJ;
|
|
|
|
uSliceInteger = (int)ceil(uSlice) - 1;
|
|
vSliceInteger = (int)ceil(vSlice) - 1;
|
|
|
|
uSliceFractional = uSlice - uSliceInteger;
|
|
vSliceFractional = vSlice - vSliceInteger;
|
|
|
|
weightTopLeft = uSliceFractional * vSliceFractional;
|
|
weightBottomLeft = uSliceFractional * (1.0f - vSliceFractional);
|
|
|
|
weightTopRight = (1.0f - uSliceFractional) * vSliceFractional;
|
|
weightBottomRight = (1.0f - uSliceFractional) * (1.0f - vSliceFractional);
|
|
|
|
// Composite one slice into the intermediate image
|
|
for (j=0, vj = 0; j < myThis->CountJ; j++, vj += vjIncrement)
|
|
{
|
|
image->Position(pixels,uSliceInteger + (vSliceInteger+j)*myThis->IntermediateWidth);
|
|
|
|
for (i=0, vi = 0; i < myThis->CountI; )
|
|
{
|
|
// Skip opaque pixels in intermediate image
|
|
skipped = image->Skip(pixels);
|
|
|
|
// Update both runs if to be aligned with intermediate pixels
|
|
if (skipped > 0)
|
|
{
|
|
i += skipped;
|
|
vi += viIncrement/*myThis->ImageSampleDistance*/ * skipped;
|
|
}
|
|
else
|
|
{
|
|
if (myThis->IntermixIntersectingGeometry)
|
|
{
|
|
float depth = myThis->IntermediateZBuffer[myThis->ImageSampleDistance * (uSliceInteger + i) + myThis->ImageSampleDistance * (vSliceInteger + j) * myThis->IntermediateWidth * myThis->ImageSampleDistance];
|
|
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
if (k*myThis->ImageSampleDistance <= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
else
|
|
{
|
|
|
|
if (k*myThis->ImageSampleDistance >= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
}
|
|
|
|
// Only process non-opaque pixels
|
|
if (pixels->Offset == 0)
|
|
{
|
|
if (myThis->IsVoxelClipped(i*myThis->ImageSampleDistance,j*myThis->ImageSampleDistance,k*myThis->ImageSampleDistance) == 1)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
|
|
vi += viIncrement;
|
|
continue;
|
|
}
|
|
|
|
|
|
|
|
oldOpacity = pixels->Opacity;
|
|
oldRed = pixels->Red;
|
|
oldGreen = pixels->Green;
|
|
oldBlue = pixels->Blue;
|
|
|
|
locationTL = vi + vj + vk;
|
|
locationTR = locationTL + viIncrement;
|
|
locationBL = locationTL + vjIncrement;
|
|
locationBR = locationBL + viIncrement;
|
|
|
|
/* switch (myThis->MajorAxis)
|
|
{
|
|
case VTK_X_AXIS:
|
|
locationTL = vj * plane + vi *dimensions[0] + vk;
|
|
|
|
|
|
locationTR = vj * plane + (vi+myThis->ImageSampleDistance) *dimensions[0] + vk;
|
|
locationBL = (vj+myThis->ImageSampleDistance) * plane + vi *dimensions[0] + vk;
|
|
locationBR = (vj+myThis->ImageSampleDistance) * plane + vi *dimensions[0] + vk;
|
|
|
|
break;
|
|
|
|
case VTK_Y_AXIS:
|
|
locationTL = vi * plane + vk *dimensions[0] + vj;
|
|
locationTR = (vi+myThis->ImageSampleDistance) * plane + vk *dimensions[0] + vj;
|
|
locationBL = vi * plane + vk *dimensions[0] + (vj+myThis->ImageSampleDistance);
|
|
locationBR = (vi+myThis->ImageSampleDistance) * plane + vk *dimensions[0] + (vj+myThis->ImageSampleDistance);
|
|
break;
|
|
|
|
case VTK_Z_AXIS:
|
|
default:
|
|
locationTL = vk * plane + vj *dimensions[0] + vi;
|
|
locationTR = vk * plane + vj *dimensions[0] + (vi+myThis->ImageSampleDistance);
|
|
locationBL = vk * plane + (vj+myThis->ImageSampleDistance) *dimensions[0] + vi;
|
|
locationBR = vk * plane + (vj+myThis->ImageSampleDistance) *dimensions[0] + (vi+myThis->ImageSampleDistance);
|
|
break;
|
|
}*/
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
|
|
sampledOpacity = 0.0f;
|
|
sampledRed = 0.0f;
|
|
sampledGreen = 0.0f;
|
|
sampledBlue = 0.0f;
|
|
|
|
value = dptr[locationTL];
|
|
sampledOpacity += SOTF[value] * weightTopLeft;
|
|
sampledRed += CTF[value*3 + 0] * weightTopLeft;
|
|
sampledGreen += CTF[value*3 + 1] * weightTopLeft;
|
|
sampledBlue += CTF[value*3 + 2] * weightTopLeft;
|
|
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationTR];
|
|
sampledOpacity += SOTF[value] * weightTopRight;
|
|
sampledRed += CTF[value*3 + 0] * weightTopRight;
|
|
sampledGreen += CTF[value*3 + 1] * weightTopRight;
|
|
sampledBlue += CTF[value*3 + 2] * weightTopRight;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
value = dptr[locationBL];
|
|
sampledOpacity += SOTF[value] * weightBottomLeft;
|
|
sampledRed += CTF[value*3 + 0] * weightBottomLeft;
|
|
sampledGreen += CTF[value*3 + 1] * weightBottomLeft;
|
|
sampledBlue += CTF[value*3 + 2] * weightBottomLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationBR];
|
|
sampledOpacity += SOTF[value] * weightBottomRight;
|
|
sampledRed += CTF[value*3 + 0] * weightBottomRight;
|
|
sampledGreen += CTF[value*3 + 1] * weightBottomRight;
|
|
sampledBlue += CTF[value*3 + 2] * weightBottomRight;
|
|
}
|
|
}
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
sampledGradientMagnitude = 0.0f;
|
|
|
|
gradientOpacity = gradientOpacityConstant;
|
|
|
|
encodedNormal = nptr[locationTL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[locationTL];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightTopLeft;
|
|
}
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationTR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[locationTR];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightTopRight;
|
|
}
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
encodedNormal = nptr[locationBL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[locationBL];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightBottomLeft;
|
|
}
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationBR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[locationBR];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightBottomRight;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
if (sampledGradientMagnitude > 255.0f)
|
|
gradientOpacity = GOTF[255];
|
|
else if (sampledGradientMagnitude < 0.0f)
|
|
gradientOpacity = GOTF[0];
|
|
else
|
|
gradientOpacity = GOTF[(unsigned char) sampledGradientMagnitude];
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
sampledOpacity *= gradientOpacity;
|
|
}
|
|
|
|
// Alpha Compositing
|
|
newRed = oldRed + sampledOpacity * sampledRed * (1.0f - oldOpacity);
|
|
newGreen = oldGreen + sampledOpacity * sampledGreen * (1.0f - oldOpacity);
|
|
newBlue = oldBlue + sampledOpacity * sampledBlue * (1.0f - oldOpacity);
|
|
newOpacity = oldOpacity + sampledOpacity * (1.0f - oldOpacity);
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[locationTL];
|
|
sampledValue += float(value) * weightTopLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationTR];
|
|
sampledValue += float(value) * weightTopRight;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
value = dptr[locationBL];
|
|
sampledValue += float(value) * weightBottomLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationBR];
|
|
sampledValue += float(value) * weightBottomRight;
|
|
}
|
|
}
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue > pixels->Value)
|
|
{
|
|
newRed = CTF[int(sampledValue)*3+0];
|
|
newGreen = CTF[int(sampledValue)*3+1];
|
|
newBlue = CTF[int(sampledValue)*3+2];
|
|
newOpacity = SOTF[int(sampledValue)];
|
|
|
|
pixels->Value = sampledValue;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
sampledRed = isoRed;
|
|
sampledGreen = isoGreen;
|
|
sampledBlue = isoBlue;
|
|
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[locationTL];
|
|
sampledValue += float(value) * weightTopLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationTR];
|
|
sampledValue += float(value) * weightTopRight;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
value = dptr[locationBL];
|
|
sampledValue += float(value) * weightBottomLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
|
|
{
|
|
value = dptr[locationBR];
|
|
sampledValue += float(value) * weightBottomRight;
|
|
}
|
|
}
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue >= myThis->IsoValue)
|
|
{
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
|
|
adjustedTL = weightTopLeft;
|
|
adjustedBL = weightBottomLeft;
|
|
adjustedTR = weightTopRight;
|
|
adjustedBR = weightBottomRight;
|
|
|
|
if (i + 1 >= myThis->CountI)
|
|
{
|
|
adjustedTL += adjustedTR;
|
|
adjustedBL += adjustedBR;
|
|
}
|
|
|
|
if (j + 1 >= myThis->CountJ)
|
|
{
|
|
adjustedTL += adjustedBL;
|
|
adjustedTR += adjustedBR;
|
|
}
|
|
|
|
encodedNormal = nptr[locationTL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationTR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
encodedNormal = nptr[locationBL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationBR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
|
|
}
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
|
|
}
|
|
|
|
|
|
newRed = sampledRed;
|
|
newGreen = sampledGreen;
|
|
newBlue = sampledBlue;
|
|
newOpacity = 1.0f;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
|
|
|
|
pixels->Red = newRed;
|
|
pixels->Green = newGreen;
|
|
pixels->Blue = newBlue;
|
|
pixels->Opacity = newOpacity;
|
|
|
|
if (newOpacity >= 1.0f)
|
|
{
|
|
// The current intermediate pixel is opaque, so exit
|
|
|
|
// loop and skip opaque pixels
|
|
pixels->Offset = 1;
|
|
|
|
}
|
|
else
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
// Lacroute's parallel projection shear-warp algorithm with runlength encoded volume using nearest neighbour interpolation
|
|
template <class T>
|
|
void CompositeIntermediateNearestRLE(vtkShearWarpRLEImage *image, vtkVolumeShearWarpMapper *myThis)
|
|
{
|
|
vtkShearWarpRLESlice<T> *slice;
|
|
vtkShearWarpRLERun<T> *top;
|
|
vtkShearWarpPixelData *pixels;
|
|
|
|
float sampledRed,sampledGreen,sampledBlue,sampledOpacity;
|
|
float sampledValue;
|
|
float sampledGradientMagnitude;
|
|
float oldRed,oldGreen,oldBlue,oldOpacity;
|
|
|
|
float newRed,newGreen,newBlue,newOpacity;
|
|
float redDiffuse,greenDiffuse,blueDiffuse;
|
|
float redSpecular,greenSpecular,blueSpecular;
|
|
float gradientOpacity;
|
|
|
|
float isoRed;
|
|
float isoGreen;
|
|
float isoBlue;
|
|
|
|
int topIndex;
|
|
|
|
int skipped;
|
|
int runLength;
|
|
|
|
int i,j,k;
|
|
int vj,vk;
|
|
|
|
float *SOTF;
|
|
float *CTF;
|
|
float *GTF;
|
|
float *GOTF;
|
|
float gradientOpacityConstant;
|
|
int gradientOpacityIsConstant;
|
|
unsigned short encodedNormal;
|
|
unsigned char gradientMagnitude;
|
|
|
|
float uSlice;
|
|
float vSlice;
|
|
|
|
int uSliceInteger;
|
|
int vSliceInteger;
|
|
|
|
T value;
|
|
|
|
int kStart;
|
|
int kEnd;
|
|
int kIncrement;
|
|
int vkStart;
|
|
int vkIncrement;
|
|
|
|
int halfDistance = myThis->ImageSampleDistance / 2;
|
|
|
|
vtkShearWarpRLEVolume<T> *encodedVolume = dynamic_cast< vtkShearWarpRLEVolume<T> * > (myThis->EncodedVolume);
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
kStart = myThis->CountK - 1;
|
|
kEnd = -1 + halfDistance;
|
|
|
|
kIncrement = -1;
|
|
|
|
vkStart = (myThis->CountK - 1) * myThis->ImageSampleDistance - halfDistance;
|
|
vkIncrement = -myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
kStart = 0;
|
|
kEnd = myThis->CountK - halfDistance;
|
|
kIncrement = 1;
|
|
|
|
vkStart = halfDistance;
|
|
vkIncrement = myThis->ImageSampleDistance;
|
|
}
|
|
|
|
SOTF = myThis->Volume->GetCorrectedScalarOpacityArray();
|
|
CTF = myThis->Volume->GetRGBArray();
|
|
GTF = myThis->Volume->GetGrayArray();
|
|
GOTF = myThis->Volume->GetGradientOpacityArray();
|
|
gradientOpacityConstant = myThis->Volume->GetGradientOpacityConstant();
|
|
|
|
if (gradientOpacityConstant > 0.0f)
|
|
gradientOpacityIsConstant = 1;
|
|
else
|
|
gradientOpacityIsConstant = 0;
|
|
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
isoRed = CTF[int(myThis->IsoValue)*3 + 0];
|
|
isoGreen = CTF[int(myThis->IsoValue)*3 + 1];
|
|
isoBlue = CTF[int(myThis->IsoValue)*3 + 2];
|
|
}
|
|
|
|
for (k = kStart, vk = vkStart; k != kEnd; k += kIncrement, vk += vkIncrement)
|
|
{
|
|
uSlice = k*myThis->ShearI + myThis->TranslationI;
|
|
vSlice = k*myThis->ShearJ + myThis->TranslationJ;
|
|
|
|
uSliceInteger = (int)ceil(uSlice) - 1;
|
|
vSliceInteger = (int)ceil(vSlice) - 1;
|
|
|
|
slice = encodedVolume->GetSlice(myThis->MajorAxis,vk);
|
|
|
|
// Composite one slice into the intermediate image
|
|
for (j=0, vj = halfDistance; j < myThis->CountJ-halfDistance; j++, vj += myThis->ImageSampleDistance)
|
|
{
|
|
top = slice->GetLineRuns(vj);
|
|
topIndex = halfDistance;
|
|
|
|
while (topIndex >= top->Length)
|
|
|
|
{
|
|
topIndex -= top->Length;
|
|
top++;
|
|
}
|
|
|
|
image->Position(pixels,uSliceInteger + (vSliceInteger+j)*myThis->IntermediateWidth);
|
|
|
|
|
|
for (i=0; i < myThis->CountI; )
|
|
{
|
|
// Update run
|
|
while (topIndex >= top->Length)
|
|
{
|
|
|
|
topIndex -= top->Length;
|
|
top++;
|
|
}
|
|
|
|
// Skip opaque pixels in intermediate image
|
|
skipped = image->Skip(pixels);
|
|
|
|
// Update both runs if to be aligned with intermediate pixels
|
|
if (skipped > 0)
|
|
{
|
|
i += skipped;
|
|
topIndex += skipped * myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
runLength = top->Length - topIndex;
|
|
|
|
// Skip transparent voxels
|
|
|
|
|
|
if (top->VoxelData == NULL)
|
|
{
|
|
while (topIndex < top->Length)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
}
|
|
|
|
}
|
|
else
|
|
{
|
|
|
|
// This loops samples voxels, performs shading and performs
|
|
// compositing into the intermediate image
|
|
while (topIndex < top->Length)//h=0; h < runLength; h+=myThis->ImageSampleDistance)
|
|
{
|
|
if (myThis->IntermixIntersectingGeometry)
|
|
{
|
|
float depth = myThis->IntermediateZBuffer[myThis->ImageSampleDistance * (uSliceInteger + i) + myThis->ImageSampleDistance * (vSliceInteger + j) * myThis->IntermediateWidth * myThis->ImageSampleDistance];
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
if (vk <= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
else
|
|
{
|
|
|
|
if (vk >= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
}
|
|
|
|
// Only process non-opaque pixels
|
|
if (pixels->Offset == 0)
|
|
{
|
|
if (myThis->IsVoxelClipped(i*myThis->ImageSampleDistance,vj,vk) == 1)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
topIndex+=myThis->ImageSampleDistance;
|
|
continue;
|
|
}
|
|
|
|
oldOpacity = pixels->Opacity;
|
|
oldRed = pixels->Red;
|
|
oldGreen = pixels->Green;
|
|
oldBlue = pixels->Blue;
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
value = top->VoxelData[topIndex].Value;
|
|
sampledOpacity = SOTF[value];
|
|
sampledRed = CTF[value*3 + 0];
|
|
sampledGreen = CTF[value*3 + 1];
|
|
sampledBlue = CTF[value*3 + 2];
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
encodedNormal = top->VoxelData[topIndex].EncodedNormal;
|
|
|
|
redDiffuse = myThis->RedDiffuseShadingTable[encodedNormal];
|
|
redSpecular = myThis->RedSpecularShadingTable[encodedNormal];
|
|
greenDiffuse = myThis->GreenDiffuseShadingTable[encodedNormal];
|
|
greenSpecular = myThis->GreenSpecularShadingTable[encodedNormal];
|
|
blueDiffuse = myThis->BlueDiffuseShadingTable[encodedNormal];
|
|
blueSpecular = myThis->BlueSpecularShadingTable[encodedNormal];
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = top->VoxelData[topIndex].GradientMagnitude;
|
|
sampledGradientMagnitude = float(gradientMagnitude);
|
|
}
|
|
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
if (sampledGradientMagnitude > 255.0f)
|
|
gradientOpacity = GOTF[255];
|
|
else if (sampledGradientMagnitude < 0.0f)
|
|
gradientOpacity = GOTF[0];
|
|
else
|
|
gradientOpacity = GOTF[(unsigned char) sampledGradientMagnitude];
|
|
}
|
|
else
|
|
gradientOpacity = gradientOpacityConstant;
|
|
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
sampledOpacity *= gradientOpacity;
|
|
}
|
|
|
|
// Alpha Compositing
|
|
newRed = oldRed + sampledOpacity * sampledRed * (1.0f - oldOpacity);
|
|
newGreen = oldGreen + sampledOpacity * sampledGreen * (1.0f - oldOpacity);
|
|
newBlue = oldBlue + sampledOpacity * sampledBlue * (1.0f - oldOpacity);
|
|
newOpacity = oldOpacity + sampledOpacity * (1.0f - oldOpacity);
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
value = top->VoxelData[topIndex].Value;
|
|
sampledValue = value;
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue > pixels->Value)
|
|
{
|
|
newRed = CTF[int(sampledValue)*3+0];
|
|
newGreen = CTF[int(sampledValue)*3+1];
|
|
newBlue = CTF[int(sampledValue)*3+2];
|
|
newOpacity = SOTF[int(sampledValue)];
|
|
pixels->Value = sampledValue;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
sampledRed = isoRed;
|
|
sampledGreen = isoGreen;
|
|
sampledBlue = isoBlue;
|
|
sampledOpacity = 1.0f;
|
|
|
|
if (myThis->Shade)
|
|
|
|
{
|
|
encodedNormal = top->VoxelData[topIndex].EncodedNormal;
|
|
|
|
redDiffuse = myThis->RedDiffuseShadingTable[encodedNormal];
|
|
redSpecular = myThis->RedSpecularShadingTable[encodedNormal];
|
|
greenDiffuse = myThis->GreenDiffuseShadingTable[encodedNormal];
|
|
greenSpecular = myThis->GreenSpecularShadingTable[encodedNormal];
|
|
blueDiffuse = myThis->BlueDiffuseShadingTable[encodedNormal];
|
|
blueSpecular = myThis->BlueSpecularShadingTable[encodedNormal];
|
|
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
}
|
|
|
|
newRed = sampledRed;
|
|
newGreen = sampledGreen;
|
|
newBlue = sampledBlue;
|
|
newOpacity = 1.0f;
|
|
|
|
|
|
}
|
|
|
|
pixels->Red = newRed;
|
|
pixels->Green = newGreen;
|
|
pixels->Blue = newBlue;
|
|
pixels->Opacity = newOpacity;
|
|
|
|
if (newOpacity >= 1.0f)
|
|
{
|
|
// The current intermediate pixel is opaque, so exit
|
|
// loop and skip opaque pixels
|
|
pixels->Offset = 1;
|
|
break;
|
|
}
|
|
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
topIndex+=myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
// The current pixel has an offset greater than zero, so
|
|
// we exit the loop and skip all opaque pixels
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Lacroute's parallel projection shear-warp algorithm with runlength encoded volume using bilinear interpolation
|
|
template <class T>
|
|
void CompositeIntermediateLinearRLE(vtkShearWarpRLEImage *image, vtkVolumeShearWarpMapper *myThis)
|
|
{
|
|
vtkShearWarpRLESlice<T> *slice;
|
|
vtkShearWarpRLERun<T> *top;
|
|
vtkShearWarpRLERun<T> *bottom;
|
|
vtkShearWarpPixelData *pixels;
|
|
|
|
float sampledRed,sampledGreen,sampledBlue,sampledOpacity;
|
|
float sampledValue;
|
|
float sampledGradientMagnitude;
|
|
float oldRed,oldGreen,oldBlue,oldOpacity;
|
|
float newRed,newGreen,newBlue,newOpacity;
|
|
float redDiffuse,greenDiffuse,blueDiffuse;
|
|
float redSpecular,greenSpecular,blueSpecular;
|
|
float gradientOpacity;
|
|
|
|
float isoRed;
|
|
float isoGreen;
|
|
float isoBlue;
|
|
|
|
int topIndex;
|
|
int bottomIndex;
|
|
int topStart;
|
|
int bottomStart;
|
|
|
|
int skipped;
|
|
int runLength;
|
|
|
|
int i,j,k,h;
|
|
int vj,vk;
|
|
|
|
float *SOTF;
|
|
float *CTF;
|
|
float *GTF;
|
|
float *GOTF;
|
|
|
|
float gradientOpacityConstant;
|
|
int gradientOpacityIsConstant;
|
|
unsigned short encodedNormal;
|
|
unsigned char gradientMagnitude;
|
|
|
|
float uSlice;
|
|
float vSlice;
|
|
|
|
int uSliceInteger;
|
|
int vSliceInteger;
|
|
|
|
|
|
float uSliceFractional;
|
|
float vSliceFractional;
|
|
|
|
float weightTopLeft;
|
|
float weightBottomLeft;
|
|
float weightTopRight;
|
|
float weightBottomRight;
|
|
|
|
|
|
float adjustedTL;
|
|
float adjustedTR;
|
|
float adjustedBL;
|
|
float adjustedBR;
|
|
|
|
T value;
|
|
|
|
int kStart;
|
|
int kEnd;
|
|
int kIncrement;
|
|
int vkIncrement;
|
|
|
|
vtkShearWarpRLEVolume<T> *encodedVolume = dynamic_cast< vtkShearWarpRLEVolume<T> * > (myThis->EncodedVolume);
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
kStart = myThis->CountK - 1;
|
|
kEnd = -1;
|
|
kIncrement = -1;
|
|
vkIncrement = -myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
kStart = 0;
|
|
kEnd = myThis->CountK;
|
|
kIncrement = 1;
|
|
vkIncrement = myThis->ImageSampleDistance;
|
|
}
|
|
|
|
SOTF = myThis->Volume->GetCorrectedScalarOpacityArray();
|
|
CTF = myThis->Volume->GetRGBArray();
|
|
GTF = myThis->Volume->GetGrayArray();
|
|
GOTF = myThis->Volume->GetGradientOpacityArray();
|
|
gradientOpacityConstant = myThis->Volume->GetGradientOpacityConstant();
|
|
|
|
if (gradientOpacityConstant > 0.0f)
|
|
gradientOpacityIsConstant = 1;
|
|
else
|
|
gradientOpacityIsConstant = 0;
|
|
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
isoRed = CTF[int(myThis->IsoValue)*3 + 0];
|
|
isoGreen = CTF[int(myThis->IsoValue)*3 + 1];
|
|
isoBlue = CTF[int(myThis->IsoValue)*3 + 2];
|
|
}
|
|
|
|
for (k = kStart, vk = kStart*myThis->ImageSampleDistance; k != kEnd; k += kIncrement, vk += vkIncrement)
|
|
{
|
|
topIndex = 0;
|
|
bottomIndex = 0;
|
|
topStart = 0;
|
|
bottomStart = 0;
|
|
|
|
uSlice = k*myThis->ShearI + myThis->TranslationI;
|
|
vSlice = k*myThis->ShearJ + myThis->TranslationJ;
|
|
|
|
uSliceInteger = (int)ceil(uSlice) - 1;
|
|
vSliceInteger = (int)ceil(vSlice) - 1;
|
|
|
|
uSliceFractional = uSlice - uSliceInteger;
|
|
vSliceFractional = vSlice - vSliceInteger;
|
|
|
|
weightTopLeft = uSliceFractional * vSliceFractional;
|
|
weightBottomLeft = uSliceFractional * (1.0f - vSliceFractional);
|
|
weightTopRight = (1.0f - uSliceFractional) * vSliceFractional;
|
|
weightBottomRight = (1.0f - uSliceFractional) * (1.0f - vSliceFractional);
|
|
|
|
slice = encodedVolume->GetSlice(myThis->MajorAxis,vk);
|
|
|
|
// Composite one slice into the intermediate image
|
|
for (j=0, vj = 0; j < myThis->CountJ; j++, vj += myThis->ImageSampleDistance)
|
|
{
|
|
top = slice->GetLineRuns(vj);
|
|
|
|
if ((j + 1) < myThis->CountJ)
|
|
bottom = slice->GetLineRuns(vj+myThis->ImageSampleDistance);
|
|
else
|
|
bottom = NULL;
|
|
|
|
topStart = 0;
|
|
bottomStart = 0;
|
|
|
|
topIndex = 0;
|
|
bottomIndex = 0;
|
|
|
|
image->Position(pixels,uSliceInteger + (vSliceInteger+j)*myThis->IntermediateWidth);
|
|
|
|
for (i=0; i < myThis->CountI; )
|
|
{
|
|
// Update runs
|
|
while (topIndex >= top->Length)
|
|
{
|
|
topIndex -= top->Length;
|
|
topStart += top->Length;
|
|
top++;
|
|
}
|
|
|
|
if (bottom != NULL)
|
|
{
|
|
while (bottomIndex >= bottom->Length)
|
|
{
|
|
bottomIndex -= bottom->Length;
|
|
bottomStart += bottom->Length;
|
|
bottom++;
|
|
}
|
|
}
|
|
|
|
// Skip opaque pixels in intermediate image
|
|
skipped = image->Skip(pixels);
|
|
|
|
// Update both runs if to be aligned with intermediate pixels
|
|
if (skipped > 0)
|
|
{
|
|
i += skipped;
|
|
|
|
topIndex += skipped * myThis->ImageSampleDistance;
|
|
bottomIndex += skipped * myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
if (bottom != NULL)
|
|
{
|
|
if ( ((topStart + top->Length) - (topStart + topIndex)) < ((bottomStart + bottom->Length) - (bottomStart + bottomIndex)) )
|
|
runLength = ((topStart + top->Length) - (topStart + topIndex));
|
|
else
|
|
runLength = ((bottomStart + bottom->Length) - (bottomStart + bottomIndex));
|
|
}
|
|
else
|
|
runLength = top->Length - topIndex;
|
|
|
|
// Skip transparent voxels in both runs
|
|
|
|
if (top->VoxelData == NULL && (((bottom != NULL) ? bottom->VoxelData : NULL) == NULL))
|
|
|
|
{
|
|
for (h = 0; h < runLength; h+=myThis->ImageSampleDistance)
|
|
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
bottomIndex += myThis->ImageSampleDistance;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
|
|
// This loops samples voxels from both runs,
|
|
// interpolates, performs shading and performs
|
|
// compositing into the intermediate image
|
|
for (h=0; h < runLength; h+= myThis->ImageSampleDistance)
|
|
{
|
|
if (myThis->IntermixIntersectingGeometry)
|
|
{
|
|
float depth = myThis->IntermediateZBuffer[myThis->ImageSampleDistance * (uSliceInteger + i) + myThis->ImageSampleDistance * (vSliceInteger + j) * myThis->IntermediateWidth * myThis->ImageSampleDistance];
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
if (vk <= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
else
|
|
{
|
|
if (vk >= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
}
|
|
|
|
// Only process non-opaque pixels
|
|
if (pixels->Offset == 0)
|
|
{
|
|
|
|
if (myThis->IsVoxelClipped(i*myThis->ImageSampleDistance,vj,vk) == 1)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
bottomIndex += myThis->ImageSampleDistance;
|
|
continue;
|
|
}
|
|
|
|
|
|
oldOpacity = pixels->Opacity;
|
|
oldRed = pixels->Red;
|
|
oldGreen = pixels->Green;
|
|
oldBlue = pixels->Blue;
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
sampledOpacity = 0.0f;
|
|
sampledRed = 0.0f;
|
|
sampledGreen = 0.0f;
|
|
sampledBlue = 0.0f;
|
|
|
|
if (top->VoxelData != NULL)
|
|
{
|
|
value = top->VoxelData[topIndex].Value;
|
|
sampledOpacity += SOTF[value] * weightTopLeft;
|
|
sampledRed += CTF[value*3 + 0] * weightTopLeft;
|
|
sampledGreen += CTF[value*3 + 1] * weightTopLeft;
|
|
sampledBlue += CTF[value*3 + 2] * weightTopLeft;
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
value = top->VoxelData[topIndex+myThis->ImageSampleDistance].Value;
|
|
sampledOpacity += SOTF[value] * weightTopRight;
|
|
sampledRed += CTF[value*3 + 0] * weightTopRight;
|
|
sampledGreen += CTF[value*3 + 1] * weightTopRight;
|
|
sampledBlue += CTF[value*3 + 2] * weightTopRight;
|
|
}
|
|
}
|
|
|
|
if (bottom != NULL)
|
|
{
|
|
if (bottom->VoxelData != NULL)
|
|
{
|
|
value = bottom->VoxelData[bottomIndex].Value;
|
|
sampledOpacity += SOTF[value] * weightBottomLeft;
|
|
sampledRed += CTF[value*3 + 0] * weightBottomLeft;
|
|
sampledGreen += CTF[value*3 + 1] * weightBottomLeft;
|
|
sampledBlue += CTF[value*3 + 2] * weightBottomLeft;
|
|
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
value = bottom->VoxelData[bottomIndex+myThis->ImageSampleDistance].Value;
|
|
sampledOpacity += SOTF[value] * weightBottomRight;
|
|
sampledRed += CTF[value*3 + 0] * weightBottomRight;
|
|
sampledGreen += CTF[value*3 + 1] * weightBottomRight;
|
|
sampledBlue += CTF[value*3 + 2] * weightBottomRight;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
sampledGradientMagnitude = 0.0f;
|
|
gradientOpacity = gradientOpacityConstant;
|
|
|
|
if (top->VoxelData != NULL)
|
|
{
|
|
encodedNormal = top->VoxelData[topIndex].EncodedNormal;
|
|
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = top->VoxelData[topIndex].GradientMagnitude;
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightTopLeft;
|
|
}
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
encodedNormal = top->VoxelData[topIndex+myThis->ImageSampleDistance].EncodedNormal;
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = top->VoxelData[topIndex+myThis->ImageSampleDistance].GradientMagnitude;
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightTopRight;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (bottom != NULL)
|
|
{
|
|
if (bottom->VoxelData != NULL)
|
|
{
|
|
encodedNormal = bottom->VoxelData[bottomIndex].EncodedNormal;
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = bottom->VoxelData[bottomIndex].GradientMagnitude;
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightBottomLeft;
|
|
}
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
encodedNormal = bottom->VoxelData[bottomIndex+myThis->ImageSampleDistance].EncodedNormal;
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
|
|
gradientMagnitude = bottom->VoxelData[bottomIndex+myThis->ImageSampleDistance].GradientMagnitude;
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightBottomRight;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
if (sampledGradientMagnitude > 255.0f)
|
|
gradientOpacity = GOTF[255];
|
|
else if (sampledGradientMagnitude < 0.0f)
|
|
gradientOpacity = GOTF[0];
|
|
else
|
|
gradientOpacity = GOTF[(unsigned char) sampledGradientMagnitude];
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
sampledOpacity *= gradientOpacity;
|
|
}
|
|
|
|
// Alpha Compositing
|
|
newRed = oldRed + sampledOpacity * sampledRed * (1.0f - oldOpacity);
|
|
newGreen = oldGreen + sampledOpacity * sampledGreen * (1.0f - oldOpacity);
|
|
newBlue = oldBlue + sampledOpacity * sampledBlue * (1.0f - oldOpacity);
|
|
newOpacity = oldOpacity + sampledOpacity * (1.0f - oldOpacity);
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
sampledValue = 0.0f;
|
|
|
|
if (top->VoxelData != NULL)
|
|
{
|
|
value = top->VoxelData[topIndex].Value;
|
|
sampledValue += float(value) * weightTopLeft;
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
value = top->VoxelData[topIndex+myThis->ImageSampleDistance].Value;
|
|
sampledValue += float(value) * weightTopRight;
|
|
}
|
|
}
|
|
|
|
if (bottom != NULL)
|
|
{
|
|
if (bottom->VoxelData != NULL)
|
|
{
|
|
value = bottom->VoxelData[bottomIndex].Value;
|
|
sampledValue += float(value) * weightBottomLeft;
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
|
|
value = bottom->VoxelData[bottomIndex+myThis->ImageSampleDistance].Value;
|
|
sampledValue += float(value) * weightBottomRight;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue > pixels->Value)
|
|
{
|
|
newRed = CTF[int(sampledValue)*3+0];
|
|
newGreen = CTF[int(sampledValue)*3+1];
|
|
newBlue = CTF[int(sampledValue)*3+2];
|
|
newOpacity = SOTF[int(sampledValue)];
|
|
pixels->Value = sampledValue;
|
|
}
|
|
else
|
|
{
|
|
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
sampledRed = isoRed;
|
|
sampledGreen = isoGreen;
|
|
|
|
sampledBlue = isoBlue;
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
|
|
adjustedTL = weightTopLeft;
|
|
adjustedBL = weightBottomLeft;
|
|
adjustedTR = weightTopRight;
|
|
adjustedBR = weightBottomRight;
|
|
|
|
|
|
if (h + myThis->ImageSampleDistance >= runLength)
|
|
{
|
|
adjustedTL += adjustedTR;
|
|
|
|
adjustedBL += adjustedBR;
|
|
}
|
|
|
|
if (top->VoxelData == NULL)
|
|
{
|
|
adjustedBL += adjustedTL;
|
|
adjustedBR += adjustedTR;
|
|
}
|
|
else if ((((bottom != NULL) ? bottom->VoxelData : NULL) == NULL))
|
|
{
|
|
adjustedTL += adjustedBL;
|
|
adjustedTR += adjustedBR;
|
|
}
|
|
|
|
if (top->VoxelData != NULL)
|
|
{
|
|
encodedNormal = top->VoxelData[topIndex].EncodedNormal;
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
encodedNormal = top->VoxelData[topIndex+myThis->ImageSampleDistance].EncodedNormal;
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
}
|
|
}
|
|
|
|
if (bottom != NULL)
|
|
{
|
|
if (bottom->VoxelData != NULL)
|
|
{
|
|
encodedNormal = bottom->VoxelData[bottomIndex].EncodedNormal;
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
|
|
if (h + myThis->ImageSampleDistance < runLength)
|
|
{
|
|
|
|
encodedNormal = bottom->VoxelData[bottomIndex+myThis->ImageSampleDistance].EncodedNormal;
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
|
|
|
|
}
|
|
|
|
newRed = sampledRed;
|
|
newGreen = sampledGreen;
|
|
newBlue = sampledBlue;
|
|
newOpacity = 1.0f;
|
|
}
|
|
|
|
pixels->Red = newRed;
|
|
pixels->Green = newGreen;
|
|
pixels->Blue = newBlue;
|
|
pixels->Opacity = newOpacity;
|
|
|
|
|
|
if (newOpacity >= 1.0f)
|
|
{
|
|
// The current intermediate pixel is opaque, so exit
|
|
// loop and skip opaque pixels
|
|
pixels->Offset = 1;
|
|
break;
|
|
}
|
|
|
|
image->Advance(pixels,1);
|
|
|
|
i++;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
bottomIndex += myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
// The current pixel has an offset greater than zero, so
|
|
// we exit the loop and skip all opaque pixels
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Lacroute's perspective projection shear-warp algorithm with runlength encoded volume using bilinear interpolation
|
|
template <class T>
|
|
void CompositeIntermediateLinearRLEPerspective(vtkShearWarpRLEImage *image, vtkVolumeShearWarpMapper *myThis)
|
|
{
|
|
vtkShearWarpRLERun<T> **line;
|
|
float *lineIndex;
|
|
vtkShearWarpRLERun<T> *currentLine;
|
|
|
|
|
|
int currentLineIndex;
|
|
vtkShearWarpRLESlice<T> *slice;
|
|
vtkShearWarpPixelData *pixels;
|
|
int voxels;
|
|
int footprint;
|
|
int left;
|
|
|
|
float sampledRed,sampledGreen,sampledBlue,sampledOpacity;
|
|
float sampledValue;
|
|
float sampledGradientMagnitude;
|
|
float oldRed,oldGreen,oldBlue,oldOpacity;
|
|
float newRed,newGreen,newBlue,newOpacity;
|
|
|
|
float redDiffuse,greenDiffuse,blueDiffuse;
|
|
float redSpecular,greenSpecular,blueSpecular;
|
|
float gradientOpacity;
|
|
|
|
float isoRed;
|
|
float isoGreen;
|
|
|
|
float isoBlue;
|
|
|
|
int skipped;
|
|
int runLength;
|
|
|
|
int i,j,k,h,g;
|
|
int vk;
|
|
float vj;
|
|
|
|
float *SOTF;
|
|
float *CTF;
|
|
float *GTF;
|
|
float *GOTF;
|
|
float gradientOpacityConstant;
|
|
int gradientOpacityIsConstant;
|
|
unsigned short encodedNormal;
|
|
unsigned char gradientMagnitude;
|
|
|
|
float uSlice;
|
|
float vSlice;
|
|
|
|
int uSliceInteger;
|
|
int vSliceInteger;
|
|
|
|
float scaleFactor;
|
|
float weight;
|
|
|
|
T value;
|
|
|
|
|
|
int kStart;
|
|
int kEnd;
|
|
int kIncrement;
|
|
int vkIncrement;
|
|
|
|
vtkShearWarpRLEVolume<T> *encodedVolume = dynamic_cast< vtkShearWarpRLEVolume<T> * > (myThis->EncodedVolume);
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
kStart = myThis->CountK - 1;
|
|
kEnd = -1;
|
|
kIncrement = -1;
|
|
vkIncrement = -myThis->ImageSampleDistance;
|
|
|
|
}
|
|
else
|
|
{
|
|
kStart = 0;
|
|
kEnd = myThis->CountK;
|
|
kIncrement = 1;
|
|
vkIncrement = myThis->ImageSampleDistance;
|
|
}
|
|
|
|
SOTF = myThis->Volume->GetCorrectedScalarOpacityArray();
|
|
CTF = myThis->Volume->GetRGBArray();
|
|
|
|
GTF = myThis->Volume->GetGrayArray();
|
|
GOTF = myThis->Volume->GetGradientOpacityArray();
|
|
gradientOpacityConstant = myThis->Volume->GetGradientOpacityConstant();
|
|
|
|
if (gradientOpacityConstant > 0.0f)
|
|
gradientOpacityIsConstant = 1;
|
|
else
|
|
gradientOpacityIsConstant = 0;
|
|
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
|
|
{
|
|
isoRed = CTF[int(myThis->IsoValue)*3 + 0];
|
|
isoGreen = CTF[int(myThis->IsoValue)*3 + 1];
|
|
isoBlue = CTF[int(myThis->IsoValue)*3 + 2];
|
|
}
|
|
|
|
float vjIncrement;
|
|
float viIncrement;
|
|
|
|
|
|
for (k = kStart, vk = kStart*myThis->ImageSampleDistance; k != kEnd; k += kIncrement, vk += vkIncrement)
|
|
{
|
|
if (myThis->ReverseOrder)
|
|
scaleFactor = 1.0f / (1.0 - float(kStart*myThis->ImageSampleDistance-vk) * myThis->Scale);
|
|
else
|
|
scaleFactor = 1.0f / (1.0 + float(vk) * myThis->Scale);
|
|
|
|
footprint = (int) (1.0f + ceil ( 1.0f / scaleFactor));
|
|
|
|
vjIncrement = myThis->ImageSampleDistance / scaleFactor;
|
|
viIncrement = myThis->ImageSampleDistance / scaleFactor;
|
|
|
|
|
|
line = new vtkShearWarpRLERun<T>*[footprint];
|
|
lineIndex = new float[footprint];
|
|
|
|
uSlice = k*myThis->ShearI + myThis->TranslationI;
|
|
vSlice = k*myThis->ShearJ + myThis->TranslationJ;
|
|
|
|
uSliceInteger = (int)ceil(uSlice) - 1;
|
|
vSliceInteger = (int)ceil(vSlice) - 1;
|
|
|
|
|
|
|
|
slice = encodedVolume->GetSlice(myThis->MajorAxis,vk);
|
|
|
|
// Composite one slice into the intermediate image
|
|
for (j=0, vj = 0; j < myThis->CountJ*scaleFactor; j++, vj += vjIncrement)
|
|
{
|
|
for (g = 0; g < footprint; g++)
|
|
{
|
|
if (j+g < (myThis->CountJ*scaleFactor))
|
|
line[g] = slice->GetLineRuns(int(vj+g*myThis->ImageSampleDistance));
|
|
else
|
|
line[g] = NULL;
|
|
|
|
lineIndex[g] = 0;
|
|
}
|
|
|
|
image->Position(pixels,uSliceInteger + (vSliceInteger+j)*myThis->IntermediateWidth);
|
|
|
|
for (i=0; i < myThis->CountI*scaleFactor; )
|
|
{
|
|
// Skip opaque pixels in intermediate image
|
|
skipped = image->Skip(pixels);
|
|
|
|
// Update both runs if to be aligned with intermediate pixels
|
|
if (skipped > 0)
|
|
{
|
|
i += skipped;
|
|
|
|
for (g = 0; g < footprint; g++)
|
|
{
|
|
if (line[g] == NULL)
|
|
break;
|
|
|
|
|
|
lineIndex[g] += skipped * viIncrement;
|
|
|
|
while (lineIndex[g] >= line[g]->Length)
|
|
|
|
{
|
|
lineIndex[g] -= line[g]->Length;
|
|
line[g]++;
|
|
}
|
|
}
|
|
|
|
}
|
|
else
|
|
{
|
|
sampledOpacity = 0.0f;
|
|
sampledRed = 0.0f;
|
|
sampledGreen = 0.0f;
|
|
sampledBlue = 0.0f;
|
|
sampledValue = 0.0f;
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
sampledGradientMagnitude = 0.0f;
|
|
gradientOpacity = gradientOpacityConstant;
|
|
|
|
oldOpacity = pixels->Opacity;
|
|
oldRed = pixels->Red;
|
|
oldGreen = pixels->Green;
|
|
oldBlue = pixels->Blue;
|
|
|
|
voxels = 0;
|
|
|
|
// This loops samples voxels from both runs,
|
|
// interpolates, performs shading and performs
|
|
// compositing into the intermediate image
|
|
for (g = 0; g < footprint; g++)
|
|
{
|
|
|
|
weight = 1.0;
|
|
|
|
if (line[g] == NULL)
|
|
break;
|
|
|
|
left = footprint * myThis->ImageSampleDistance;
|
|
|
|
if (i + left >= myThis->CountI*scaleFactor)
|
|
left = (int)(myThis->CountI*scaleFactor - i) * myThis->ImageSampleDistance;
|
|
|
|
currentLine = line[g];
|
|
currentLineIndex = (int) lineIndex[g];
|
|
|
|
while (left > 0)
|
|
{
|
|
if (currentLineIndex >= currentLine->Length)
|
|
{
|
|
currentLineIndex -= currentLine->Length;
|
|
currentLine++;
|
|
}
|
|
|
|
runLength = (left < currentLine->Length) ? left : currentLine->Length;
|
|
left -= runLength;
|
|
|
|
if (currentLine->VoxelData == NULL)
|
|
{
|
|
currentLineIndex += runLength;
|
|
}
|
|
else
|
|
{
|
|
for (h = 0; h < runLength; h+= myThis->ImageSampleDistance)
|
|
{
|
|
voxels++;
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
value = currentLine->VoxelData[(int)currentLineIndex].Value;
|
|
if(value > 16000)
|
|
continue;
|
|
sampledOpacity += SOTF[value] * weight;
|
|
sampledRed += CTF[value*3 + 0] * weight;
|
|
sampledGreen += CTF[value*3 + 1] * weight;
|
|
sampledBlue += CTF[value*3 + 2] * weight;
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
encodedNormal = currentLine->VoxelData[(int)currentLineIndex].EncodedNormal;
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weight;
|
|
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weight;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weight;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weight;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = currentLine->VoxelData[(int)currentLineIndex].GradientMagnitude;
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weight;
|
|
}
|
|
}
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
value = currentLine->VoxelData[(int)currentLineIndex].Value;
|
|
sampledValue += float(value) * weight;
|
|
|
|
}
|
|
else
|
|
{
|
|
value = currentLine->VoxelData[(int)currentLineIndex].Value;
|
|
sampledValue += float(value) * weight;
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
|
|
encodedNormal = currentLine->VoxelData[(int)currentLineIndex].EncodedNormal;
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weight;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weight;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weight;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weight;
|
|
}
|
|
}
|
|
|
|
currentLineIndex += myThis->ImageSampleDistance;
|
|
}
|
|
}
|
|
}
|
|
|
|
lineIndex[g] += viIncrement;
|
|
|
|
if (lineIndex[g] >= line[g]->Length)
|
|
{
|
|
lineIndex[g] -= line[g]->Length;
|
|
|
|
line[g]++;
|
|
}
|
|
|
|
}
|
|
|
|
if (voxels > 0)
|
|
{
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
sampledRed /= (float)voxels;
|
|
sampledGreen /= (float)voxels;
|
|
sampledBlue /= (float)voxels;
|
|
sampledOpacity /= (float)voxels;
|
|
|
|
// Alpha Compositing
|
|
if (myThis->Shade)
|
|
|
|
{
|
|
sampledGradientMagnitude /= voxels;
|
|
|
|
if (sampledGradientMagnitude > 255.0f)
|
|
gradientOpacity = GOTF[255];
|
|
|
|
else if (sampledGradientMagnitude < 0.0f)
|
|
gradientOpacity = GOTF[0];
|
|
else
|
|
gradientOpacity = GOTF[(unsigned char) sampledGradientMagnitude];
|
|
|
|
redDiffuse /= (float)voxels;
|
|
redSpecular /= (float)voxels;
|
|
greenDiffuse /= (float)voxels;
|
|
greenSpecular /= (float)voxels;
|
|
blueDiffuse /= (float)voxels;
|
|
blueSpecular /= (float)voxels;
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
sampledOpacity *= gradientOpacity;
|
|
|
|
}
|
|
|
|
newRed = oldRed + sampledOpacity * sampledRed * (1.0f - oldOpacity);
|
|
newGreen = oldGreen + sampledOpacity * sampledGreen * (1.0f - oldOpacity);
|
|
newBlue = oldBlue + sampledOpacity * sampledBlue * (1.0f - oldOpacity);
|
|
newOpacity = oldOpacity + sampledOpacity * (1.0f - oldOpacity);
|
|
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
sampledValue /= (float)voxels;
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue > pixels->Value)
|
|
{
|
|
newRed = CTF[int(sampledValue)*3+0];
|
|
newGreen = CTF[int(sampledValue)*3+1];
|
|
newBlue = CTF[int(sampledValue)*3+2];
|
|
newOpacity = SOTF[int(sampledValue)];
|
|
pixels->Value = sampledValue;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
sampledValue /= (float)voxels;
|
|
|
|
// Isosurface
|
|
if (sampledValue > myThis->IsoValue)
|
|
{
|
|
sampledRed = isoRed;
|
|
sampledGreen = isoGreen;
|
|
sampledBlue = isoBlue;
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse /= (float)voxels;
|
|
redSpecular /= (float)voxels;
|
|
greenDiffuse /= (float)voxels;
|
|
greenSpecular /= (float)voxels;
|
|
blueDiffuse /= (float)voxels;
|
|
blueSpecular /= (float)voxels;
|
|
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
}
|
|
|
|
newRed = sampledRed;
|
|
newGreen = sampledGreen;
|
|
newBlue = sampledBlue;
|
|
newOpacity = 1.0f;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
|
|
pixels->Red = newRed;
|
|
pixels->Green = newGreen;
|
|
pixels->Blue = newBlue;
|
|
pixels->Opacity = newOpacity;
|
|
|
|
if (newOpacity >= 1.0f)
|
|
{
|
|
// The current intermediate pixel is opaque, so exit
|
|
// loop and skip opaque pixels
|
|
pixels->Offset = 1;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
}
|
|
}
|
|
}
|
|
|
|
delete[] lineIndex;
|
|
delete[] line;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// Parallel projection shear-warp fast classification algorithm using nearest neighbour interpolation
|
|
template <class T>
|
|
void CompositeIntermediateNearestUnclassified(vtkShearWarpRLEImage *image, vtkVolumeShearWarpMapper *myThis)
|
|
{
|
|
vtkShearWarpPixelData *pixels = 0;
|
|
vtkShearWarpOctreeRun *runs = new vtkShearWarpOctreeRun[myThis->CountJ];
|
|
vtkShearWarpOctreeRun *top;
|
|
|
|
int topIndex;
|
|
float sampledRed,sampledGreen,sampledBlue,sampledOpacity;
|
|
float sampledValue;
|
|
float sampledGradientMagnitude;
|
|
float oldRed,oldGreen,oldBlue,oldOpacity;
|
|
float newRed,newGreen,newBlue,newOpacity;
|
|
float redDiffuse,greenDiffuse,blueDiffuse;
|
|
float redSpecular,greenSpecular,blueSpecular;
|
|
float gradientOpacity;
|
|
|
|
float isoRed;
|
|
float isoGreen;
|
|
float isoBlue;
|
|
|
|
int skipped;
|
|
|
|
int i,j,k;
|
|
int vi,vj,vk;
|
|
|
|
float *SOTF;
|
|
float *CTF;
|
|
float *GTF;
|
|
float *GOTF;
|
|
float gradientOpacityConstant;
|
|
int gradientOpacityIsConstant;
|
|
unsigned short encodedNormal;
|
|
unsigned char gradientMagnitude;
|
|
|
|
float uSlice;
|
|
float vSlice;
|
|
|
|
int uSliceInteger;
|
|
int vSliceInteger;
|
|
|
|
T value;
|
|
|
|
int kStart;
|
|
int vkStart;
|
|
int kEnd;
|
|
int kIncrement;
|
|
int viIncrement;
|
|
int vjIncrement;
|
|
int vkIncrement;
|
|
int size;
|
|
|
|
T *dptr = (T*) myThis->GetInput()->GetScalarPointer();
|
|
|
|
unsigned short *nptr = myThis->GradientEstimator->GetEncodedNormals();
|
|
unsigned char *gptr = myThis->GradientEstimator->GetGradientMagnitudes();
|
|
|
|
int *dimensions = myThis->GetInput()->GetDimensions();
|
|
int location;
|
|
int plane = dimensions[0]*dimensions[1];
|
|
|
|
vtkShearWarpOctree<T> *octree = dynamic_cast< vtkShearWarpOctree<T> * > (myThis->Octree);
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
kStart = myThis->CountK - 1;
|
|
kEnd = -1;
|
|
kIncrement = -1;
|
|
}
|
|
else
|
|
{
|
|
kStart = 0;
|
|
kEnd = myThis->CountK;
|
|
kIncrement = 1;
|
|
}
|
|
|
|
SOTF = myThis->Volume->GetCorrectedScalarOpacityArray();
|
|
CTF = myThis->Volume->GetRGBArray();
|
|
GTF = myThis->Volume->GetGrayArray();
|
|
GOTF = myThis->Volume->GetGradientOpacityArray();
|
|
gradientOpacityConstant = myThis->Volume->GetGradientOpacityConstant();
|
|
|
|
if (gradientOpacityConstant > 0.0f)
|
|
gradientOpacityIsConstant = 1;
|
|
else
|
|
gradientOpacityIsConstant = 0;
|
|
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
isoRed = CTF[int(myThis->IsoValue)*3 + 0];
|
|
isoGreen = CTF[int(myThis->IsoValue)*3 + 1];
|
|
isoBlue = CTF[int(myThis->IsoValue)*3 + 2];
|
|
}
|
|
|
|
switch (myThis->MajorAxis)
|
|
{
|
|
case VTK_X_AXIS:
|
|
viIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vjIncrement = plane * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * myThis->ImageSampleDistance;
|
|
vkStart = kStart * myThis->ImageSampleDistance;
|
|
break;
|
|
|
|
case VTK_Y_AXIS:
|
|
viIncrement = plane * myThis->ImageSampleDistance;
|
|
vjIncrement = myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * dimensions[0] * myThis->ImageSampleDistance;
|
|
vkStart = kStart * dimensions[0] * myThis->ImageSampleDistance;
|
|
break;
|
|
|
|
case VTK_Z_AXIS:
|
|
default:
|
|
viIncrement = myThis->ImageSampleDistance;
|
|
vjIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * plane * myThis->ImageSampleDistance;
|
|
vkStart = kStart * plane * myThis->ImageSampleDistance;
|
|
break;
|
|
}
|
|
|
|
|
|
for (k = kStart, vk = vkStart; k != kEnd; k += kIncrement, vk += vkIncrement)
|
|
{
|
|
uSlice = k*myThis->ShearI + myThis->TranslationI;
|
|
vSlice = k*myThis->ShearJ + myThis->TranslationJ;
|
|
|
|
|
|
uSliceInteger = (int)ceil(uSlice) - 1;
|
|
vSliceInteger = (int)ceil(vSlice) - 1;
|
|
|
|
|
|
size = 0;
|
|
|
|
// Composite one slice into the intermediate image
|
|
for (j=0, vj = 0; j < myThis->CountJ; j++, vj += vjIncrement)
|
|
{
|
|
size -= 2 * myThis->ImageSampleDistance;
|
|
|
|
if (size <= 0)
|
|
size = octree->GetLineRuns(runs,myThis->MajorAxis, k*myThis->ImageSampleDistance, j*myThis->ImageSampleDistance);
|
|
|
|
top = runs;
|
|
topIndex = 0;
|
|
|
|
image->Position(pixels,uSliceInteger + (vSliceInteger+j)*myThis->IntermediateWidth);
|
|
|
|
for (i=0, vi = 0; i < myThis->CountI; )
|
|
{
|
|
// Update runs
|
|
while (topIndex >= top->Length)
|
|
{
|
|
topIndex -= top->Length;
|
|
top++;
|
|
}
|
|
|
|
// Skip opaque pixels in intermediate image
|
|
skipped = image->Skip(pixels);
|
|
|
|
// Update both runs if to be aligned with intermediate pixels
|
|
if (skipped > 0)
|
|
{
|
|
i += skipped;
|
|
vi += skipped * viIncrement;
|
|
topIndex += skipped * myThis->ImageSampleDistance;
|
|
}
|
|
|
|
else
|
|
{
|
|
// Skip transparent voxels
|
|
|
|
|
|
if (top->Type == VTK_SHEAR_WARP_OCTREE_TRANSPARENT)
|
|
{
|
|
while (topIndex < top->Length)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
|
|
// This loops samples voxels, performs shading and performs
|
|
// compositing into the intermediate image
|
|
while (topIndex < top->Length)
|
|
{
|
|
if (myThis->IntermixIntersectingGeometry)
|
|
{
|
|
float depth = myThis->IntermediateZBuffer[myThis->ImageSampleDistance * (uSliceInteger + i) + myThis->ImageSampleDistance * (vSliceInteger + j) * myThis->IntermediateWidth * myThis->ImageSampleDistance];
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
if (k*myThis->ImageSampleDistance <= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
else
|
|
{
|
|
if (k*myThis->ImageSampleDistance >= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
}
|
|
|
|
// Only process non-opaque pixels
|
|
if (pixels->Offset == 0)
|
|
{
|
|
if (myThis->IsVoxelClipped(i*myThis->ImageSampleDistance,j*myThis->ImageSampleDistance,k*myThis->ImageSampleDistance) == 1)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
continue;
|
|
}
|
|
|
|
oldOpacity = pixels->Opacity;
|
|
oldRed = pixels->Red;
|
|
oldGreen = pixels->Green;
|
|
oldBlue = pixels->Blue;
|
|
|
|
location = vi + vj + vk;
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
sampledOpacity = 0.0f;
|
|
sampledRed = 0.0f;
|
|
sampledGreen = 0.0f;
|
|
sampledBlue = 0.0f;
|
|
|
|
value = dptr[location];
|
|
sampledOpacity += SOTF[value];
|
|
sampledRed += CTF[value*3 + 0];
|
|
|
|
|
|
sampledGreen += CTF[value*3 + 1];
|
|
sampledBlue += CTF[value*3 + 2];
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
sampledGradientMagnitude = 0.0f;
|
|
gradientOpacity = gradientOpacityConstant;
|
|
|
|
encodedNormal = nptr[location];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal];
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal];
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal];
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal];
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal];
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal];
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[location];
|
|
sampledGradientMagnitude += float(gradientMagnitude);
|
|
|
|
if (sampledGradientMagnitude > 255.0f)
|
|
gradientOpacity = GOTF[255];
|
|
else if (sampledGradientMagnitude < 0.0f)
|
|
gradientOpacity = GOTF[0];
|
|
else
|
|
gradientOpacity = GOTF[(unsigned char) sampledGradientMagnitude];
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
sampledOpacity *= gradientOpacity;
|
|
|
|
}
|
|
|
|
// Alpha Compositing
|
|
newRed = oldRed + sampledOpacity * sampledRed * (1.0f - oldOpacity);
|
|
newGreen = oldGreen + sampledOpacity * sampledGreen * (1.0f - oldOpacity);
|
|
newBlue = oldBlue + sampledOpacity * sampledBlue * (1.0f - oldOpacity);
|
|
newOpacity = oldOpacity + sampledOpacity * (1.0f - oldOpacity);
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[location];
|
|
sampledValue += float(value);
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue > pixels->Value)
|
|
{
|
|
|
|
newRed = CTF[int(sampledValue)*3+0];
|
|
newGreen = CTF[int(sampledValue)*3+1];
|
|
newBlue = CTF[int(sampledValue)*3+2];
|
|
newOpacity = SOTF[int(sampledValue)];
|
|
pixels->Value = sampledValue;
|
|
}
|
|
else
|
|
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[location];
|
|
sampledValue += float(value);
|
|
|
|
|
|
if (sampledValue >= myThis->IsoValue)
|
|
{
|
|
sampledRed = isoRed;
|
|
|
|
sampledGreen = isoGreen;
|
|
sampledBlue = isoBlue;
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
|
|
encodedNormal = nptr[location];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal];
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal];
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal];
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal];
|
|
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal];
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal];
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
}
|
|
|
|
newRed = sampledRed;
|
|
newGreen = sampledGreen;
|
|
newBlue = sampledBlue;
|
|
newOpacity = 1.0f;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
|
|
pixels->Red = newRed;
|
|
|
|
pixels->Green = newGreen;
|
|
pixels->Blue = newBlue;
|
|
pixels->Opacity = newOpacity;
|
|
|
|
if (newOpacity >= 1.0f)
|
|
{
|
|
// The current intermediate pixel is opaque, so exit
|
|
// loop and skip opaque pixels
|
|
pixels->Offset = 1;
|
|
break;
|
|
}
|
|
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
// The current pixel has an offset greater than zero, so
|
|
// we exit the loop and skip all opaque pixels
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
delete[] runs;
|
|
}
|
|
|
|
|
|
// Parallel projection shear-warp fast classification algorithm using bilinear interpolation
|
|
template <class T>
|
|
void CompositeIntermediateLinearUnclassified(vtkShearWarpRLEImage *image, vtkVolumeShearWarpMapper *myThis)
|
|
{
|
|
vtkShearWarpPixelData *pixels = 0;
|
|
vtkShearWarpOctreeRun *runs = new vtkShearWarpOctreeRun[myThis->CountJ];
|
|
vtkShearWarpOctreeRun *top,*bottom;
|
|
|
|
int topIndex, bottomIndex, topStart, bottomStart;
|
|
|
|
float sampledRed,sampledGreen,sampledBlue,sampledOpacity;
|
|
float sampledValue;
|
|
float sampledGradientMagnitude;
|
|
float oldRed,oldGreen,oldBlue,oldOpacity;
|
|
float newRed,newGreen,newBlue,newOpacity;
|
|
float redDiffuse,greenDiffuse,blueDiffuse;
|
|
float redSpecular,greenSpecular,blueSpecular;
|
|
float gradientOpacity;
|
|
|
|
float isoRed;
|
|
float isoGreen;
|
|
|
|
float isoBlue;
|
|
|
|
int skipped;
|
|
int runLength;
|
|
|
|
int i,j,k,h;
|
|
int vi,vj,vk;
|
|
|
|
float *SOTF;
|
|
float *CTF;
|
|
float *GTF;
|
|
float *GOTF;
|
|
float gradientOpacityConstant;
|
|
int gradientOpacityIsConstant;
|
|
unsigned short encodedNormal;
|
|
unsigned char gradientMagnitude;
|
|
|
|
float uSlice;
|
|
float vSlice;
|
|
|
|
int uSliceInteger;
|
|
int vSliceInteger;
|
|
|
|
float uSliceFractional;
|
|
float vSliceFractional;
|
|
|
|
float weightTopLeft;
|
|
float weightBottomLeft;
|
|
float weightTopRight;
|
|
float weightBottomRight;
|
|
|
|
float adjustedTL;
|
|
float adjustedTR;
|
|
float adjustedBL;
|
|
float adjustedBR;
|
|
|
|
T value;
|
|
|
|
int kStart;
|
|
int vkStart;
|
|
int kEnd;
|
|
int kIncrement;
|
|
int viIncrement;
|
|
int vjIncrement;
|
|
int vkIncrement;
|
|
int size;
|
|
|
|
T *dptr = (T*) myThis->GetInput()->GetScalarPointer();
|
|
|
|
unsigned short *nptr = myThis->GradientEstimator->GetEncodedNormals();
|
|
unsigned char *gptr = myThis->GradientEstimator->GetGradientMagnitudes();
|
|
|
|
int *dimensions = myThis->GetInput()->GetDimensions();
|
|
int locationTL,locationTR,locationBL,locationBR;
|
|
int plane = dimensions[0]*dimensions[1];
|
|
|
|
vtkShearWarpOctree<T> *octree = dynamic_cast< vtkShearWarpOctree<T> * > (myThis->Octree);
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
kStart = myThis->CountK - 1;
|
|
kEnd = -1;
|
|
kIncrement = -1;
|
|
}
|
|
else
|
|
{
|
|
kStart = 0;
|
|
kEnd = myThis->CountK;
|
|
kIncrement = 1;
|
|
}
|
|
|
|
SOTF = myThis->Volume->GetCorrectedScalarOpacityArray();
|
|
CTF = myThis->Volume->GetRGBArray();
|
|
GTF = myThis->Volume->GetGrayArray();
|
|
GOTF = myThis->Volume->GetGradientOpacityArray();
|
|
gradientOpacityConstant = myThis->Volume->GetGradientOpacityConstant();
|
|
|
|
if (gradientOpacityConstant > 0.0f)
|
|
gradientOpacityIsConstant = 1;
|
|
else
|
|
gradientOpacityIsConstant = 0;
|
|
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
isoRed = CTF[int(myThis->IsoValue)*3 + 0];
|
|
isoGreen = CTF[int(myThis->IsoValue)*3 + 1];
|
|
isoBlue = CTF[int(myThis->IsoValue)*3 + 2];
|
|
}
|
|
|
|
switch (myThis->MajorAxis)
|
|
{
|
|
case VTK_X_AXIS:
|
|
viIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vjIncrement = plane * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * myThis->ImageSampleDistance;
|
|
vkStart = kStart * myThis->ImageSampleDistance;
|
|
break;
|
|
|
|
case VTK_Y_AXIS:
|
|
viIncrement = plane * myThis->ImageSampleDistance;
|
|
vjIncrement = myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * dimensions[0] * myThis->ImageSampleDistance;
|
|
vkStart = kStart * dimensions[0] * myThis->ImageSampleDistance;
|
|
break;
|
|
|
|
case VTK_Z_AXIS:
|
|
default:
|
|
viIncrement = myThis->ImageSampleDistance;
|
|
vjIncrement = dimensions[0] * myThis->ImageSampleDistance;
|
|
vkIncrement = kIncrement * plane * myThis->ImageSampleDistance;
|
|
vkStart = kStart * plane * myThis->ImageSampleDistance;
|
|
break;
|
|
}
|
|
|
|
|
|
for (k = kStart, vk = vkStart; k != kEnd; k += kIncrement, vk += vkIncrement)
|
|
{
|
|
uSlice = k*myThis->ShearI + myThis->TranslationI;
|
|
vSlice = k*myThis->ShearJ + myThis->TranslationJ;
|
|
|
|
uSliceInteger = (int)ceil(uSlice) - 1;
|
|
vSliceInteger = (int)ceil(vSlice) - 1;
|
|
|
|
uSliceFractional = uSlice - uSliceInteger;
|
|
vSliceFractional = vSlice - vSliceInteger;
|
|
|
|
weightTopLeft = uSliceFractional * vSliceFractional;
|
|
weightBottomLeft = uSliceFractional * (1.0f - vSliceFractional);
|
|
weightTopRight = (1.0f - uSliceFractional) * vSliceFractional;
|
|
weightBottomRight = (1.0f - uSliceFractional) * (1.0f - vSliceFractional);
|
|
|
|
size = 0;
|
|
|
|
// Composite one slice into the intermediate image
|
|
for (j=0, vj = 0; j < myThis->CountJ; j++, vj += vjIncrement)
|
|
{
|
|
size -= myThis->ImageSampleDistance;
|
|
|
|
if (size <= 0)
|
|
{
|
|
size = octree->GetLineRuns(runs,myThis->MajorAxis, k*myThis->ImageSampleDistance, j*myThis->ImageSampleDistance);
|
|
size -= myThis->ImageSampleDistance;
|
|
}
|
|
|
|
top = runs;
|
|
|
|
if ((j + 1) < myThis->CountJ)
|
|
{
|
|
size -= myThis->ImageSampleDistance;
|
|
|
|
if (size <= 0)
|
|
{
|
|
size = octree->GetLineRuns(runs, myThis->MajorAxis, k*myThis->ImageSampleDistance, (j+1)*myThis->ImageSampleDistance);
|
|
size -= myThis->ImageSampleDistance;
|
|
}
|
|
|
|
bottom = runs;
|
|
}
|
|
else
|
|
bottom = NULL;
|
|
|
|
topStart = 0;
|
|
bottomStart = 0;
|
|
|
|
topIndex = 0;
|
|
bottomIndex = 0;
|
|
|
|
image->Position(pixels,uSliceInteger + (vSliceInteger+j)*myThis->IntermediateWidth);
|
|
|
|
for (i=0, vi = 0; i < myThis->CountI; )
|
|
{
|
|
// Update runs
|
|
while (topIndex >= top->Length)
|
|
{
|
|
topIndex -= top->Length;
|
|
topStart += top->Length;
|
|
top++;
|
|
}
|
|
|
|
if (bottom != NULL)
|
|
{
|
|
while (bottomIndex >= bottom->Length)
|
|
{
|
|
bottomIndex -= bottom->Length;
|
|
bottomStart += bottom->Length;
|
|
bottom++;
|
|
}
|
|
}
|
|
|
|
// Skip opaque pixels in intermediate image
|
|
|
|
skipped = image->Skip(pixels);
|
|
|
|
// Update both runs if to be aligned with intermediate pixels
|
|
if (skipped > 0)
|
|
{
|
|
i += skipped;
|
|
vi += skipped * viIncrement;
|
|
topIndex += skipped * myThis->ImageSampleDistance;
|
|
bottomIndex += skipped * myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
if (bottom != NULL)
|
|
{
|
|
if ( ((topStart + top->Length) - (topStart + topIndex)) < ((bottomStart + bottom->Length) - (bottomStart + bottomIndex)) )
|
|
runLength = ((topStart + top->Length) - (topStart + topIndex));
|
|
else
|
|
runLength = ((bottomStart + bottom->Length) - (bottomStart + bottomIndex));
|
|
}
|
|
else
|
|
runLength = top->Length - topIndex;
|
|
|
|
// Skip transp1arent voxels in both runs
|
|
|
|
if (top->Type == VTK_SHEAR_WARP_OCTREE_TRANSPARENT && (((bottom != NULL) ? bottom->Type : VTK_SHEAR_WARP_OCTREE_TRANSPARENT) == VTK_SHEAR_WARP_OCTREE_TRANSPARENT))
|
|
{
|
|
for (h = 0; h < runLength; h+=myThis->ImageSampleDistance)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
bottomIndex += myThis->ImageSampleDistance;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// This loops samples voxels from both runs,
|
|
// interpolates, performs shading and performs
|
|
// compositing into the intermediate image
|
|
for (h=0; h < runLength; h+= myThis->ImageSampleDistance)
|
|
{
|
|
if (myThis->IntermixIntersectingGeometry)
|
|
{
|
|
float depth = myThis->IntermediateZBuffer[myThis->ImageSampleDistance * (uSliceInteger + i) + myThis->ImageSampleDistance * (vSliceInteger + j) * myThis->IntermediateWidth * myThis->ImageSampleDistance];
|
|
|
|
if (myThis->ReverseOrder)
|
|
{
|
|
if (k*myThis->ImageSampleDistance <= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
else
|
|
{
|
|
if (k*myThis->ImageSampleDistance >= depth)
|
|
pixels->Offset = 1;
|
|
}
|
|
|
|
}
|
|
|
|
// Only process non-opaque pixels
|
|
if (pixels->Offset == 0)
|
|
{
|
|
if (myThis->IsVoxelClipped(i*myThis->ImageSampleDistance,j*myThis->ImageSampleDistance,k*myThis->ImageSampleDistance) == 1)
|
|
{
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
|
|
topIndex += myThis->ImageSampleDistance;
|
|
bottomIndex += myThis->ImageSampleDistance;
|
|
continue;
|
|
}
|
|
|
|
oldOpacity = pixels->Opacity;
|
|
oldRed = pixels->Red;
|
|
oldGreen = pixels->Green;
|
|
oldBlue = pixels->Blue;
|
|
|
|
locationTL = vi + vj + vk;
|
|
locationTR = locationTL + viIncrement;
|
|
locationBL = locationTL + vjIncrement;
|
|
locationBR = locationBL + viIncrement;
|
|
|
|
if (myThis->FunctionType == VTK_SHEAR_WARP_COMPOSITE_FUNCTION)
|
|
{
|
|
sampledOpacity = 0.0f;
|
|
sampledRed = 0.0f;
|
|
sampledGreen = 0.0f;
|
|
sampledBlue = 0.0f;
|
|
|
|
value = dptr[locationTL];
|
|
sampledOpacity += SOTF[value] * weightTopLeft;
|
|
sampledRed += CTF[value*3 + 0] * weightTopLeft;
|
|
sampledGreen += CTF[value*3 + 1] * weightTopLeft;
|
|
sampledBlue += CTF[value*3 + 2] * weightTopLeft;
|
|
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationTR];
|
|
sampledOpacity += SOTF[value] * weightTopRight;
|
|
sampledRed += CTF[value*3 + 0] * weightTopRight;
|
|
sampledGreen += CTF[value*3 + 1] * weightTopRight;
|
|
sampledBlue += CTF[value*3 + 2] * weightTopRight;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
|
|
value = dptr[locationBL];
|
|
sampledOpacity += SOTF[value] * weightBottomLeft;
|
|
sampledRed += CTF[value*3 + 0] * weightBottomLeft;
|
|
sampledGreen += CTF[value*3 + 1] * weightBottomLeft;
|
|
sampledBlue += CTF[value*3 + 2] * weightBottomLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationBR];
|
|
sampledOpacity += SOTF[value] * weightBottomRight;
|
|
|
|
sampledRed += CTF[value*3 + 0] * weightBottomRight;
|
|
sampledGreen += CTF[value*3 + 1] * weightBottomRight;
|
|
sampledBlue += CTF[value*3 + 2] * weightBottomRight;
|
|
}
|
|
}
|
|
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
sampledGradientMagnitude = 0.0f;
|
|
gradientOpacity = gradientOpacityConstant;
|
|
|
|
encodedNormal = nptr[locationTL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightTopLeft;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightTopLeft;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[locationTL];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightTopLeft;
|
|
}
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationTR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightTopRight;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightTopRight;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[locationTR];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightTopRight;
|
|
}
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
encodedNormal = nptr[locationBL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightBottomLeft;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightBottomLeft;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
gradientMagnitude = gptr[locationBL];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightBottomLeft;
|
|
}
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationBR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * weightBottomRight;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * weightBottomRight;
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
|
|
gradientMagnitude = gptr[locationBR];
|
|
sampledGradientMagnitude += float(gradientMagnitude) * weightBottomRight;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (!gradientOpacityIsConstant)
|
|
{
|
|
if (sampledGradientMagnitude > 255.0f)
|
|
gradientOpacity = GOTF[255];
|
|
else if (sampledGradientMagnitude < 0.0f)
|
|
gradientOpacity = GOTF[0];
|
|
else
|
|
gradientOpacity = GOTF[(unsigned char) sampledGradientMagnitude];
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
sampledOpacity *= gradientOpacity;
|
|
}
|
|
|
|
// Alpha Compositing
|
|
newRed = oldRed + sampledOpacity * sampledRed * (1.0f - oldOpacity);
|
|
newGreen = oldGreen + sampledOpacity * sampledGreen * (1.0f - oldOpacity);
|
|
newBlue = oldBlue + sampledOpacity * sampledBlue * (1.0f - oldOpacity);
|
|
newOpacity = oldOpacity + sampledOpacity * (1.0f - oldOpacity);
|
|
}
|
|
else if (myThis->FunctionType == VTK_SHEAR_WARP_MIP_FUNCTION)
|
|
{
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[locationTL];
|
|
sampledValue += float(value) * weightTopLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationTR];
|
|
sampledValue += float(value) * weightTopRight;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
value = dptr[locationBL];
|
|
sampledValue += float(value) * weightBottomLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationBR];
|
|
sampledValue += float(value) * weightBottomRight;
|
|
|
|
}
|
|
}
|
|
|
|
// Maximum Intensity Projection
|
|
if (sampledValue > pixels->Value)
|
|
{
|
|
newRed = CTF[int(sampledValue)*3+0];
|
|
|
|
newGreen = CTF[int(sampledValue)*3+1];
|
|
newBlue = CTF[int(sampledValue)*3+2];
|
|
newOpacity = SOTF[int(sampledValue)];
|
|
pixels->Value = sampledValue;
|
|
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
|
|
sampledRed = isoRed;
|
|
sampledGreen = isoGreen;
|
|
sampledBlue = isoBlue;
|
|
|
|
sampledValue = 0.0f;
|
|
|
|
value = dptr[locationTL];
|
|
sampledValue += float(value) * weightTopLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationTR];
|
|
sampledValue += float(value) * weightTopRight;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
value = dptr[locationBL];
|
|
sampledValue += float(value) * weightBottomLeft;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
value = dptr[locationBR];
|
|
sampledValue += float(value) * weightBottomRight;
|
|
}
|
|
}
|
|
|
|
|
|
// Isosurface
|
|
if (sampledValue >= myThis->IsoValue)
|
|
{
|
|
if (myThis->Shade)
|
|
{
|
|
redDiffuse = 0.0f;
|
|
redSpecular = 0.0f;
|
|
blueDiffuse = 0.0f;
|
|
blueSpecular = 0.0f;
|
|
greenDiffuse = 0.0f;
|
|
greenSpecular = 0.0f;
|
|
|
|
adjustedTL = weightTopLeft;
|
|
adjustedBL = weightBottomLeft;
|
|
adjustedTR = weightTopRight;
|
|
adjustedBR = weightBottomRight;
|
|
|
|
if (i + 1 >= myThis->CountI)
|
|
{
|
|
adjustedTL += adjustedTR;
|
|
adjustedBL += adjustedBR;
|
|
}
|
|
|
|
if (j + 1 >= myThis->CountJ)
|
|
{
|
|
adjustedTL += adjustedBL;
|
|
adjustedTR += adjustedBR;
|
|
|
|
}
|
|
|
|
|
|
encodedNormal = nptr[locationTL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedTL;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedTL;
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationTR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedTR;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedTR;
|
|
}
|
|
|
|
if (j + 1 < myThis->CountJ)
|
|
{
|
|
encodedNormal = nptr[locationBL];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedBL;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedBL;
|
|
|
|
|
|
if (i + 1 < myThis->CountI)
|
|
{
|
|
encodedNormal = nptr[locationBR];
|
|
|
|
redDiffuse += myThis->RedDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
redSpecular += myThis->RedSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
greenDiffuse += myThis->GreenDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
greenSpecular += myThis->GreenSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
blueDiffuse += myThis->BlueDiffuseShadingTable[encodedNormal] * adjustedBR;
|
|
blueSpecular += myThis->BlueSpecularShadingTable[encodedNormal] * adjustedBR;
|
|
|
|
}
|
|
}
|
|
|
|
sampledRed *= redDiffuse + redSpecular;
|
|
sampledGreen *= greenDiffuse + greenSpecular;
|
|
sampledBlue *= blueDiffuse + blueSpecular;
|
|
|
|
|
|
}
|
|
|
|
|
|
newRed = sampledRed;
|
|
newGreen = sampledGreen;
|
|
newBlue = sampledBlue;
|
|
newOpacity = 1.0f;
|
|
}
|
|
else
|
|
{
|
|
newRed = oldRed;
|
|
newGreen = oldGreen;
|
|
newBlue = oldBlue;
|
|
newOpacity = oldOpacity;
|
|
}
|
|
|
|
}
|
|
|
|
pixels->Red = newRed;
|
|
pixels->Green = newGreen;
|
|
pixels->Blue = newBlue;
|
|
pixels->Opacity = newOpacity;
|
|
|
|
if (newOpacity >= 1.0f)
|
|
{
|
|
// The current intermediate pixel is opaque, so exit
|
|
// loop and skip opaque pixels
|
|
pixels->Offset = 1;
|
|
break;
|
|
}
|
|
|
|
|
|
image->Advance(pixels,1);
|
|
i++;
|
|
vi += viIncrement;
|
|
topIndex += myThis->ImageSampleDistance;
|
|
bottomIndex += myThis->ImageSampleDistance;
|
|
}
|
|
else
|
|
{
|
|
// The current pixel has an offset greater than zero, so
|
|
// we exit the loop and skip all opaque pixels
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
delete[] runs;
|
|
}
|
|
|
|
|
|
vtkVolumeShearWarpMapper* vtkVolumeShearWarpMapper::New()
|
|
{
|
|
// First try to create the object from the vtkObjectFactory
|
|
vtkObject* ret=vtkGraphicsFactory::CreateInstance("vtkVolumeShearWarpMapper");
|
|
return (vtkVolumeShearWarpMapper*)ret;
|
|
}
|
|
|
|
|
|
vtkVolumeShearWarpMapper::vtkVolumeShearWarpMapper()
|
|
{
|
|
this->ImageSampleDistance = 1;
|
|
this->MinimumImageSampleDistance = 1;
|
|
this->MaximumImageSampleDistance = 4;
|
|
this->AutoAdjustSampleDistances = 1;
|
|
|
|
this->PerspectiveMatrix = vtkMatrix4x4::New();
|
|
this->ViewToWorldMatrix = vtkMatrix4x4::New();
|
|
this->ViewToVoxelsMatrix = vtkMatrix4x4::New();
|
|
this->VoxelsToViewMatrix = vtkMatrix4x4::New();
|
|
this->WorldToVoxelsMatrix = vtkMatrix4x4::New();
|
|
|
|
this->VoxelsToWorldMatrix = vtkMatrix4x4::New();
|
|
this->VoxelTransformMatrix = vtkMatrix4x4::New();
|
|
this->ViewportMatrix = vtkMatrix4x4::New();
|
|
this->ShearMatrix = vtkMatrix4x4::New();
|
|
this->WarpMatrix = vtkMatrix4x4::New();
|
|
this->PermutationMatrix = vtkMatrix4x4::New();
|
|
this->PermutedViewToVoxelsMatrix = vtkMatrix4x4::New();
|
|
this->PermutedVoxelsToViewMatrix = vtkMatrix4x4::New();
|
|
this->PerspectiveTransform = vtkTransform::New();
|
|
|
|
this->EncodedVolume = NULL;
|
|
this->Octree = NULL;
|
|
|
|
this->ImageData = NULL;
|
|
this->ImageWidth = 0;
|
|
this->ImageHeight = 0;
|
|
this->AllocatedSize = 0;
|
|
|
|
this->IntemediateImage = NULL;
|
|
this->IntermediateWidth = 0;
|
|
this->IntermediateHeight = 0;
|
|
this->ScalarOpacityMTime = 0;
|
|
|
|
this->RenderTimeTable = NULL;
|
|
this->RenderVolumeTable = NULL;
|
|
this->RenderRendererTable = NULL;
|
|
this->RenderTableSize = 0;
|
|
this->RenderTableEntries = 0;
|
|
|
|
this->GradientEstimator = vtkFiniteDifferenceGradientEstimator::New();
|
|
this->GradientShader = vtkEncodedGradientShader::New();
|
|
|
|
this->IsoValue = 0.0f;
|
|
|
|
|
|
this->RunlengthEncoding = 0;
|
|
this->FastClassification = 0;
|
|
this->Shade = 0;
|
|
this->ParallelProjection = 0;
|
|
this->MyPerspectiveProjection = 0;
|
|
|
|
this->ZBuffer = NULL;
|
|
this->IntermediateZBuffer = NULL;
|
|
this->IntermixIntersectingGeometry = 0;
|
|
|
|
/*
|
|
this->FrontDepth = 1.0;
|
|
this->DepthDensity = 2.0;
|
|
this->DeltaDepth = 1.0 / 255.0;
|
|
this->DepthTableSize = 0;
|
|
|
|
ComputeDepthTable(0,511);
|
|
*/
|
|
this->Debug = 0;
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::SetGradientEstimator(
|
|
vtkEncodedGradientEstimator *gradest )
|
|
{
|
|
|
|
// If we are setting it to its current value, don't do anything
|
|
if ( this->GradientEstimator == gradest )
|
|
{
|
|
return;
|
|
}
|
|
|
|
// If we already have a gradient estimator, unregister it.
|
|
if ( this->GradientEstimator )
|
|
{
|
|
this->GradientEstimator->UnRegister(this);
|
|
this->GradientEstimator = NULL;
|
|
}
|
|
|
|
// If we are passing in a non-NULL estimator, register it
|
|
if ( gradest )
|
|
{
|
|
gradest->Register( this );
|
|
}
|
|
|
|
// Actually set the estimator, and consider the object Modified
|
|
this->GradientEstimator = gradest;
|
|
this->Modified();
|
|
}
|
|
|
|
vtkVolumeShearWarpMapper::~vtkVolumeShearWarpMapper()
|
|
{
|
|
|
|
this->PerspectiveMatrix->Delete();
|
|
this->ViewToWorldMatrix->Delete();
|
|
this->ViewToVoxelsMatrix->Delete();
|
|
this->VoxelsToViewMatrix->Delete();
|
|
this->WorldToVoxelsMatrix->Delete();
|
|
this->VoxelsToWorldMatrix->Delete();
|
|
this->VoxelTransformMatrix->Delete();
|
|
this->ViewportMatrix->Delete();
|
|
this->ShearMatrix->Delete();
|
|
this->WarpMatrix->Delete();
|
|
this->PermutationMatrix->Delete();
|
|
this->PermutedViewToVoxelsMatrix->Delete();
|
|
this->PerspectiveTransform->Delete();
|
|
|
|
this->SetGradientEstimator( NULL );
|
|
this->GradientShader->Delete();
|
|
|
|
if (this->EncodedVolume != NULL)
|
|
|
|
delete this->EncodedVolume;
|
|
|
|
if (this->IntemediateImage != NULL)
|
|
delete this->IntemediateImage;
|
|
|
|
if (this->ImageData != NULL)
|
|
delete[] this->ImageData;
|
|
|
|
if (this->Octree != NULL)
|
|
delete this->Octree;
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::Update()
|
|
{
|
|
if ( this->GetInput() )
|
|
{
|
|
this->GetInput()->UpdateInformation();
|
|
this->GetInput()->SetUpdateExtentToWholeExtent();
|
|
this->GetInput()->Update();
|
|
}
|
|
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::Render( vtkRenderer *ren,
|
|
vtkVolume *vol )
|
|
{
|
|
|
|
this->Volume = vol;
|
|
|
|
// make sure that we have scalar input and update the scalar input
|
|
if ( this->GetInput() == NULL )
|
|
{
|
|
vtkErrorMacro(<< "No Input!");
|
|
|
|
return;
|
|
}
|
|
else
|
|
{
|
|
this->GetInput()->UpdateInformation();
|
|
this->GetInput()->SetUpdateExtentToWholeExtent();
|
|
this->GetInput()->Update();
|
|
}
|
|
|
|
// Start timing now. We didn't want to capture the update of the
|
|
// input data in the times
|
|
this->Timer->StartTimer();
|
|
|
|
vol->UpdateTransferFunctions( ren );
|
|
|
|
ren->ComputeAspect();
|
|
double *aspect = ren->GetAspect();
|
|
vtkCamera *cam = ren->GetActiveCamera();
|
|
|
|
// Keep track of the projection matrix - we'll need it in a couple of places
|
|
// Get the projection matrix. The method is called perspective, but
|
|
// the matrix is valid for perspective and parallel viewing transforms.
|
|
// Don't replace this with the GetCompositePerspectiveTransformMatrix because that
|
|
|
|
// turns off stereo rendering!!!
|
|
this->PerspectiveTransform->SetMatrix(cam->GetPerspectiveTransformMatrix(aspect[0]/aspect[1],
|
|
0.0, 1.0 ));
|
|
this->PerspectiveTransform->Concatenate(cam->GetViewTransformMatrix());
|
|
this->PerspectiveMatrix->DeepCopy(this->PerspectiveTransform->GetMatrix());
|
|
|
|
|
|
// Compute some matrices from voxels to view and vice versa based
|
|
// on the whole input
|
|
|
|
this->VoxelTransformMatrix->DeepCopy( vol->GetMatrix() );
|
|
|
|
|
|
|
|
this->ComputeMatrices( this->GetInput(), vol );
|
|
this->ParallelProjection = cam->GetParallelProjection();
|
|
|
|
// How big is the viewport in pixels?
|
|
double *viewport = ren->GetViewport();
|
|
|
|
int *renWinSize = ren->GetRenderWindow()->GetSize();
|
|
|
|
// Save this so that we can restore it if the image is cancelled
|
|
int oldImageSampleDistance = this->ImageSampleDistance;
|
|
|
|
// If we are automatically adjusting the size to achieve a desired frame
|
|
// rate, then do that adjustment here. Base the new image sample distance
|
|
// on the previous one and the previous render time. Don't let
|
|
// the adjusted image sample distance be less than the minimum image sample
|
|
// distance or more than the maximum image sample distance.
|
|
if ( this->AutoAdjustSampleDistances )
|
|
{
|
|
float oldTime = this->RetrieveRenderTime( ren, vol );
|
|
float newTime = vol->GetAllocatedRenderTime();
|
|
this->ImageSampleDistance = (this->ImageSampleDistance * int(0.5+(sqrt(oldTime / newTime))));//int(rint(sqrt(oldTime / newTime))));
|
|
this->ImageSampleDistance =
|
|
(this->ImageSampleDistance>this->MaximumImageSampleDistance)?
|
|
(this->MaximumImageSampleDistance):(this->ImageSampleDistance);
|
|
this->ImageSampleDistance =
|
|
(this->ImageSampleDistance<this->MinimumImageSampleDistance)?
|
|
|
|
(this->MinimumImageSampleDistance):(this->ImageSampleDistance);
|
|
}
|
|
|
|
vol->UpdateScalarOpacityforSampleSize( ren, this->ImageSampleDistance );
|
|
|
|
// The full image fills the viewport. First, compute the actual viewport
|
|
// size, then divide by the ImageSampleDistance to find the full image
|
|
// size in pixels
|
|
this->ImageViewportSize[0] = static_cast<int>
|
|
(static_cast<float>(renWinSize[0]) * (viewport[2]-viewport[0]));
|
|
this->ImageViewportSize[1] = static_cast<int>
|
|
(static_cast<float>(renWinSize[1]) * (viewport[3]-viewport[1]));
|
|
|
|
this->ImageViewportSize[0] =
|
|
static_cast<int>( static_cast<float>(this->ImageViewportSize[0]) );
|
|
this->ImageViewportSize[1] =
|
|
static_cast<int>( static_cast<float>(this->ImageViewportSize[1]) );
|
|
|
|
/*
|
|
this->ImageViewportSize[0] =
|
|
static_cast<int>( static_cast<float>(this->ImageViewportSize[0]) /
|
|
static_cast<float>(this->ImageSampleDistance) );
|
|
this->ImageViewportSize[1] =
|
|
static_cast<int>( static_cast<float>(this->ImageViewportSize[1]) /
|
|
static_cast<float>(this->ImageSampleDistance) );
|
|
*/
|
|
vol->UpdateTransferFunctions( ren );
|
|
|
|
vol->UpdateScalarOpacityforSampleSize( ren, this->ImageSampleDistance );
|
|
|
|
this->Shade = vol->GetProperty()->GetShade();
|
|
|
|
this->GradientEstimator->SetInput( this->GetInput() );
|
|
|
|
if ( this->Shade )
|
|
{
|
|
this->GradientShader->UpdateShadingTable( ren, vol,
|
|
this->GradientEstimator );
|
|
this->EncodedNormals =
|
|
this->GradientEstimator->GetEncodedNormals();
|
|
|
|
this->RedDiffuseShadingTable =
|
|
this->GradientShader->GetRedDiffuseShadingTable(vol);
|
|
this->GreenDiffuseShadingTable =
|
|
this->GradientShader->GetGreenDiffuseShadingTable(vol);
|
|
this->BlueDiffuseShadingTable =
|
|
this->GradientShader->GetBlueDiffuseShadingTable(vol);
|
|
|
|
this->RedSpecularShadingTable =
|
|
this->GradientShader->GetRedSpecularShadingTable(vol);
|
|
this->GreenSpecularShadingTable =
|
|
this->GradientShader->GetGreenSpecularShadingTable(vol);
|
|
this->BlueSpecularShadingTable =
|
|
this->GradientShader->GetBlueSpecularShadingTable(vol);
|
|
}
|
|
else
|
|
{
|
|
|
|
this->EncodedNormals = NULL;
|
|
this->RedDiffuseShadingTable = NULL;
|
|
this->GreenDiffuseShadingTable = NULL;
|
|
this->BlueDiffuseShadingTable = NULL;
|
|
this->RedSpecularShadingTable = NULL;
|
|
this->GreenSpecularShadingTable = NULL;
|
|
this->BlueSpecularShadingTable = NULL;
|
|
}
|
|
|
|
// If we have non-constant opacity on the gradient magnitudes,
|
|
// we need to use the gradient magnitudes to look up the opacity
|
|
if ( vol->GetGradientOpacityConstant() == -1.0 )
|
|
{
|
|
this->GradientMagnitudes =
|
|
this->GradientEstimator->GetGradientMagnitudes();
|
|
}
|
|
else
|
|
{
|
|
this->GradientMagnitudes = NULL;
|
|
}
|
|
|
|
float bounds[6];
|
|
int dim[3];
|
|
|
|
|
|
this->GetInput()->GetDimensions(dim);
|
|
bounds[0] = bounds[2] = bounds[4] = 0.0;
|
|
bounds[1] = static_cast<float>(dim[0]-1);
|
|
bounds[3] = static_cast<float>(dim[1]-1);
|
|
bounds[5] = static_cast<float>(dim[2]-1);
|
|
|
|
float voxelPoint[3];
|
|
float viewPoint[8][4];
|
|
int i, j, k;
|
|
// unsigned char *ucptr;
|
|
float minX, minY, maxX, maxY, minZ, maxZ;
|
|
|
|
minX = 1.0;
|
|
minY = 1.0;
|
|
maxX = -1.0;
|
|
maxY = -1.0;
|
|
minZ = 1.0;
|
|
maxZ = 0.0;
|
|
|
|
double camPos[3];
|
|
double worldBounds[6];
|
|
vol->GetBounds( worldBounds );
|
|
int insideFlag = 0;
|
|
ren->GetActiveCamera()->GetPosition( camPos );
|
|
if ( camPos[0] >= worldBounds[0] &&
|
|
camPos[0] <= worldBounds[1] &&
|
|
camPos[1] >= worldBounds[2] &&
|
|
camPos[1] <= worldBounds[3] &&
|
|
camPos[2] >= worldBounds[4] &&
|
|
camPos[2] <= worldBounds[5] )
|
|
{
|
|
insideFlag = 1;
|
|
}
|
|
|
|
|
|
|
|
// Copy the voxelsToView matrix to 16 floats
|
|
float voxelsToViewMatrix[16];
|
|
for ( j = 0; j < 4; j++ )
|
|
{
|
|
for ( i = 0; i < 4; i++ )
|
|
{
|
|
voxelsToViewMatrix[j*4+i] =
|
|
static_cast<float>(this->VoxelsToViewMatrix->GetElement(j,i));
|
|
}
|
|
}
|
|
|
|
// Convert the voxel bounds to view coordinates to find out the
|
|
// size and location of the image we need to generate.
|
|
int idx = 0;
|
|
if ( insideFlag )
|
|
{
|
|
minX = -1.0;
|
|
maxX = 1.0;
|
|
minY = -1.0;
|
|
maxY = 1.0;
|
|
minZ = 0.001;
|
|
maxZ = 0.001;
|
|
}
|
|
else
|
|
|
|
{
|
|
for ( k = 0; k < 2; k++ )
|
|
{
|
|
voxelPoint[2] = bounds[4+k];
|
|
for ( j = 0; j < 2; j++ )
|
|
{
|
|
voxelPoint[1] = bounds[2+j];
|
|
for ( i = 0; i < 2; i++ )
|
|
{
|
|
voxelPoint[0] = bounds[i];
|
|
vtkVSWMultiplyPointMacro( voxelPoint, viewPoint[idx],
|
|
voxelsToViewMatrix );
|
|
|
|
minX = (viewPoint[idx][0]<minX)?(viewPoint[idx][0]):(minX);
|
|
minY = (viewPoint[idx][1]<minY)?(viewPoint[idx][1]):(minY);
|
|
maxX = (viewPoint[idx][0]>maxX)?(viewPoint[idx][0]):(maxX);
|
|
maxY = (viewPoint[idx][1]>maxY)?(viewPoint[idx][1]):(maxY);
|
|
minZ = (viewPoint[idx][2]<minZ)?(viewPoint[idx][2]):(minZ);
|
|
maxZ = (viewPoint[idx][2]>maxZ)?(viewPoint[idx][2]):(maxZ);
|
|
idx++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if ( minZ < 0.001 || maxZ > 0.9999 )
|
|
{
|
|
minX = -1.0;
|
|
maxX = 1.0;
|
|
minY = -1.0;
|
|
maxY = 1.0;
|
|
insideFlag = 1;
|
|
}
|
|
|
|
this->MinimumViewDistance =
|
|
(minZ<0.001)?(0.001):((minZ>0.999)?(0.999):(minZ));
|
|
|
|
if ( !ren->GetRenderWindow()->GetAbortRender() )
|
|
this->FactorViewMatrix();
|
|
|
|
if ( this->ClippingPlanes )
|
|
this->InitializeClippingPlanes( this->ClippingPlanes );
|
|
|
|
if (this->IntermixIntersectingGeometry == 1)
|
|
{
|
|
if ( !ren->GetRenderWindow()->GetAbortRender() )
|
|
this->ExtractZBuffer(ren,vol);
|
|
}
|
|
|
|
|
|
if ( !ren->GetRenderWindow()->GetAbortRender() )
|
|
this->CompositeIntermediate();
|
|
|
|
if ( !ren->GetRenderWindow()->GetAbortRender() )
|
|
this->RenderTexture(ren,vol);
|
|
|
|
if ( !ren->GetRenderWindow()->GetAbortRender() )
|
|
{
|
|
this->Timer->StopTimer();
|
|
this->TimeToDraw = this->Timer->GetElapsedTime();
|
|
this->StoreRenderTime( ren, vol, this->TimeToDraw );
|
|
}
|
|
// Restore the image sample distance so that automatic adjustment
|
|
// will work correctly
|
|
else
|
|
{
|
|
this->ImageSampleDistance = oldImageSampleDistance;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
#if 0
|
|
vtkVolumeShearWarpMapper *vtkVolumeShearWarpMapper::New()
|
|
{
|
|
|
|
// First try to create the object from the vtkObjectFactory
|
|
vtkObject* ret =
|
|
vtkGraphicsFactory::CreateInstance("vtkVolumeShearWarpMapper");
|
|
|
|
if (ret != NULL)
|
|
return (vtkVolumeShearWarpMapper*)ret;
|
|
|
|
return vtkOpenGLVolumeShearWarpMapper::New();
|
|
}
|
|
#endif
|
|
|
|
void vtkVolumeShearWarpMapper::ComputeMatrices( vtkImageData *data,
|
|
vtkVolume *vol )
|
|
{
|
|
// Get the data spacing. This scaling is not accounted for in
|
|
// the volume's matrix, so we must add it in.
|
|
double volumeSpacing[3];
|
|
data->GetSpacing( volumeSpacing );
|
|
|
|
// Get the origin of the data. This translation is not accounted for in
|
|
// the volume's matrix, so we must add it in.
|
|
double volumeOrigin[3];
|
|
data->GetOrigin( volumeOrigin );
|
|
|
|
// Get the dimensions of the data.
|
|
int volumeDimensions[3];
|
|
data->GetDimensions( volumeDimensions );
|
|
|
|
// Create some transform objects that we will need later
|
|
vtkTransform *voxelsTransform = vtkTransform::New();
|
|
vtkTransform *voxelsToViewTransform = vtkTransform::New();
|
|
|
|
// Get the volume matrix. This is a volume to world matrix right now.
|
|
// We'll need to invert it, translate by the origin and scale by the
|
|
// spacing to change it to a world to voxels matrix.
|
|
vtkMatrix4x4 *volMatrix = vtkMatrix4x4::New();
|
|
volMatrix->DeepCopy( vol->GetMatrix() );
|
|
|
|
voxelsToViewTransform->SetMatrix( volMatrix );
|
|
|
|
// Create a transform that will account for the scaling and translation of
|
|
// the scalar data. The is the volume to voxels matrix.
|
|
voxelsTransform->Identity();
|
|
voxelsTransform->Translate(volumeOrigin[0],
|
|
volumeOrigin[1],
|
|
volumeOrigin[2] );
|
|
|
|
voxelsTransform->Scale( volumeSpacing[0],
|
|
|
|
volumeSpacing[1],
|
|
volumeSpacing[2] );
|
|
|
|
|
|
|
|
// Now concatenate the volume's matrix with this scalar data matrix
|
|
voxelsToViewTransform->PreMultiply();
|
|
voxelsToViewTransform->Concatenate( voxelsTransform->GetMatrix() );
|
|
|
|
// Now we actually have the world to voxels matrix - copy it out
|
|
this->WorldToVoxelsMatrix->DeepCopy( voxelsToViewTransform->GetMatrix() );
|
|
this->WorldToVoxelsMatrix->Invert();
|
|
|
|
// We also want to invert this to get voxels to world
|
|
this->VoxelsToWorldMatrix->DeepCopy( voxelsToViewTransform->GetMatrix() );
|
|
|
|
|
|
// Compute the voxels to view transform by concatenating the
|
|
// voxels to world matrix with the projection matrix (world to view)
|
|
|
|
voxelsToViewTransform->PostMultiply();
|
|
voxelsToViewTransform->Concatenate( this->PerspectiveMatrix );
|
|
|
|
|
|
|
|
this->VoxelsToViewMatrix->DeepCopy( voxelsToViewTransform->GetMatrix() );
|
|
|
|
this->ViewToVoxelsMatrix->DeepCopy( this->VoxelsToViewMatrix );
|
|
this->ViewToVoxelsMatrix->Invert();
|
|
|
|
voxelsToViewTransform->Delete();
|
|
voxelsTransform->Delete();
|
|
volMatrix->Delete();
|
|
}
|
|
|
|
float vtkVolumeShearWarpMapper::RetrieveRenderTime( vtkRenderer *ren,
|
|
vtkVolume *vol )
|
|
{
|
|
int i;
|
|
|
|
for ( i = 0; i < this->RenderTableEntries; i++ )
|
|
{
|
|
if ( this->RenderVolumeTable[i] == vol &&
|
|
this->RenderRendererTable[i] == ren )
|
|
{
|
|
return this->RenderTimeTable[i];
|
|
}
|
|
}
|
|
|
|
return 0.0;
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::StoreRenderTime( vtkRenderer *ren,
|
|
vtkVolume *vol,
|
|
float time )
|
|
{
|
|
int i;
|
|
for ( i = 0; i < this->RenderTableEntries; i++ )
|
|
|
|
{
|
|
if ( this->RenderVolumeTable[i] == vol &&
|
|
this->RenderRendererTable[i] == ren )
|
|
{
|
|
this->RenderTimeTable[i] = time;
|
|
return;
|
|
}
|
|
}
|
|
|
|
|
|
// Need to increase size
|
|
if ( this->RenderTableEntries >= this->RenderTableSize )
|
|
{
|
|
if ( this->RenderTableSize == 0 )
|
|
{
|
|
this->RenderTableSize = 10;
|
|
}
|
|
else
|
|
{
|
|
this->RenderTableSize *= 2;
|
|
}
|
|
|
|
float *oldTimePtr = this->RenderTimeTable;
|
|
vtkVolume **oldVolumePtr = this->RenderVolumeTable;
|
|
|
|
vtkRenderer **oldRendererPtr = this->RenderRendererTable;
|
|
|
|
this->RenderTimeTable = new float [this->RenderTableSize];
|
|
this->RenderVolumeTable = new vtkVolume *[this->RenderTableSize];
|
|
this->RenderRendererTable = new vtkRenderer *[this->RenderTableSize];
|
|
|
|
for (i = 0; i < this->RenderTableEntries; i++ )
|
|
{
|
|
this->RenderTimeTable[i] = oldTimePtr[i];
|
|
this->RenderVolumeTable[i] = oldVolumePtr[i];
|
|
this->RenderRendererTable[i] = oldRendererPtr[i];
|
|
}
|
|
|
|
delete [] oldTimePtr;
|
|
delete [] oldVolumePtr;
|
|
delete [] oldRendererPtr;
|
|
}
|
|
|
|
this->RenderTimeTable[this->RenderTableEntries] = time;
|
|
this->RenderVolumeTable[this->RenderTableEntries] = vol;
|
|
this->RenderRendererTable[this->RenderTableEntries] = ren;
|
|
|
|
this->RenderTableEntries++;
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::CompositeIntermediate()
|
|
{
|
|
int scalarType = this->GetInput()->GetScalarType();
|
|
int interpolationType = this->Volume->GetProperty()->GetInterpolationType();
|
|
|
|
|
|
this->ImageWidth = 32;
|
|
this->ImageHeight = 32;
|
|
|
|
while (this->ImageWidth < this->MaximumIntermediateDimension)
|
|
this->ImageWidth = this->ImageWidth << 1;
|
|
|
|
while (this->ImageHeight < this->MaximumIntermediateDimension)
|
|
this->ImageHeight = this->ImageHeight << 1;
|
|
|
|
int imageSize = this->ImageWidth * this->ImageHeight;
|
|
|
|
if (imageSize > this->AllocatedSize)
|
|
{
|
|
this->AllocatedSize = imageSize;
|
|
|
|
if (this->ImageData != NULL)
|
|
delete[] this->ImageData;
|
|
|
|
this->ImageData = new unsigned char[imageSize*4];
|
|
|
|
if (this->IntemediateImage != NULL)
|
|
delete this->IntemediateImage;
|
|
|
|
this->IntemediateImage = new vtkShearWarpRLEImage(imageSize);
|
|
}
|
|
else
|
|
this->IntemediateImage->Clear();
|
|
|
|
if (this->RunlengthEncoding == 1)
|
|
{
|
|
if (this->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
if (this->EncodedVolume != NULL)
|
|
{
|
|
if (this->EncodedVolume->IsScalarEncoded())
|
|
{
|
|
if (this->EncodedVolume->GetIsoValue() != this->IsoValue)
|
|
{
|
|
delete this->EncodedVolume;
|
|
this->EncodedVolume = NULL;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
delete this->EncodedVolume;
|
|
this->EncodedVolume = NULL;
|
|
}
|
|
}
|
|
|
|
|
|
if (this->EncodedVolume == NULL)
|
|
{
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
{
|
|
vtkShearWarpRLEVolume<unsigned char> *encodedVolume = new vtkShearWarpRLEVolume<unsigned char>();
|
|
encodedVolume->encodeScalar(this->GetInput(),this->Volume,this->GradientEstimator,this->IsoValue);
|
|
this->EncodedVolume = encodedVolume;
|
|
break;
|
|
}
|
|
case VTK_UNSIGNED_SHORT:
|
|
{
|
|
vtkShearWarpRLEVolume<unsigned short> *encodedVolume = new vtkShearWarpRLEVolume<unsigned short>();
|
|
encodedVolume->encodeScalar(this->GetInput(),this->Volume,this->GradientEstimator,this->IsoValue);
|
|
this->EncodedVolume = encodedVolume;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
if ( interpolationType == VTK_NEAREST_INTERPOLATION )
|
|
CompositeIntermediateNearestRLE<unsigned char>(this->IntemediateImage,this);
|
|
else
|
|
CompositeIntermediateLinearRLE<unsigned char>(this->IntemediateImage,this);
|
|
|
|
break;
|
|
|
|
case VTK_UNSIGNED_SHORT:
|
|
if ( interpolationType == VTK_NEAREST_INTERPOLATION )
|
|
|
|
CompositeIntermediateNearestRLE<unsigned short>(this->IntemediateImage,this);
|
|
else
|
|
// if(ParallelProjection)
|
|
CompositeIntermediateLinearRLE<unsigned short>(this->IntemediateImage,this);
|
|
// else
|
|
// CompositeIntermediateLinearRLEPerspective<unsigned short>(this->IntemediateImage,this);
|
|
|
|
break;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// If the scalar opacity transfer function has been modified
|
|
// the runlength encoding has to be redone
|
|
unsigned long scalarOpacityMTime = this->Volume->GetProperty()->GetScalarOpacity()->GetMTime();
|
|
|
|
if (this->EncodedVolume != NULL)
|
|
{
|
|
if ( scalarOpacityMTime > this->ScalarOpacityMTime)
|
|
{
|
|
delete this->EncodedVolume;
|
|
this->EncodedVolume = NULL;
|
|
}
|
|
}
|
|
|
|
if (this->EncodedVolume == NULL)
|
|
{
|
|
this->ScalarOpacityMTime = scalarOpacityMTime;
|
|
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
{
|
|
vtkShearWarpRLEVolume<unsigned char> *encodedVolume = new vtkShearWarpRLEVolume<unsigned char>();
|
|
encodedVolume->encodeOpacity(this->GetInput(),this->Volume,this->GradientEstimator,0.0f);
|
|
this->EncodedVolume = encodedVolume;
|
|
break;
|
|
}
|
|
|
|
|
|
case VTK_UNSIGNED_SHORT:
|
|
{
|
|
vtkShearWarpRLEVolume<unsigned short> *encodedVolume = new vtkShearWarpRLEVolume<unsigned short>();
|
|
|
|
encodedVolume->encodeOpacity(this->GetInput(),this->Volume,this->GradientEstimator,0.0f);
|
|
this->EncodedVolume = encodedVolume;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
if ( interpolationType == VTK_NEAREST_INTERPOLATION )
|
|
CompositeIntermediateNearestRLE<unsigned char>(this->IntemediateImage,this);
|
|
else
|
|
CompositeIntermediateLinearRLE<unsigned char>(this->IntemediateImage,this);
|
|
|
|
break;
|
|
|
|
case VTK_UNSIGNED_SHORT:
|
|
if ( interpolationType == VTK_NEAREST_INTERPOLATION )
|
|
CompositeIntermediateNearestRLE<unsigned short>(this->IntemediateImage,this);
|
|
else
|
|
// if(ParallelProjection)
|
|
CompositeIntermediateLinearRLE<unsigned short>(this->IntemediateImage,this);
|
|
// else
|
|
// CompositeIntermediateLinearRLEPerspective<unsigned short>(this->IntemediateImage,this);
|
|
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
else if (this->FastClassification == 1)
|
|
{
|
|
if (this->EncodedVolume != NULL)
|
|
{
|
|
delete this->EncodedVolume;
|
|
this->EncodedVolume = NULL;
|
|
}
|
|
|
|
if (this->Octree == NULL)
|
|
{
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
this->Octree = new vtkShearWarpOctree<unsigned char>();
|
|
dynamic_cast< vtkShearWarpOctree<unsigned char> *>(this->Octree)->build(this->GetInput());
|
|
break;
|
|
case VTK_UNSIGNED_SHORT:
|
|
this->Octree = new vtkShearWarpOctree<unsigned short>();
|
|
dynamic_cast< vtkShearWarpOctree<unsigned short> *>(this->Octree)->build(this->GetInput());
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (this->FunctionType == VTK_SHEAR_WARP_ISOSURFACE_FUNCTION)
|
|
{
|
|
if (!this->Octree->IsScalarEncoded() || this->Octree->GetIsoValue() != this->IsoValue)
|
|
{
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
dynamic_cast< vtkShearWarpOctree<unsigned char> *>(this->Octree)->classifyScalar((unsigned char)this->IsoValue);
|
|
break;
|
|
case VTK_UNSIGNED_SHORT:
|
|
dynamic_cast< vtkShearWarpOctree<unsigned short> *>(this->Octree)->classifyScalar((unsigned short)this->IsoValue);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
unsigned long scalarOpacityMTime = this->Volume->GetProperty()->GetScalarOpacity()->GetMTime();
|
|
|
|
if ( scalarOpacityMTime > this->ScalarOpacityMTime)
|
|
{
|
|
this->ScalarOpacityMTime = scalarOpacityMTime;
|
|
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
dynamic_cast< vtkShearWarpOctree<unsigned char> *>(this->Octree)->classifyOpacity(this->Volume);
|
|
break;
|
|
case VTK_UNSIGNED_SHORT:
|
|
dynamic_cast< vtkShearWarpOctree<unsigned short> *>(this->Octree)->classifyOpacity(this->Volume);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if ( interpolationType == VTK_NEAREST_INTERPOLATION )
|
|
{
|
|
// Nearest neighbor
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
CompositeIntermediateNearestUnclassified<unsigned char>(this->IntemediateImage,this);
|
|
break;
|
|
case VTK_UNSIGNED_SHORT:
|
|
CompositeIntermediateNearestUnclassified<unsigned short>(this->IntemediateImage,this);
|
|
break;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Linear interpolation
|
|
switch ( scalarType )
|
|
{
|
|
|
|
case VTK_UNSIGNED_CHAR:
|
|
CompositeIntermediateLinearUnclassified<unsigned char>(this->IntemediateImage,this);
|
|
break;
|
|
case VTK_UNSIGNED_SHORT:
|
|
CompositeIntermediateLinearUnclassified<unsigned short>(this->IntemediateImage,this);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if ( interpolationType == VTK_NEAREST_INTERPOLATION )
|
|
{
|
|
// Nearest neighbor
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
CompositeIntermediateNearestSimple<unsigned char>(this->IntemediateImage,this);
|
|
break;
|
|
case VTK_UNSIGNED_SHORT:
|
|
CompositeIntermediateNearestSimple<unsigned short>(this->IntemediateImage,this);
|
|
break;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Linear interpolation
|
|
switch ( scalarType )
|
|
{
|
|
case VTK_UNSIGNED_CHAR:
|
|
CompositeIntermediateLinearSimple<unsigned char>(this->IntemediateImage,this);
|
|
break;
|
|
case VTK_UNSIGNED_SHORT:
|
|
|
|
CompositeIntermediateLinearSimple<unsigned short>(this->IntemediateImage,this);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
BuildImage(this->ImageData,this->IntemediateImage->GetPixelData());
|
|
|
|
bool Debug_IntermediateImage = true;
|
|
if(Debug_IntermediateImage)
|
|
{
|
|
FILE *fp;
|
|
fp = fopen("c:\\temp\\output.raw", "wb+");
|
|
fwrite(this->ImageData, 4, this->ImageWidth * this->ImageHeight, fp);
|
|
fclose(fp);
|
|
}
|
|
|
|
}
|
|
|
|
|
|
void vtkVolumeShearWarpMapper::BuildImage(unsigned char *id, vtkShearWarpPixelData *im)
|
|
{
|
|
int i,j;
|
|
float red,green,blue,opacity;
|
|
|
|
int startX = 0;//int(float(this->ImageWidth) * 0.5f - float(this->IntermediateWidth) * 0.5f);
|
|
int startY = 0;//int(float(this->ImageHeight) * 0.5f - float(this->IntermediateHeight) * 0.5f);
|
|
|
|
for (j = 0; j < this->ImageHeight; j++)
|
|
|
|
{
|
|
for (i = 0; i < this->ImageWidth; i++)
|
|
{
|
|
red = 0.0f;
|
|
|
|
green = 0.0f;
|
|
blue = 0.0f;
|
|
opacity = 0.0f;
|
|
|
|
if (/*i >= startX && j >= startY && */i < startX+this->IntermediateWidth && j < startY+this->IntermediateHeight)
|
|
{
|
|
red = im->Red;
|
|
green = im->Green;
|
|
blue = im->Blue;
|
|
opacity = im->Opacity;
|
|
|
|
if (red > 1.0f)
|
|
red = 1.0f;
|
|
else if (red < 0.0f)
|
|
red = 0.0f;
|
|
|
|
if (green > 1.0f)
|
|
|
|
green = 1.0f;
|
|
else if (green < 0.0f)
|
|
green = 0.0f;
|
|
|
|
if (blue > 1.0f)
|
|
blue = 1.0f;
|
|
else if (blue < 0.0f)
|
|
blue = 0.0f;
|
|
|
|
if (opacity > 1.0f)
|
|
opacity = 1.0f;
|
|
else if (opacity < 0.0f)
|
|
opacity = 0.0f;
|
|
|
|
im++;
|
|
}
|
|
/*
|
|
else
|
|
{
|
|
red = 1.0f;
|
|
green = 0.0f;
|
|
blue = 0.0f;
|
|
opacity = 1.0f;
|
|
}
|
|
*/
|
|
id[0] = (unsigned char) (255.0f * red);
|
|
id[1] = (unsigned char) (255.0f * green);
|
|
id[2] = (unsigned char) (255.0f * blue);
|
|
id[3] = (unsigned char) (255.0f * opacity);
|
|
|
|
id += 4;
|
|
}
|
|
}
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::ExtractZBuffer(vtkRenderer *ren, vtkVolume *vol)
|
|
{
|
|
float iposition[4][2];
|
|
float itranslation[2];
|
|
float isx, isy;
|
|
float ix,iy;
|
|
float w00,w01,w10,w11;
|
|
|
|
int *renWinSize = ren->GetRenderWindow()->GetSize();
|
|
|
|
// float *viewport = ren->GetViewport();
|
|
|
|
// The coefficients of the 2D warp matrix
|
|
w00 = WarpMatrix->Element[0][0];
|
|
w01 = WarpMatrix->Element[0][1];
|
|
w10 = WarpMatrix->Element[1][0];
|
|
w11 = WarpMatrix->Element[1][1];
|
|
|
|
ix = (float) this->vtkVolumeShearWarpMapper::IntermediateWidth * this->ImageSampleDistance;
|
|
iy = (float) this->vtkVolumeShearWarpMapper::IntermediateHeight * this->ImageSampleDistance;
|
|
|
|
|
|
// Intermediate
|
|
iposition[0][0] = 0.0f * w00 + 0.0f * w01;
|
|
iposition[0][1] = 0.0f * w10 + 0.0f * w11;
|
|
|
|
iposition[1][0] = ix * w00 + 0.0f * w01;
|
|
iposition[1][1] = ix * w10 + 0.0f * w11;
|
|
|
|
iposition[2][0] = ix * w00 + iy * w01;
|
|
iposition[2][1] = ix * w10 + iy * w11;
|
|
|
|
iposition[3][0] = 0.0f * w00 + iy * w01;
|
|
iposition[3][1] = 0.0f * w10 + iy * w11;
|
|
|
|
// Intermediate
|
|
itranslation[0] = ix * 0.5f * w00 + iy * 0.5f * w01;
|
|
itranslation[1] = ix * 0.5f * w10 + iy * 0.5f * w11;
|
|
|
|
// Intermediate
|
|
isx = 1;//(float) this->vtkVolumeShearWarpMapper::ImageSampleDistance;
|
|
isy = 1;//(float) this->vtkVolumeShearWarpMapper::ImageSampleDistance;
|
|
|
|
double *t = vol->GetCenter();
|
|
float a[4] = {t[0],t[1],t[2],1.0f};
|
|
float b[4];
|
|
|
|
|
|
this->PerspectiveMatrix->MultiplyPoint(a,b);
|
|
b[2] = 0.0f;
|
|
|
|
float x1 = isx*(iposition[0][0]-itranslation[0]) + b[0] * (float) renWinSize[0] * 0.5f;
|
|
float x2 = isx*(iposition[1][0]-itranslation[0]) + b[0] * (float) renWinSize[0] * 0.5f;
|
|
float x3 = isx*(iposition[2][0]-itranslation[0]) + b[0] * (float) renWinSize[0] * 0.5f;
|
|
float x4 = isx*(iposition[3][0]-itranslation[0]) + b[0] * (float) renWinSize[0] * 0.5f;
|
|
|
|
float y1 = isy*(iposition[0][1]-itranslation[1]) + b[1] * (float) renWinSize[1] * 0.5f;
|
|
float y2 = isy*(iposition[1][1]-itranslation[1]) + b[1] * (float) renWinSize[1] * 0.5f;
|
|
float y3 = isy*(iposition[2][1]-itranslation[1]) + b[1] * (float) renWinSize[1] * 0.5f;
|
|
float y4 = isy*(iposition[3][1]-itranslation[1]) + b[1] * (float) renWinSize[1] * 0.5f;
|
|
|
|
|
|
|
|
x1 = x1 + renWinSize[0] * 0.5f;
|
|
y1 = y1 + renWinSize[1] * 0.5f;
|
|
|
|
|
|
x2 = x2 + renWinSize[0] * 0.5f;
|
|
y2 = y2 + renWinSize[1] * 0.5f;
|
|
|
|
x3 = x3 + renWinSize[0] * 0.5f;
|
|
y3 = y3 + renWinSize[1] * 0.5f;
|
|
|
|
|
|
x4 = x4 + renWinSize[0] * 0.5f;
|
|
y4 = y4 + renWinSize[1] * 0.5f;
|
|
|
|
float minx = x1;
|
|
float miny = y1;
|
|
float maxx = x1;
|
|
float maxy = y1;
|
|
|
|
if (x2 < minx)
|
|
minx = x2;
|
|
|
|
if (y2 < miny)
|
|
miny = y2;
|
|
|
|
if (x3 < minx)
|
|
minx = x3;
|
|
|
|
if (y3 < miny)
|
|
miny = y3;
|
|
|
|
if (x4 < minx)
|
|
minx = x4;
|
|
|
|
if (y4 < miny)
|
|
miny = y4;
|
|
|
|
|
|
if (x2 > maxx)
|
|
|
|
maxx = x2;
|
|
|
|
if (y2 > maxy)
|
|
maxy = y2;
|
|
|
|
if (x3 > maxx)
|
|
maxx = x3;
|
|
|
|
if (y3 > maxy)
|
|
maxy = y3;
|
|
|
|
if (x4 > maxx)
|
|
maxx = x4;
|
|
|
|
if (y4 > maxy)
|
|
maxy = y4;
|
|
|
|
int left = 0;
|
|
int top = 0;
|
|
|
|
if (minx < 0)
|
|
{
|
|
left = (int) -minx;
|
|
minx = 0;
|
|
}
|
|
|
|
if (miny < 0)
|
|
{
|
|
top = (int) -miny;
|
|
miny = 0;
|
|
}
|
|
|
|
if (maxx > renWinSize[0] - 1)
|
|
maxx = renWinSize[0] - 1;
|
|
|
|
if (maxy > renWinSize[1] - 1)
|
|
maxy = renWinSize[1] - 1;
|
|
|
|
int zx1 = int(minx + 0.5f);
|
|
int zy1 = int(miny + 0.5f);
|
|
int zx2 = int(maxx - 0.5f);
|
|
int zy2 = int(maxy - 0.5f);
|
|
|
|
|
|
this->ZBuffer = ren->GetRenderWindow()->GetZbufferData(zx1,zy1,zx2,zy2);
|
|
this->ZBufferSize[0] = zx2 - zx1 + 1;
|
|
this->ZBufferSize[1] = zy2 - zy1 + 1;
|
|
|
|
this->IntermediateZBuffer = new float[(this->ImageSampleDistance*this->IntermediateWidth) * (this->ImageSampleDistance*this->IntermediateHeight)];
|
|
|
|
// this->Unwarp(this->IntermediateZBuffer,this->IntermediateWidth,this->IntermediateHeight,this->ZBuffer,this->ZBufferSize[0],this->ZBufferSize[1],this->WarpMatrix);
|
|
this->Unwarp(this->IntermediateZBuffer, this->ImageSampleDistance*this->IntermediateWidth, this->ImageSampleDistance*this->IntermediateHeight, this->ZBuffer, left, top, this->ZBufferSize[0],this->ZBufferSize[1],this->WarpMatrix);
|
|
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::Unwarp(float *destination, int dWidth, int dHeight, float *source, int left, int top, int sWidth, int sHeight, vtkMatrix4x4* w)
|
|
|
|
{
|
|
float xs, ys; // source image coordinates
|
|
float xsMin,ysMin;
|
|
int i, j; // counters
|
|
float xd, yd; // precomputed destination values
|
|
float pc; // perspective correction
|
|
float inv00, inv01, inv03; // elements of 1st row of inverted warp matrix
|
|
float inv10, inv11, inv13; // elements of 2nd row of inverted warp matrix
|
|
float inv30, inv31, inv33; // elements of 3rd row of inverted warp matrix
|
|
|
|
inv00 = w->Element[0][0];
|
|
inv01 = w->Element[0][1];
|
|
inv03 = w->Element[0][3];
|
|
inv10 = w->Element[1][0];
|
|
inv11 = w->Element[1][1];
|
|
inv13 = w->Element[1][3];
|
|
inv30 = w->Element[3][0];
|
|
inv31 = w->Element[3][1];
|
|
inv33 = w->Element[3][3];
|
|
|
|
xsMin = 4096;
|
|
ysMin = 4096;
|
|
|
|
for (j=0; j < dHeight; j++)
|
|
{
|
|
yd = (float)j;
|
|
for (i=0; i < dWidth; i++)
|
|
{
|
|
xd = (float)i;
|
|
pc = xd * inv30 + yd * inv31 + inv33;
|
|
xs = ((xd * inv00 + yd * inv01 + inv03)) / pc;
|
|
ys = ((xd * inv10 + yd * inv11 + inv13)) / pc;
|
|
|
|
if (xs < xsMin)
|
|
xsMin = xs;
|
|
|
|
|
|
if (ys < ysMin)
|
|
ysMin = ys;
|
|
}
|
|
}
|
|
|
|
for (j=0; j < dHeight; j++)
|
|
{
|
|
yd = (float)j;
|
|
for (i=0; i < dWidth; i++)
|
|
{
|
|
xd = (float)i;
|
|
pc = xd * inv30 + yd * inv31 + inv33;
|
|
|
|
xs = ((xd * inv00 + yd * inv01 + inv03)) / pc;
|
|
ys = ((xd * inv10 + yd * inv11 + inv13)) / pc;
|
|
|
|
xs -= xsMin;
|
|
ys -= ysMin;
|
|
|
|
|
|
xs -= left;
|
|
ys -= top;
|
|
|
|
float depth;
|
|
|
|
// Check if pixel is inside image
|
|
if (xs > sWidth - 1 || ys > sHeight - 1 || xs < 0 || ys < 0)
|
|
destination[i + j * dWidth] = 0.0f;
|
|
else
|
|
{
|
|
depth = source[int(xs) + int(ys) * sWidth];
|
|
// memcpy(destination + (i + j * dWidth), source + (int(xs) + int(ys) * sWidth), sizeof(float));
|
|
// memcpy(destination + (i + j * dWidth), source + (int(xs) + int(ys) * sWidth), sizeof(float));
|
|
depth = depth * this->PermutedViewToVoxelsMatrix->Element[2][2] + this->PermutedViewToVoxelsMatrix->Element[2][3];
|
|
destination[i + j * dWidth] = depth;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void vtkVolumeShearWarpMapper::InitializeClippingPlanes( vtkPlaneCollection *planes )
|
|
{
|
|
vtkPlane *onePlane;
|
|
double planePoint[4];
|
|
double normalPoint[4];
|
|
float *clippingPlane;
|
|
float d;
|
|
int i;
|
|
|
|
this->ClippingPlaneCount = planes->GetNumberOfItems();
|
|
|
|
if (this->ClippingPlaneCount == 0)
|
|
return;
|
|
|
|
// loop through all the clipping planes
|
|
for ( i = 0; i < this->ClippingPlaneCount; i++ )
|
|
{
|
|
onePlane = (vtkPlane *)planes->GetItemAsObject(i);
|
|
|
|
onePlane->GetOrigin(planePoint);
|
|
onePlane->GetNormal(normalPoint);
|
|
normalPoint[0] += planePoint[0];
|
|
normalPoint[1] += planePoint[1];
|
|
normalPoint[2] += planePoint[2];
|
|
planePoint[3] = 1.0;
|
|
normalPoint[3] = 1.0;
|
|
|
|
this->WorldToVoxelsMatrix->MultiplyPoint(planePoint,planePoint);
|
|
this->WorldToVoxelsMatrix->MultiplyPoint(normalPoint,normalPoint);
|
|
|
|
this->PermutationMatrix->MultiplyPoint(planePoint,planePoint);
|
|
this->PermutationMatrix->MultiplyPoint(normalPoint,normalPoint);
|
|
|
|
clippingPlane = this->ClippingPlane + 4*i;
|
|
clippingPlane[0] = normalPoint[0] - planePoint[0];
|
|
clippingPlane[1] = normalPoint[1] - planePoint[1];
|
|
clippingPlane[2] = normalPoint[2] - planePoint[2];
|
|
|
|
d = sqrt(clippingPlane[0] * clippingPlane[0] + clippingPlane[1] * clippingPlane[1] + clippingPlane[2] * clippingPlane[2]);
|
|
clippingPlane[0] /= d;
|
|
clippingPlane[1] /= d;
|
|
clippingPlane[2] /= d;
|
|
|
|
clippingPlane[3] = clippingPlane[0] * planePoint[0] + clippingPlane[1] * planePoint[1] + clippingPlane[2] * planePoint[2];
|
|
}
|
|
}
|
|
|
|
int vtkVolumeShearWarpMapper::IsVoxelClipped(int x, int y, int z)
|
|
{
|
|
float *clippingPlane;
|
|
|
|
for (int i=0; i<this->ClippingPlaneCount; i++)
|
|
{
|
|
clippingPlane = this->ClippingPlane + 4*i;
|
|
|
|
if (clippingPlane[0] * (float)x +
|
|
clippingPlane[1] * (float)y +
|
|
clippingPlane[2] * (float)z < clippingPlane[3])
|
|
return 1;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// Print the vtkVolumeShearWarpMapper
|
|
void vtkVolumeShearWarpMapper::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
if( this->GradientShader )
|
|
{
|
|
os << indent << "GradientShader:\n";
|
|
this->GradientShader->PrintSelf(os,indent.GetNextIndent());
|
|
}
|
|
else
|
|
{
|
|
os << indent << "GradientShader: (none)\n";
|
|
}
|
|
os << indent << "FunctionType: " << this->FunctionType << "\n";
|
|
os << indent << "MinimumImageSampleDistance: " <<
|
|
this->MinimumImageSampleDistance << "\n";
|
|
if( this->GradientEstimator )
|
|
{
|
|
os << indent << "GradientEstimator:\n";
|
|
this->GradientEstimator->PrintSelf(os,indent.GetNextIndent());
|
|
}
|
|
else
|
|
{
|
|
os << indent << "GradientEstimator: (none)\n";
|
|
}
|
|
os << indent << "AutoAdjustSampleDistances: " <<
|
|
this->AutoAdjustSampleDistances << "\n";
|
|
os << indent << "ParallelProjection: " << this->ParallelProjection <<
|
|
"\n";
|
|
os << indent << "MaximumImageSampleDistance: " <<
|
|
this->MaximumImageSampleDistance << "\n";
|
|
os << indent << "IntermixIntersectingGeometry: " <<
|
|
this->IntermixIntersectingGeometry << "\n";
|
|
os << indent << "IsoValue: " << this->IsoValue << "\n";
|
|
os << indent << "MyPerspectiveProjection: " <<
|
|
this->MyPerspectiveProjection << "\n";
|
|
os << indent << "FastClassification: " << this->FastClassification <<
|
|
"\n";
|
|
os << indent << "RunlengthEncoding: " << this->RunlengthEncoding << "\n";
|
|
os << indent << "ImageSampleDistance: " << this->ImageSampleDistance
|
|
<< "\n";
|
|
}
|
|
|