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
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 *¤tTriangle,
|
|
vtkIdType ¤tTetra,
|
|
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 *¤tTriangle,
|
|
vtkIdType ¤tTetra,
|
|
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
|
|
}
|
|
|
|
|
|
|