Cloned library of VTK-5.0.0 with extra build files for internal package management.
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.
 
 
 
 
 
 

1167 lines
36 KiB

/*=========================================================================
Program: Visualization Toolkit
Module: $RCSfile: vtkUnstructuredGridBunykRayCastFunction.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 "vtkUnstructuredGridBunykRayCastFunction.h"
#include "vtkObjectFactory.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridVolumeRayCastMapper.h"
#include "vtkVolume.h"
#include "vtkRenderer.h"
#include "vtkCamera.h"
#include "vtkMatrix4x4.h"
#include "vtkTransform.h"
#include "vtkCell.h"
#include "vtkCellType.h"
#include "vtkMath.h"
#include "vtkPointData.h"
#include "vtkCellArray.h"
#include "vtkDoubleArray.h"
#include "vtkIdList.h"
#include "vtkPiecewiseFunction.h"
#include "vtkColorTransferFunction.h"
#include "vtkVolumeProperty.h"
#include "vtkUnstructuredGridVolumeRayCastIterator.h"
vtkCxxRevisionMacro(vtkUnstructuredGridBunykRayCastFunction, "$Revision: 1.1 $");
vtkStandardNewMacro(vtkUnstructuredGridBunykRayCastFunction);
#define VTK_BUNYKRCF_NUMLISTS 100000
#define VTK_BUNYKRCF_MAX_COMPONENTS 4
template <class T>
vtkIdType TemplateCastRay(
const T *scalars,
vtkUnstructuredGridBunykRayCastFunction *self,
int numComponents,
int x, int y,
double farClipZ,
vtkUnstructuredGridBunykRayCastFunction::Intersection *&intersectionPtr,
vtkUnstructuredGridBunykRayCastFunction::Triangle *&currentTriangle,
vtkIdType &currentTetra,
vtkIdType *intersectedCells,
double *intersectionLengths,
T *nearIntersections,
T *farIntersections,
int maxNumIntersections);
//-----------------------------------------------------------------------------
// This is an internal hidden class.
class vtkUnstructuredGridBunykRayCastIterator
: public vtkUnstructuredGridVolumeRayCastIterator
{
public:
vtkTypeRevisionMacro(vtkUnstructuredGridBunykRayCastIterator,
vtkUnstructuredGridVolumeRayCastIterator);
static vtkUnstructuredGridBunykRayCastIterator *New();
void Initialize(int x, int y);
vtkIdType GetNextIntersections(vtkIdList *intersectedCells,
vtkDoubleArray *intersectionLengths,
vtkDataArray *scalars,
vtkDataArray *nearIntersections,
vtkDataArray *farIntersections);
vtkSetObjectMacro(RayCastFunction, vtkUnstructuredGridBunykRayCastFunction);
vtkGetObjectMacro(RayCastFunction, vtkUnstructuredGridBunykRayCastFunction);
protected:
vtkUnstructuredGridBunykRayCastIterator();
~vtkUnstructuredGridBunykRayCastIterator();
int RayPosition[2];
vtkUnstructuredGridBunykRayCastFunction *RayCastFunction;
vtkUnstructuredGridBunykRayCastFunction::Intersection *IntersectionPtr;
vtkUnstructuredGridBunykRayCastFunction::Triangle *CurrentTriangle;
vtkIdType CurrentTetra;
private:
vtkUnstructuredGridBunykRayCastIterator(const vtkUnstructuredGridBunykRayCastIterator&); // Not implemented
void operator=(const vtkUnstructuredGridBunykRayCastIterator&); // Not implemented
};
vtkCxxRevisionMacro(vtkUnstructuredGridBunykRayCastIterator, "$Revision: 1.1 $");
vtkStandardNewMacro(vtkUnstructuredGridBunykRayCastIterator);
vtkUnstructuredGridBunykRayCastIterator::vtkUnstructuredGridBunykRayCastIterator()
{
this->RayCastFunction = NULL;
}
vtkUnstructuredGridBunykRayCastIterator::~vtkUnstructuredGridBunykRayCastIterator()
{
this->SetRayCastFunction(NULL);
}
void vtkUnstructuredGridBunykRayCastIterator::Initialize(int x, int y)
{
this->RayPosition[0] = x; this->RayPosition[1] = y;
this->IntersectionPtr
= this->RayCastFunction->GetIntersectionList(this->RayPosition[0],
this->RayPosition[1]);
this->CurrentTriangle = NULL;
this->CurrentTetra = -1;
// Intersect cells until we get to Bounds[0] (the near clip plane).
while(TemplateCastRay((const float *)NULL,
this->RayCastFunction, 0,
this->RayPosition[0], this->RayPosition[1],
this->Bounds[0],
this->IntersectionPtr,
this->CurrentTriangle,
this->CurrentTetra,
(vtkIdType *)NULL,
(double *)NULL,
(float *)NULL,
(float *)NULL,
this->MaxNumberOfIntersections) > 0);
}
vtkIdType vtkUnstructuredGridBunykRayCastIterator::GetNextIntersections(
vtkIdList *intersectedCells,
vtkDoubleArray *intersectionLengths,
vtkDataArray *scalars,
vtkDataArray *nearIntersections,
vtkDataArray *farIntersections)
{
if (intersectedCells)
{
intersectedCells->SetNumberOfIds(this->MaxNumberOfIntersections);
}
if (intersectionLengths)
{
intersectionLengths->SetNumberOfComponents(1);
intersectionLengths->SetNumberOfTuples(this->MaxNumberOfIntersections);
}
vtkIdType numIntersections = 0;
if (!scalars)
{
numIntersections = TemplateCastRay
((const float *)NULL, this->RayCastFunction, 0,
this->RayPosition[0], this->RayPosition[1], this->Bounds[1],
this->IntersectionPtr, this->CurrentTriangle, this->CurrentTetra,
(intersectedCells ? intersectedCells->GetPointer(0) : NULL),
(intersectionLengths ? intersectionLengths->GetPointer(0) : NULL),
(float *)NULL, (float *)NULL, this->MaxNumberOfIntersections);
}
else
{
if ( (scalars->GetDataType() != nearIntersections->GetDataType())
|| (scalars->GetDataType() != farIntersections->GetDataType()) )
{
vtkErrorMacro(<< "Data types for scalars do not match up.");
}
nearIntersections->SetNumberOfComponents(scalars->GetNumberOfComponents());
nearIntersections->SetNumberOfTuples(this->MaxNumberOfIntersections);
farIntersections->SetNumberOfComponents(scalars->GetNumberOfComponents());
farIntersections->SetNumberOfTuples(this->MaxNumberOfIntersections);
switch (scalars->GetDataType())
{
vtkTemplateMacro
(numIntersections = TemplateCastRay
((const VTK_TT *)scalars->GetVoidPointer(0),
this->RayCastFunction, scalars->GetNumberOfComponents(),
this->RayPosition[0], this->RayPosition[1], this->Bounds[1],
this->IntersectionPtr, this->CurrentTriangle, this->CurrentTetra,
(intersectedCells ? intersectedCells->GetPointer(0) : NULL),
(intersectionLengths ? intersectionLengths->GetPointer(0) : NULL),
(VTK_TT *)nearIntersections->GetVoidPointer(0),
(VTK_TT *)farIntersections->GetVoidPointer(0),
this->MaxNumberOfIntersections));
}
nearIntersections->SetNumberOfTuples(numIntersections);
farIntersections->SetNumberOfTuples(numIntersections);
}
if (intersectedCells)
{
intersectedCells->SetNumberOfIds(numIntersections);
}
if (intersectionLengths)
{
intersectionLengths->SetNumberOfTuples(numIntersections);
}
return numIntersections;
}
//-----------------------------------------------------------------------------
// Constructor - initially everything to null, and create a matrix for use later
vtkUnstructuredGridBunykRayCastFunction::vtkUnstructuredGridBunykRayCastFunction()
{
this->Renderer = NULL;
this->Volume = NULL;
this->Mapper = NULL;
this->Valid = 0;
this->Points = NULL;
this->Image = NULL;
this->TriangleList = NULL;
this->TetraTriangles = NULL;
this->NumberOfPoints = 0;
this->ImageSize[0] = 0;
this->ImageSize[1] = 0;
this->ViewToWorldMatrix = vtkMatrix4x4::New();
for (int i = 0; i < VTK_BUNYKRCF_MAX_ARRAYS; i++ )
{
this->IntersectionBuffer[i] = NULL;
this->IntersectionBufferCount[i] = 0;
}
this->SavedTriangleListInput = NULL;
}
// Destructor - release all memory
vtkUnstructuredGridBunykRayCastFunction::~vtkUnstructuredGridBunykRayCastFunction()
{
delete [] this->Points;
this->ClearImage();
delete [] this->Image;
this->Image = NULL;
delete [] this->TetraTriangles;
int i;
for (i = 0; i < 20; i++ )
{
delete [] this->IntersectionBuffer[i];
}
while ( this->TriangleList )
{
Triangle *tmp;
tmp = this->TriangleList->Next;
delete this->TriangleList;
this->TriangleList = tmp;
}
this->ViewToWorldMatrix->Delete();
}
// Clear the intersection image. This does NOT release memory -
// it just sets the link pointers to NULL. The memory is
// contained in the IntersectionBuffer arrays.
void vtkUnstructuredGridBunykRayCastFunction::ClearImage()
{
int i;
if ( this->Image )
{
for ( i = 0; i < this->ImageSize[0]*this->ImageSize[1]; i++ )
{
this->Image[i] = NULL;
}
}
for ( i = 0; i < VTK_BUNYKRCF_MAX_ARRAYS; i++ )
{
this->IntersectionBufferCount[i] = 0;
}
}
// Since we are managing the memory ourself for these intersections,
// we need a new method. In this method we return an unused
// intersection element from our storage arrays. If we don't
// have one, we create a new storage array (unless we have run
// out of memory). The memory can never shrink, and will only be
// deleted when the class is destructed.
void *vtkUnstructuredGridBunykRayCastFunction::NewIntersection()
{
// Look for the first buffer that has enough space, or the
// first one that has not yet been allocated
int i;
for ( i = 0; i < VTK_BUNYKRCF_MAX_ARRAYS; i++ )
{
if ( !this->IntersectionBuffer[i] ||
this->IntersectionBufferCount[i] < VTK_BUNYKRCF_ARRAY_SIZE )
{
break;
}
}
// We have run out of space - return NULL
if ( i == VTK_BUNYKRCF_MAX_ARRAYS )
{
vtkErrorMacro("Out of space for intersections!");
return NULL;
}
// We need another array - allocate it and set its count to 0 indicating
// that we have not used any elements yet
if ( !this->IntersectionBuffer[i] )
{
this->IntersectionBuffer[i] = new Intersection[VTK_BUNYKRCF_ARRAY_SIZE];
this->IntersectionBufferCount[i] = 0;
}
// Return the first unused element
return (this->IntersectionBuffer[i] + (this->IntersectionBufferCount[i]++));
}
// The Intialize method is called from the ray caster at the start of
// rendering. In this method we check if the render is valid (there is
// a renderer, a volume, a mapper, input, etc). We build the basic
// structured if necessary. Then we compute the view dependent information
// such as plane equations and barycentric coordinates per triangle,
// tranfromed points in view space, and the intersection list per pixel.
void vtkUnstructuredGridBunykRayCastFunction::Initialize( vtkRenderer *ren,
vtkVolume *vol )
{
// Check if this is a valid render - we have all the required info
// such as the volume, renderer, mapper, input, etc.
this->Valid = this->CheckValidity( ren, vol );
if ( !this->Valid )
{
return;
}
// Cache some objects for later use during rendering
this->Mapper = vtkUnstructuredGridVolumeRayCastMapper::SafeDownCast( vol->GetMapper() );
this->Renderer = ren;
this->Volume = vol;
vtkUnstructuredGrid *input = this->Mapper->GetInput();
int numPoints = input->GetNumberOfPoints();
// If the number of points have changed, recreate the structure
if ( numPoints != this->NumberOfPoints )
{
delete [] this->Points;
this->Points = new double[3*numPoints];
this->NumberOfPoints = numPoints;
}
// Get the image size from the ray cast mapper. The ImageViewportSize is
// the size of the whole viewport (this does not necessary equal pixel
// size since we may be over / undersampling on the image plane). The size
// (which will be stored in ImageSize) and the ImageOrigin represent the
// subregion of the whole image that we will be considering.
int size[2];
this->Mapper->GetImageInUseSize( size );
this->Mapper->GetImageOrigin( this->ImageOrigin );
this->Mapper->GetImageViewportSize( this->ImageViewportSize );
// If our intersection image is not the right size, recreate it.
// Clear out any old intersections
this->ClearImage();
if ( this->ImageSize[0]*this->ImageSize[1] != size[0]*size[1] )
{
delete [] this->Image;
this->Image = new Intersection *[size[0]*size[1]];
this->ImageSize[0] = size[0];
this->ImageSize[1] = size[1];
this->ClearImage();
}
// Tranform the points. As a by product, compute the ViewToWorldMatrix
// that will be used later
this->TransformPoints();
// If it has not yet been built, or the data has changed in
// some way, we will need to recreate the triangle list. This is
// view independent - although we will leave space in the structure
// for the view dependent info
this->UpdateTriangleList();
// For each triangle store the plane equation and barycentric
// coefficients to be used to speed up rendering
this->ComputeViewDependentInfo();
// Project each boundary triangle onto the image and store intersections
// sorted by depth
this->ComputePixelIntersections();
}
int vtkUnstructuredGridBunykRayCastFunction::CheckValidity( vtkRenderer *ren,
vtkVolume *vol )
{
// We must have a renderer
if ( !ren )
{
vtkErrorMacro("No Renderer");
return 0;
}
// We must have a volume
if ( !vol )
{
vtkErrorMacro("No Volume");
return 0;
}
// We must have a mapper of the correct type
vtkUnstructuredGridVolumeRayCastMapper *mapper =
vtkUnstructuredGridVolumeRayCastMapper::SafeDownCast( vol->GetMapper() );
if ( !mapper )
{
vtkErrorMacro("No mapper or wrong type");
return 0;
}
// The mapper must have input
vtkUnstructuredGrid *input = mapper->GetInput();
if ( !input )
{
vtkErrorMacro("No input to mapper");
return 0;
}
// The input must have some points. This is a silent
// error - just render nothing if it occurs.
int numPoints = input->GetNumberOfPoints();
if ( numPoints == 0 )
{
this->Valid = 0;
return 0;
}
return 1;
}
// This is done once per render - transform the points into view coordinate.
// We also compute the ViewToWorldMatrix here (by inverting the matrix
// we use to project to view coordinates) so that later on in the
// rendering process we can convert points back to world coordinates.
void vtkUnstructuredGridBunykRayCastFunction::TransformPoints()
{
vtkRenderer *ren = this->Renderer;
vtkVolume *vol = this->Volume;
ren->ComputeAspect();
double *aspect = ren->GetAspect();
vtkTransform *perspectiveTransform = vtkTransform::New();
vtkMatrix4x4 *perspectiveMatrix = vtkMatrix4x4::New();
// Get the view matrix in two steps - there is a one step method in camera
// but it turns off stereo so we do not want to use that one
vtkCamera *cam = ren->GetActiveCamera();
perspectiveTransform->Identity();
perspectiveTransform->
Concatenate(cam->GetPerspectiveTransformMatrix(aspect[0]/
aspect[1], 0.0, 1.0 ));
perspectiveTransform->Concatenate(cam->GetViewTransformMatrix());
perspectiveTransform->Concatenate(vol->GetMatrix());
perspectiveMatrix->DeepCopy(perspectiveTransform->GetMatrix());
// Invert this project matrix and store for later use
this->ViewToWorldMatrix->DeepCopy(perspectiveTransform->GetMatrix());
this->ViewToWorldMatrix->Invert();
double *origPtr;
double *transformedPtr = this->Points;
double in[4], out[4];
in[3] = 1.0;
vtkUnstructuredGrid *input = this->Mapper->GetInput();
int numPoints = input->GetNumberOfPoints();
// Loop through all the points and transform them
for ( int i = 0; i < numPoints; i++ )
{
origPtr = input->GetPoint( i );
in[0] = origPtr[0];
in[1] = origPtr[1];
in[2] = origPtr[2];
perspectiveMatrix->MultiplyPoint( in, out );
transformedPtr[0] = (out[0]/out[3] + 1.0)/2.0 *
(double)this->ImageViewportSize[0] - this->ImageOrigin[0];
transformedPtr[1] = (out[1]/out[3] + 1.0)/2.0 *
(double)this->ImageViewportSize[1] - this->ImageOrigin[1];
transformedPtr[2] = out[2]/out[3];
transformedPtr += 3;
}
perspectiveTransform->Delete();
perspectiveMatrix->Delete();
}
// This is done once per change in the data - build a list of
// enumerated triangles (up to four per tetra). Don't store
// duplicates, so we'll have to search for them.
void vtkUnstructuredGridBunykRayCastFunction::UpdateTriangleList()
{
int needsUpdate = 0;
// If we have never created the list, we need updating
if ( !this->TriangleList )
{
needsUpdate = 1;
}
// If the data has changed in some way then we need to update
vtkUnstructuredGrid *input = this->Mapper->GetInput();
if ( this->SavedTriangleListInput != input ||
input->GetMTime() > this->SavedTriangleListMTime.GetMTime() )
{
needsUpdate = 1;
}
// If we don't need updating, return
if ( !needsUpdate )
{
return;
}
// Clear out the old triangle list
while ( this->TriangleList )
{
Triangle *tmp;
tmp = this->TriangleList->Next;
delete this->TriangleList;
this->TriangleList = tmp;
}
this->TriangleList = NULL;
// A temporary structure to reduce search time - VTK_BUNYKRCF_NUMLISTS small
// lists instead of one big one
Triangle *tmpList[VTK_BUNYKRCF_NUMLISTS];
vtkIdType i;
for ( i = 0; i < VTK_BUNYKRCF_NUMLISTS; i++ )
{
tmpList[i] = NULL;
}
int numCells = input->GetNumberOfCells();
// Provide a warning if we find anything other than tetra
int warningNeeded = 0;
// Create a set of links from each tetra to the four triangles
// This is redundant information, but saves time during rendering
this->TetraTriangles = new Triangle *[4 * numCells];
// Loop through all the cells
for ( i = 0; i < numCells; i++ )
{
// We only handle tetra
if ( input->GetCellType(i) != VTK_TETRA )
{
warningNeeded = 1;
continue;
}
// Get the cell
vtkCell *cell = input->GetCell(i);
// Get the four points
vtkIdType pts[4];
pts[0] = cell->GetPointId(0);
pts[1] = cell->GetPointId(1);
pts[2] = cell->GetPointId(2);
pts[3] = cell->GetPointId(3);
// Build each of the four triangles
int ii, jj;
for ( jj = 0; jj < 4; jj++ )
{
vtkIdType tri[3];
int idx = 0;
for ( ii = 0; ii < 4; ii++ )
{
if ( ii != jj )
{
tri[idx++] = pts[ii];
}
}
if ( tri[0] > tri[1] )
{
vtkIdType tmptri = tri[0];
tri[0] = tri[1];
tri[1] = tmptri;
}
if ( tri[1] > tri[2] )
{
vtkIdType tmptri = tri[1];
tri[1] = tri[2];
tri[2] = tmptri;
}
if ( tri[0] > tri[1] )
{
vtkIdType tmptri = tri[0];
tri[0] = tri[1];
tri[1] = tmptri;
}
// Do we have this triangle already?
Triangle *triPtr = tmpList[tri[0]%VTK_BUNYKRCF_NUMLISTS];
while ( triPtr )
{
if ( triPtr->PointIndex[0] == tri[0] &&
triPtr->PointIndex[1] == tri[1] &&
triPtr->PointIndex[2] == tri[2] )
{
break;
}
triPtr = triPtr->Next;
}
if ( triPtr )
{
if ( triPtr->ReferredByTetra[1] != -1 )
{
vtkErrorMacro("Degenerate topology - cell face used more than twice");
}
triPtr->ReferredByTetra[1] = i;
this->TetraTriangles[i*4+jj] = triPtr;
}
else
{
Triangle *next = new Triangle;
next->PointIndex[0] = tri[0];
next->PointIndex[1] = tri[1];
next->PointIndex[2] = tri[2];
next->ReferredByTetra[0] = i;
next->ReferredByTetra[1] = -1;
next->Next = tmpList[tri[0]%VTK_BUNYKRCF_NUMLISTS];
tmpList[tri[0]%VTK_BUNYKRCF_NUMLISTS] = next;
this->TetraTriangles[i*4+jj] = next;
}
}
}
if ( warningNeeded )
{
vtkWarningMacro("Input contains more than tetrahedra - only tetrahedra are supported");
}
// Put the list together
for ( i = 0; i < VTK_BUNYKRCF_NUMLISTS; i++ )
{
if ( tmpList[i] )
{
Triangle *last = tmpList[i];
while ( last->Next )
{
last = last->Next;
}
last->Next = this->TriangleList;
this->TriangleList = tmpList[i];
}
}
this->SavedTriangleListInput = input;
this->SavedTriangleListMTime.Modified();
}
void vtkUnstructuredGridBunykRayCastFunction::ComputeViewDependentInfo()
{
Triangle *triPtr = this->TriangleList;
while ( triPtr )
{
double P1[3], P2[3];
double A[3], B[3], C[3];
A[0] = this->Points[3*triPtr->PointIndex[0]];
A[1] = this->Points[3*triPtr->PointIndex[0]+1];
A[2] = this->Points[3*triPtr->PointIndex[0]+2];
B[0] = this->Points[3*triPtr->PointIndex[1]];
B[1] = this->Points[3*triPtr->PointIndex[1]+1];
B[2] = this->Points[3*triPtr->PointIndex[1]+2];
C[0] = this->Points[3*triPtr->PointIndex[2]];
C[1] = this->Points[3*triPtr->PointIndex[2]+1];
C[2] = this->Points[3*triPtr->PointIndex[2]+2];
P1[0] = B[0] - A[0];
P1[1] = B[1] - A[1];
P1[2] = B[2] - A[2];
P2[0] = C[0] - A[0];
P2[1] = C[1] - A[1];
P2[2] = C[2] - A[2];
triPtr->Denominator = P1[0]*P2[1] - P2[0]*P1[1];
if ( triPtr->Denominator < 0 )
{
double T[3];
triPtr->Denominator = -triPtr->Denominator;
T[0] = P1[0];
T[1] = P1[1];
T[2] = P1[2];
P1[0] = P2[0];
P1[1] = P2[1];
P1[2] = P2[2];
P2[0] = T[0];
P2[1] = T[1];
P2[2] = T[2];
vtkIdType tmpIndex = triPtr->PointIndex[1];
triPtr->PointIndex[1] = triPtr->PointIndex[2];
triPtr->PointIndex[2] = tmpIndex;
}
triPtr->P1X = P1[0];
triPtr->P1Y = P1[1];
triPtr->P2X = P2[0];
triPtr->P2Y = P2[1];
double result[3];
vtkMath::Cross( P1, P2, result );
triPtr->A = result[0];
triPtr->B = result[1];
triPtr->C = result[2];
triPtr->D = -(A[0]*result[0] + A[1]*result[1] + A[2]*result[2]);
triPtr = triPtr->Next;
}
}
void vtkUnstructuredGridBunykRayCastFunction::ComputePixelIntersections()
{
Triangle *triPtr = this->TriangleList;
while ( triPtr )
{
if ( triPtr->ReferredByTetra[1] == -1 )
{
if ( this->IsTriangleFrontFacing( triPtr, triPtr->ReferredByTetra[0] ) )
{
int minX = static_cast<int>(this->Points[3*triPtr->PointIndex[0]]);
int maxX = minX+1;
int minY = static_cast<int>(this->Points[3*triPtr->PointIndex[0]+1]);
int maxY = minY+1;
int tmp;
tmp = static_cast<int>(this->Points[3*triPtr->PointIndex[1]]);
minX = (tmp<minX)?(tmp):(minX);
maxX = ((tmp+1)>maxX)?(tmp+1):(maxX);
tmp = static_cast<int>(this->Points[3*triPtr->PointIndex[1]+1]);
minY = (tmp<minY)?(tmp):(minY);
maxY = ((tmp+1)>maxY)?(tmp+1):(maxY);
tmp = static_cast<int>(this->Points[3*triPtr->PointIndex[2]]);
minX = (tmp<minX)?(tmp):(minX);
maxX = ((tmp+1)>maxX)?(tmp+1):(maxX);
tmp = static_cast<int>(this->Points[3*triPtr->PointIndex[2]+1]);
minY = (tmp<minY)?(tmp):(minY);
maxY = ((tmp+1)>maxY)?(tmp+1):(maxY);
double minZ = this->Points[3*triPtr->PointIndex[0]+2];
double ftmp;
ftmp = this->Points[3*triPtr->PointIndex[1]+2];
minZ = (ftmp<minZ)?(ftmp):(minZ);
ftmp = this->Points[3*triPtr->PointIndex[2]+2];
minZ = (ftmp<minZ)?(ftmp):(minZ);
if ( minX < this->ImageSize[0] - 1 &&
minY < this->ImageSize[1] - 1 &&
maxX >= 0 && maxY >= 0 && minZ > 0.0 )
{
minX = (minX<0)?(0):(minX);
maxX = (maxX>(this->ImageSize[0]-1))?(this->ImageSize[0]-1):(maxX);
minY = (minY<0)?(0):(minY);
maxY = (maxY>(this->ImageSize[1]-1))?(this->ImageSize[1]-1):(maxY);
int x, y;
double ax, ay, az;
ax = this->Points[3*triPtr->PointIndex[0]];
ay = this->Points[3*triPtr->PointIndex[0]+1];
az = this->Points[3*triPtr->PointIndex[0]+2];
for ( y = minY; y <= maxY; y++ )
{
double qy = (double)y - ay;
for ( x = minX; x <= maxX; x++ )
{
double qx = (double)x - ax;
if ( this->InTriangle( qx, qy, triPtr ) )
{
Intersection *intersect = (Intersection *)this->NewIntersection();
if ( intersect )
{
intersect->TriPtr = triPtr;
intersect->Z = az;
intersect->Next = NULL;
if ( !this->Image[y*this->ImageSize[0] + x] ||
intersect->Z < this->Image[y*this->ImageSize[0] + x]->Z )
{
intersect->Next = this->Image[y*this->ImageSize[0] + x];
this->Image[y*this->ImageSize[0] + x] = intersect;
}
else
{
Intersection *test = this->Image[y*this->ImageSize[0] + x];
while ( test->Next && intersect->Z > test->Next->Z )
{
test = test->Next;
}
Intersection *tmpNext = test->Next;
test->Next = intersect;
intersect->Next = tmpNext;
}
}
}
}
}
}
}
}
triPtr = triPtr->Next;
}
}
// Taken from equation on bottom of left column of page 3 - but note that the
// equation in the paper has a mistake: (q1+q2) must be less than 1 (not denom as
// stated in the paper).
int vtkUnstructuredGridBunykRayCastFunction::InTriangle( double x, double y, Triangle *triPtr )
{
double q1, q2;
q1 = (x*triPtr->P2Y - y*triPtr->P2X) / triPtr->Denominator;
q2 = (y*triPtr->P1X - x*triPtr->P1Y) / triPtr->Denominator;
if ( q1 >= 0 && q2 >= 0 && (q1+q2) <= 1.0 )
{
return 1;
}
else
{
return 0;
}
}
int vtkUnstructuredGridBunykRayCastFunction::IsTriangleFrontFacing( Triangle *triPtr, vtkIdType tetraIndex )
{
vtkCell *cell = this->Mapper->GetInput()->GetCell(tetraIndex);
vtkIdType pts[4];
pts[0] = cell->GetPointId(0);
pts[1] = cell->GetPointId(1);
pts[2] = cell->GetPointId(2);
pts[3] = cell->GetPointId(3);
int i;
for( i = 0; i < 4; i++ )
{
if ( pts[i] != triPtr->PointIndex[0] &&
pts[i] != triPtr->PointIndex[1] &&
pts[i] != triPtr->PointIndex[2] )
{
break;
}
}
double d =
triPtr->A*this->Points[3*pts[i]] +
triPtr->B*this->Points[3*pts[i]+1] +
triPtr->C*this->Points[3*pts[i]+2] +
triPtr->D;
return (d>0);
}
template <class T>
vtkIdType TemplateCastRay(
const T *scalars,
vtkUnstructuredGridBunykRayCastFunction *self,
int numComponents,
int x, int y,
double farClipZ,
vtkUnstructuredGridBunykRayCastFunction::Intersection *&intersectionPtr,
vtkUnstructuredGridBunykRayCastFunction::Triangle *&currentTriangle,
vtkIdType &currentTetra,
vtkIdType *intersectedCells,
double *intersectionLengths,
T *nearIntersections,
T *farIntersections,
int maxNumIntersections)
{
int imageViewportSize[2];
self->GetImageViewportSize( imageViewportSize );
int origin[2];
self->GetImageOrigin( origin );
float fx = x - origin[0];
float fy = y - origin[1];
double *points = self->GetPoints();
vtkUnstructuredGridBunykRayCastFunction::Triangle
**triangles = self->GetTetraTriangles();
vtkMatrix4x4 *viewToWorld = self->GetViewToWorldMatrix();
vtkUnstructuredGridBunykRayCastFunction::Triangle *nextTriangle;
vtkIdType nextTetra;
vtkIdType numIntersections = 0;
double nearZ = VTK_FLOAT_MIN;
double nearPoint[4];
double viewCoords[4];
viewCoords[0] = ((float)x / (float)(imageViewportSize[0]-1)) * 2.0 - 1.0;
viewCoords[1] = ((float)y / (float)(imageViewportSize[1]-1)) * 2.0 - 1.0;
// viewCoords[2] set when an intersection is found.
viewCoords[3] = 1.0;
if (currentTriangle)
{
// Find intersection in currentTriangle (the entry point).
nearZ = -( fx*currentTriangle->A + fy*currentTriangle->B +
currentTriangle->D) / currentTriangle->C;
viewCoords[2] = nearZ;
viewToWorld->MultiplyPoint( viewCoords, nearPoint );
nearPoint[0] /= nearPoint[3];
nearPoint[1] /= nearPoint[3];
nearPoint[2] /= nearPoint[3];
}
while (numIntersections < maxNumIntersections)
{
// If we have exited the mesh (or are entering it for the first time,
// find the next intersection with an external face (which has already
// been found with rasterization).
if (!currentTriangle)
{
if (!intersectionPtr)
{
break; // No more intersections.
}
currentTriangle = intersectionPtr->TriPtr;
currentTetra = intersectionPtr->TriPtr->ReferredByTetra[0];
intersectionPtr = intersectionPtr->Next;
// Find intersection in currentTriangle (the entry point).
nearZ = -( fx*currentTriangle->A + fy*currentTriangle->B +
currentTriangle->D) / currentTriangle->C;
viewCoords[2] = nearZ;
viewToWorld->MultiplyPoint( viewCoords, nearPoint );
nearPoint[0] /= nearPoint[3];
nearPoint[1] /= nearPoint[3];
nearPoint[2] /= nearPoint[3];
}
// Find all triangles that the ray may exit.
vtkUnstructuredGridBunykRayCastFunction::Triangle *candidate[3];
int index = 0;
int i;
for ( i = 0; i < 4; i++ )
{
if ( triangles[currentTetra*4+i] != currentTriangle )
{
if ( index == 3 )
{
vtkGenericWarningMacro( "Ugh - found too many triangles!" );
}
else
{
candidate[index++] = triangles[currentTetra*4+i];
}
}
}
double farZ = VTK_FLOAT_MAX;
int minIdx = -1;
// Determine which face the ray exits the cell from.
for ( i = 0; i < 3; i++ )
{
// Far intersection is the nearest intersectation that is farther
// than nearZ.
double tmpZ = 1.0;
if (candidate[i]->C != 0.0)
{
tmpZ =
-( fx*candidate[i]->A +
fy*candidate[i]->B +
candidate[i]->D) / candidate[i]->C;
}
if (tmpZ > nearZ && tmpZ < farZ)
{
farZ = tmpZ;
minIdx = i;
}
}
if (farZ > farClipZ)
{
// Exit happened after point of interest. Bail out now (in case
// we wish to restart).
return numIntersections;
}
if (minIdx == -1)
{
// The ray never exited the cell? Perhaps numerical inaccuracies
// got us here. Just bail out as if we exited the mesh.
nextTriangle = NULL;
nextTetra = -1;
}
else
{
if (intersectedCells)
{
intersectedCells[numIntersections] = currentTetra;
}
nextTriangle = candidate[minIdx];
// Compute intersection with exiting face.
double farPoint[4];
viewCoords[2] = farZ;
viewToWorld->MultiplyPoint( viewCoords, farPoint );
farPoint[0] /= farPoint[3];
farPoint[1] /= farPoint[3];
farPoint[2] /= farPoint[3];
double dist
= sqrt( (nearPoint[0]-farPoint[0])*(nearPoint[0]-farPoint[0])
+ (nearPoint[1]-farPoint[1])*(nearPoint[1]-farPoint[1])
+ (nearPoint[2]-farPoint[2])*(nearPoint[2]-farPoint[2]) );
if (intersectionLengths)
{
intersectionLengths[numIntersections] = dist;
}
// compute the barycentric weights
float ax, ay;
double a1, b1, c1;
ax = points[3*currentTriangle->PointIndex[0]];
ay = points[3*currentTriangle->PointIndex[0]+1];
b1 = ((fx-ax)*currentTriangle->P2Y - (fy-ay)*currentTriangle->P2X) / currentTriangle->Denominator;
c1 = ((fy-ay)*currentTriangle->P1X - (fx-ax)*currentTriangle->P1Y) / currentTriangle->Denominator;
a1 = 1.0 - b1 - c1;
double a2, b2, c2;
ax = points[3*nextTriangle->PointIndex[0]];
ay = points[3*nextTriangle->PointIndex[0]+1];
b2 = ((fx-ax)*nextTriangle->P2Y - (fy-ay)*nextTriangle->P2X) / nextTriangle->Denominator;
c2 = ((fy-ay)*nextTriangle->P1X - (fx-ax)*nextTriangle->P1Y) / nextTriangle->Denominator;
a2 = 1.0 - b2 - c2;
if (nearIntersections)
{
for (int c = 0; c < numComponents; c++)
{
double A, B, C;
A = *(scalars + numComponents*currentTriangle->PointIndex[0] + c);
B = *(scalars + numComponents*currentTriangle->PointIndex[1] + c);
C = *(scalars + numComponents*currentTriangle->PointIndex[2] + c);
nearIntersections[numComponents*numIntersections + c]
= static_cast<T>(a1 * A + b1 * B + c1 * C);
}
}
if (farIntersections)
{
for (int c = 0; c < numComponents; c++)
{
double A, B, C;
A = *(scalars + numComponents*nextTriangle->PointIndex[0] + c);
B = *(scalars + numComponents*nextTriangle->PointIndex[1] + c);
C = *(scalars + numComponents*nextTriangle->PointIndex[2] + c);
farIntersections[numComponents*numIntersections + c]
= static_cast<T>(a2 * A + b2 * B + c2 * C);
}
}
numIntersections++;
// The far triangle has one or two tetras in its referred list.
// If one, return -1 for next tetra and NULL for next triangle
// since we are exiting. If two, return the one that isn't the
// current one.
if ( (nextTriangle)->ReferredByTetra[1] == -1 )
{
nextTetra = -1;
nextTriangle = NULL;
}
else
{
if ( nextTriangle->ReferredByTetra[0] == currentTetra )
{
nextTetra = nextTriangle->ReferredByTetra[1];
}
else
{
nextTetra = nextTriangle->ReferredByTetra[0];
}
}
nearZ = farZ;
nearPoint[0] = farPoint[0];
nearPoint[1] = farPoint[1];
nearPoint[2] = farPoint[2];
nearPoint[3] = farPoint[3];
}
currentTriangle = nextTriangle;
currentTetra = nextTetra;
}
return numIntersections;
}
vtkUnstructuredGridVolumeRayCastIterator
*vtkUnstructuredGridBunykRayCastFunction::NewIterator()
{
if (!this->Valid) return NULL;
vtkUnstructuredGridBunykRayCastIterator *iterator
= vtkUnstructuredGridBunykRayCastIterator::New();
iterator->SetRayCastFunction(this);
return iterator;
}
void vtkUnstructuredGridBunykRayCastFunction::Finalize( )
{
this->Renderer = NULL;
this->Volume = NULL;
this->Mapper = NULL;
this->Valid = 0;
}
//----------------------------------------------------------------------------
void vtkUnstructuredGridBunykRayCastFunction::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
// Do not want to print this->ViewToWorldMatrix , this->ImageViewportSize
// this->ScalarOpacityUnitDistance , or this->ImageOrigin - these are
// internal ivar and not part of the public API for this class
}