/*========================================================================= Program: Visualization Toolkit Module: vtkCellPicker.cxx 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 "vtkCellPicker.h" #include "vtkObjectFactory.h" #include "vtkCommand.h" #include "vtkMath.h" #include "vtkBox.h" #include "vtkPiecewiseFunction.h" #include "vtkPlaneCollection.h" #include "vtkTransform.h" #include "vtkDoubleArray.h" #include "vtkPoints.h" #include "vtkPolygon.h" #include "vtkVoxel.h" #include "vtkGenericCell.h" #include "vtkPointData.h" #include "vtkImageData.h" #include "vtkActor.h" #include "vtkMapper.h" #include "vtkTexture.h" #include "vtkVolume.h" #include "vtkAbstractVolumeMapper.h" #include "vtkVolumeProperty.h" #include "vtkLODProp3D.h" #include "vtkImageMapper3D.h" #include "vtkRenderer.h" #include "vtkCamera.h" #include "vtkAbstractCellLocator.h" vtkStandardNewMacro(vtkCellPicker); //---------------------------------------------------------------------------- vtkCellPicker::vtkCellPicker() { // List of locators for accelerating polydata picking this->Locators = vtkCollection::New(); // For polydata picking this->Cell = vtkGenericCell::New(); this->PointIds = vtkIdList::New(); // For interpolation of volume gradients this->Gradients = vtkDoubleArray::New(); this->Gradients->SetNumberOfComponents(3); this->Gradients->SetNumberOfTuples(8); // Miscellaneous ivars this->Tolerance = 1e-6; this->VolumeOpacityIsovalue = 0.05; this->UseVolumeGradientOpacity = 0; this->PickClippingPlanes = 0; this->PickTextureData = 0; // Clear all info returned by the pick this->ResetCellPickerInfo(); } //---------------------------------------------------------------------------- vtkCellPicker::~vtkCellPicker() { this->Gradients->Delete(); this->Cell->Delete(); this->PointIds->Delete(); this->Locators->Delete(); } //---------------------------------------------------------------------------- void vtkCellPicker::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os,indent); os << indent << "MapperNormal: (" << this->MapperNormal[0] << "," << this->MapperNormal[1] << "," << this->MapperNormal[2] << ")\n"; os << indent << "PickNormal: (" << this->PickNormal[0] << "," << this->PickNormal[1] << "," << this->PickNormal[2] << ")\n"; if ( this->Texture ) { os << indent << "Texture: " << this->Texture << "\n"; } else { os << indent << "Texture: (none)"; } os << indent << "PickTextureData: " << (this->PickTextureData ? "On" : "Off") << "\n"; os << indent << "PointId: " << this->PointId << "\n"; os << indent << "CellId: " << this->CellId << "\n"; os << indent << "SubId: " << this->SubId << "\n"; os << indent << "PCoords: (" << this->PCoords[0] << ", " << this->PCoords[1] << ", " << this->PCoords[2] << ")\n"; os << indent << "PointIJK: (" << this->PointIJK[0] << ", " << this->PointIJK[1] << ", " << this->PointIJK[2] << ")\n"; os << indent << "CellIJK: (" << this->CellIJK[0] << ", " << this->CellIJK[1] << ", " << this->CellIJK[2] << ")\n"; os << indent << "ClippingPlaneId: " << this->ClippingPlaneId << "\n"; os << indent << "PickClippingPlanes: " << (this->PickClippingPlanes ? "On" : "Off") << "\n"; os << indent << "VolumeOpacityIsovalue: " << this->VolumeOpacityIsovalue << "\n"; os << indent << "UseVolumeGradientOpacity: " << (this->UseVolumeGradientOpacity ? "On" : "Off") << "\n"; } //---------------------------------------------------------------------------- void vtkCellPicker::Initialize() { this->ResetPickInfo(); this->Superclass::Initialize(); } //---------------------------------------------------------------------------- void vtkCellPicker::ResetPickInfo() { // First, reset information from the superclass, since // vtkPicker does not have a ResetPickInfo method this->DataSet = 0; this->Mapper = 0; // Reset all the information specific to this class this->ResetCellPickerInfo(); } //---------------------------------------------------------------------------- void vtkCellPicker::ResetCellPickerInfo() { this->Texture = 0; this->ClippingPlaneId = -1; this->PointId = -1; this->CellId = -1; this->SubId = -1; this->PCoords[0] = 0.0; this->PCoords[1] = 0.0; this->PCoords[2] = 0.0; this->CellIJK[0] = 0; this->CellIJK[1] = 0; this->CellIJK[2] = 0; this->PointIJK[0] = 0; this->PointIJK[1] = 0; this->PointIJK[2] = 0; this->MapperNormal[0] = 0.0; this->MapperNormal[1] = 0.0; this->MapperNormal[2] = 1.0; this->PickNormal[0] = 0.0; this->PickNormal[1] = 0.0; this->PickNormal[2] = 1.0; } //---------------------------------------------------------------------------- void vtkCellPicker::AddLocator(vtkAbstractCellLocator *locator) { if (!this->Locators->IsItemPresent(locator)) { this->Locators->AddItem(locator); } } //---------------------------------------------------------------------------- void vtkCellPicker::RemoveLocator(vtkAbstractCellLocator *locator) { this->Locators->RemoveItem(locator); } //---------------------------------------------------------------------------- void vtkCellPicker::RemoveAllLocators() { this->Locators->RemoveAllItems(); } //---------------------------------------------------------------------------- int vtkCellPicker::Pick(double selectionX, double selectionY, double selectionZ, vtkRenderer *renderer) { int pickResult = 0; if ( (pickResult = this->Superclass::Pick(selectionX, selectionY, selectionZ, renderer)) == 0) { // If no pick, set the PickNormal so that it points at the camera vtkCamera *camera = renderer->GetActiveCamera(); double cameraPos[3]; camera->GetPosition(cameraPos); if (camera->GetParallelProjection()) { // For parallel projection, use -ve direction of projection double cameraFocus[3]; camera->GetFocalPoint(cameraFocus); this->PickNormal[0] = cameraPos[0] - cameraFocus[0]; this->PickNormal[1] = cameraPos[1] - cameraFocus[1]; this->PickNormal[2] = cameraPos[2] - cameraFocus[2]; } else { // Get the vector from pick position to the camera this->PickNormal[0] = cameraPos[0] - this->PickPosition[0]; this->PickNormal[1] = cameraPos[1] - this->PickPosition[1]; this->PickNormal[2] = cameraPos[2] - this->PickPosition[2]; } vtkMath::Normalize(this->PickNormal); } return pickResult; } //---------------------------------------------------------------------------- // Tolerance for parametric coordinate matching an intersection with a plane #define VTKCELLPICKER_PLANE_TOL 1e-14 double vtkCellPicker::IntersectWithLine(double p1[3], double p2[3], double tol, vtkAssemblyPath *path, vtkProp3D *prop, vtkAbstractMapper3D *m) { vtkMapper *mapper = 0; vtkAbstractVolumeMapper *volumeMapper = 0; vtkImageMapper3D *imageMapper = 0; double tMin = VTK_DOUBLE_MAX; double t1 = 0.0; double t2 = 1.0; // Clip the ray with the mapper's ClippingPlanes and adjust t1, t2. // This limits the pick search to the inside of the clipped region. int clippingPlaneId = -1; if (m && !this->ClipLineWithPlanes(m, this->Transform->GetMatrix(), p1, p2, t1, t2, clippingPlaneId)) { return VTK_DOUBLE_MAX; } // Initialize the pick position to the frontmost clipping plane if (this->PickClippingPlanes && clippingPlaneId >= 0) { tMin = t1; } // Volume else if ( (volumeMapper = vtkAbstractVolumeMapper::SafeDownCast(m)) ) { tMin = this->IntersectVolumeWithLine(p1, p2, t1, t2, prop, volumeMapper); } // Image else if ( (imageMapper = vtkImageMapper3D::SafeDownCast(m)) ) { tMin = this->IntersectImageWithLine(p1, p2, t1, t2, prop, imageMapper); } // Actor else if ( (mapper = vtkMapper::SafeDownCast(m)) ) { tMin = this->IntersectActorWithLine(p1, p2, t1, t2, tol, prop, mapper); } // Unidentified Prop3D else { tMin = this->IntersectProp3DWithLine(p1, p2, t1, t2, tol, prop, m); } if (tMin < this->GlobalTMin) { this->GlobalTMin = tMin; this->SetPath(path); this->ClippingPlaneId = clippingPlaneId; // If tMin == t1, the pick didn't go past the first clipping plane, // so the position and normal will be set from the clipping plane. if (fabs(tMin - t1) < VTKCELLPICKER_PLANE_TOL && clippingPlaneId >= 0) { this->MapperPosition[0] = p1[0]*(1.0-t1) + p2[0]*t1; this->MapperPosition[1] = p1[1]*(1.0-t1) + p2[1]*t1; this->MapperPosition[2] = p1[2]*(1.0-t1) + p2[2]*t1; double plane[4]; m->GetClippingPlaneInDataCoords( this->Transform->GetMatrix(), clippingPlaneId, plane); vtkMath::Normalize(plane); // Want normal outward from the planes, not inward this->MapperNormal[0] = -plane[0]; this->MapperNormal[1] = -plane[1]; this->MapperNormal[2] = -plane[2]; } // The position comes from the data, so put it into world coordinates this->Transform->TransformPoint(this->MapperPosition, this->PickPosition); this->Transform->TransformNormal(this->MapperNormal, this->PickNormal); } return tMin; } //---------------------------------------------------------------------------- double vtkCellPicker::IntersectActorWithLine(const double p1[3], const double p2[3], double t1, double t2, double tol, vtkProp3D *prop, vtkMapper *mapper) { // This code was taken from the original CellPicker with almost no // modification except for the locator and texture additions. // Intersect each cell with ray. Keep track of one closest to // the eye (within the tolerance tol) and within the clipping range). // Note that we fudge the "closest to" (tMin+this->Tolerance) a little and // keep track of the cell with the best pick based on parametric // coordinate (pick the minimum, maximum parametric distance). This // breaks ties in a reasonable way when cells are the same distance // from the eye (like cells lying on a 2D plane). vtkDataSet *data = mapper->GetInput(); double tMin = VTK_DOUBLE_MAX; double minPCoords[3]; double pDistMin = VTK_DOUBLE_MAX; vtkIdType minCellId = -1; int minSubId = -1; double minXYZ[3]; minXYZ[0] = minXYZ[1] = minXYZ[2] = 0.0; // Polydata has no 3D cells int isPolyData = data->IsA("vtkPolyData"); vtkCollectionSimpleIterator iter; vtkAbstractCellLocator *locator = 0; this->Locators->InitTraversal(iter); while ( (locator = static_cast( this->Locators->GetNextItemAsObject(iter))) ) { if (locator->GetDataSet() == data) { break; } } // Make a new p1 and p2 using the clipped t1 and t2 double q1[3], q2[3]; q1[0] = p1[0]; q1[1] = p1[1]; q1[2] = p1[2]; q2[0] = p2[0]; q2[1] = p2[1]; q2[2] = p2[2]; if (t1 != 0.0 || t2 != 1.0) { for (int j = 0; j < 3; j++) { q1[j] = p1[j]*(1.0 - t1) + p2[j]*t1; q2[j] = p1[j]*(1.0 - t2) + p2[j]*t2; } } // Use the locator if one exists for this data if (locator) { if (!locator->IntersectWithLine(q1, q2, tol, tMin, minXYZ, minPCoords, minSubId, minCellId, this->Cell)) { return VTK_DOUBLE_MAX; } // Stretch tMin out to the original range if (t1 != 0.0 || t2 != 1.0) { tMin = t1*(1.0 - tMin) + t2*tMin; } // If cell is a strip, then replace cell with a sub-cell this->SubCellFromCell(this->Cell, minSubId); } else { vtkIdList *pointIds = this->PointIds; vtkIdType numCells = data->GetNumberOfCells(); for (vtkIdType cellId = 0; cellId < numCells; cellId++) { double t; double x[3]; double pcoords[3]; pcoords[0] = pcoords[1] = pcoords[2] = 0; int newSubId = -1; int numSubIds = 1; // If it is a strip, we need to iterate over the subIds int cellType = data->GetCellType(cellId); int useSubCells = this->HasSubCells(cellType); if (useSubCells) { // Get the pointIds for the strip and the length of the strip data->GetCellPoints(cellId, pointIds); numSubIds = this->GetNumberOfSubCells(pointIds, cellType); } // This will only loop once unless we need to deal with a strip for (int subId = 0; subId < numSubIds; subId++) { if (useSubCells) { // Get a sub-cell from a the strip this->GetSubCell(data, pointIds, subId, cellType, this->Cell); } else { data->GetCell(cellId, this->Cell); } int cellPicked = 0; if (isPolyData) { // Polydata can always be picked with original endpoints cellPicked = this->Cell->IntersectWithLine( const_cast(p1), const_cast(p2), tol, t, x, pcoords, newSubId); } else { // Any 3D cells need to be intersected with a line segment that // has been clipped with the clipping planes, in case one end is // actually inside the cell. cellPicked = this->Cell->IntersectWithLine( q1, q2, tol, t, x, pcoords, newSubId); // Stretch t out to the original range if (t1 != 0.0 || t2 != 1.0) { t = t1*(1.0 - t) + t2*t; } } if (cellPicked && t <= (tMin + this->Tolerance) && t >= t1 && t <= t2) { double pDist = this->Cell->GetParametricDistance(pcoords); if (pDist < pDistMin || (pDist == pDistMin && t < tMin)) { tMin = t; pDistMin = pDist; // save all of these minCellId = cellId; minSubId = newSubId; if (useSubCells) { minSubId = subId; } for (int k = 0; k < 3; k++) { minXYZ[k] = x[k]; minPCoords[k] = pcoords[k]; } } // for all subIds } // if minimum, maximum } // if a close cell } // for all cells } // Do this if a cell was intersected if (minCellId >= 0 && tMin < this->GlobalTMin) { this->ResetPickInfo(); // Get the cell, convert to triangle if it is a strip vtkGenericCell *cell = this->Cell; // If we used a locator, we already have the picked cell if (!locator) { int cellType = data->GetCellType(minCellId); if (this->HasSubCells(cellType)) { data->GetCellPoints(minCellId, this->PointIds); this->GetSubCell(data, this->PointIds, minSubId, cellType, cell); } else { data->GetCell(minCellId, cell); } } // Get the cell weights vtkIdType numPoints = cell->GetNumberOfPoints(); double *weights = new double[numPoints]; for (vtkIdType i = 0; i < numPoints; i++) { weights[i] = 0; } // Get the interpolation weights (point is thrown away) double point[3]; cell->EvaluateLocation(minSubId, minPCoords, point, weights); this->Mapper = mapper; // Get the texture from the actor or the LOD vtkActor *actor = 0; vtkLODProp3D *lodActor = 0; if ( (actor = vtkActor::SafeDownCast(prop)) ) { this->Texture = actor->GetTexture(); } else if ( (lodActor = vtkLODProp3D::SafeDownCast(prop)) ) { int lodId = lodActor->GetPickLODID(); lodActor->GetLODTexture(lodId, &this->Texture); } if (this->PickTextureData && this->Texture) { // Return the texture's image data to the user vtkImageData *image = this->Texture->GetInput(); this->DataSet = image; // Get and check the image dimensions int extent[6]; image->GetExtent(extent); int dimensionsAreValid = 1; int dimensions[3]; for (int i = 0; i < 3; i++) { dimensions[i] = extent[2*i + 1] - extent[2*i] + 1; dimensionsAreValid = (dimensionsAreValid && dimensions[i] > 0); } // Use the texture coord to set the information double tcoord[3]; if (dimensionsAreValid && this->ComputeSurfaceTCoord(data, cell, weights, tcoord)) { // Take the border into account when computing coordinates double x[3]; x[0] = extent[0] + tcoord[0]*dimensions[0] - 0.5; x[1] = extent[2] + tcoord[1]*dimensions[1] - 0.5; x[2] = extent[4] + tcoord[2]*dimensions[2] - 0.5; this->SetImageDataPickInfo(x, extent); } } else { // Return the polydata to the user this->DataSet = data; this->CellId = minCellId; this->SubId = minSubId; this->PCoords[0] = minPCoords[0]; this->PCoords[1] = minPCoords[1]; this->PCoords[2] = minPCoords[2]; // Find the point with the maximum weight double maxWeight = 0; vtkIdType iMaxWeight = -1; for (vtkIdType i = 0; i < numPoints; i++) { if (weights[i] > maxWeight) { iMaxWeight = i; } } // If maximum weight is found, use it to get the PointId if (iMaxWeight != -1) { this->PointId = cell->PointIds->GetId(iMaxWeight); } } // Set the mapper position this->MapperPosition[0] = minXYZ[0]; this->MapperPosition[1] = minXYZ[1]; this->MapperPosition[2] = minXYZ[2]; // Compute the normal if (!this->ComputeSurfaceNormal(data, cell, weights, this->MapperNormal)) { // By default, the normal points back along view ray this->MapperNormal[0] = p1[0] - p2[0]; this->MapperNormal[1] = p1[1] - p2[1]; this->MapperNormal[2] = p1[2] - p2[2]; vtkMath::Normalize(this->MapperNormal); } delete [] weights; } return tMin; } //---------------------------------------------------------------------------- // Intersect a vtkVolume with a line by ray casting. // For algorithm stability: choose a tolerance that is larger than // the expected roundoff error in computing the voxel indices from "t" #define VTKCELLPICKER_VOXEL_TOL 1e-6 double vtkCellPicker::IntersectVolumeWithLine(const double p1[3], const double p2[3], double t1, double t2, vtkProp3D *prop, vtkAbstractVolumeMapper *mapper) { vtkImageData *data = vtkImageData::SafeDownCast(mapper->GetDataSetInput()); if (data == 0) { // This picker only works with image inputs return VTK_DOUBLE_MAX; } // Convert ray to structured coordinates double spacing[3], origin[3]; int extent[6]; data->GetSpacing(spacing); data->GetOrigin(origin); data->GetExtent(extent); double x1[3], x2[3]; for (int i = 0; i < 3; i++) { x1[i] = (p1[i] - origin[i])/spacing[i]; x2[i] = (p2[i] - origin[i])/spacing[i]; } // Clip the ray with the extent, results go in s1 and s2 int planeId; double s1, s2; if (!this->ClipLineWithExtent(extent, x1, x2, s1, s2, planeId)) { return VTK_DOUBLE_MAX; } if (s1 >= t1) { t1 = s1; } if (s2 <= t2) { t2 = s2; } // Sanity check if (t2 < t1) { return VTK_DOUBLE_MAX; } // Get the property from the volume or the LOD vtkVolumeProperty *property = 0; vtkVolume *volume = 0; vtkLODProp3D *lodVolume = 0; if ( (volume = vtkVolume::SafeDownCast(prop)) ) { property = volume->GetProperty(); } else if ( (lodVolume = vtkLODProp3D::SafeDownCast(prop)) ) { int lodId = lodVolume->GetPickLODID(); lodVolume->GetLODProperty(lodId, &property); } // Get the theshold for the opacity double opacityThreshold = this->VolumeOpacityIsovalue; // Compute the length of the line intersecting the volume double rayLength = sqrt(vtkMath::Distance2BetweenPoints(x1, x2))*(t2 - t1); // This is the minimum increment that will be allowed double tTol = VTKCELLPICKER_VOXEL_TOL/rayLength*(t2 - t1); // Find out whether there are multiple components in the volume int numComponents = data->GetNumberOfScalarComponents(); int independentComponents = property->GetIndependentComponents(); int numIndependentComponents = 1; if (independentComponents) { numIndependentComponents = numComponents; } // Create a scalar array, it will be needed later vtkDataArray *scalars = vtkDataArray::CreateDataArray(data->GetScalarType()); scalars->SetNumberOfComponents(numComponents); vtkIdType scalarArraySize = numComponents*data->GetNumberOfPoints(); int scalarSize = data->GetScalarSize(); void *scalarPtr = data->GetScalarPointer(); // Go through each volume component separately double tMin = VTK_DOUBLE_MAX; for (int component = 0; component < numIndependentComponents; component++) { vtkPiecewiseFunction *scalarOpacity = property->GetScalarOpacity(component); int disableGradientOpacity = property->GetDisableGradientOpacity(component); vtkPiecewiseFunction *gradientOpacity = 0; if (!disableGradientOpacity && this->UseVolumeGradientOpacity) { gradientOpacity = property->GetGradientOpacity(component); } // This is the component used to compute the opacity int oComponent = component; if (!independentComponents) { oComponent = numComponents - 1; } // Make a new array, shifted to the desired component scalars->SetVoidArray(static_cast(static_cast(scalarPtr) + scalarSize*oComponent), scalarArraySize, 1); // Do a ray cast with linear interpolation. double opacity = 0.0; double lastOpacity = 0.0; double lastT = t1; double x[3]; double pcoords[3]; int xi[3]; // Ray cast loop double t = t1; while (t <= t2) { for (int j = 0; j < 3; j++) { // "t" is the fractional distance between endpoints x1 and x2 x[j] = x1[j]*(1.0 - t) + x2[j]*t; // Paranoia bounds check if (x[j] < extent[2*j]) { x[j] = extent[2*j]; } else if (x[j] > extent[2*j + 1]) { x[j] = extent[2*j+1]; } xi[j] = vtkMath::Floor(x[j]); pcoords[j] = x[j] - xi[j]; } opacity = this->ComputeVolumeOpacity(xi, pcoords, data, scalars, scalarOpacity, gradientOpacity); // If the ray has crossed the isosurface, then terminate the loop if (opacity > opacityThreshold) { break; } lastT = t; lastOpacity = opacity; // Compute the next "t" value that crosses a voxel boundary t = 1.0; for (int k = 0; k < 3; k++) { // Skip dimension "k" if it is perpendicular to ray if (fabs(x2[k] - x1[k]) > VTKCELLPICKER_VOXEL_TOL*rayLength) { // Compute the previous coord along dimension "k" double lastX = x1[k]*(1.0 - lastT) + x2[k]*lastT; // Increment to next slice boundary along dimension "k", // including a tolerance value for stability in cases // where lastX is just less than an integer value. double nextX = 0; if (x2[k] > x1[k]) { nextX = vtkMath::Floor(lastX + VTKCELLPICKER_VOXEL_TOL) + 1; } else { nextX = vtkMath::Ceil(lastX - VTKCELLPICKER_VOXEL_TOL) - 1; } // Compute the "t" value for this slice boundary double ttry = lastT + (nextX - lastX)/(x2[k] - x1[k]); if (ttry > lastT + tTol && ttry < t) { t = ttry; } } } } // End of "while (t <= t2)" // If the ray hit the isosurface, compute the isosurface position if (opacity > opacityThreshold) { // Backtrack to the actual surface position unless this was first step if (t > t1) { double f = (opacityThreshold - lastOpacity)/(opacity - lastOpacity); t = lastT*(1.0 - f) + t*f; for (int j = 0; j < 3; j++) { x[j] = x1[j]*(1.0 - t) + x2[j]*t; if (x[j] < extent[2*j]) { x[j] = extent[2*j]; } else if (x[j] > extent[2*j + 1]) { x[j] = extent[2*j+1]; } xi[j] = vtkMath::Floor(x[j]); pcoords[j] = x[j] - xi[j]; } } // Check to see if this is the new global minimum if (t < tMin && t < this->GlobalTMin) { this->ResetPickInfo(); tMin = t; this->Mapper = mapper; this->DataSet = data; this->SetImageDataPickInfo(x, extent); this->MapperPosition[0] = x[0]*spacing[0] + origin[0]; this->MapperPosition[1] = x[1]*spacing[1] + origin[1]; this->MapperPosition[2] = x[2]*spacing[2] + origin[2]; // Default the normal to the view-plane normal. This default // will be used if the gradient cannot be computed any other way. this->MapperNormal[0] = p1[0] - p2[0]; this->MapperNormal[1] = p1[1] - p2[1]; this->MapperNormal[2] = p1[2] - p2[2]; vtkMath::Normalize(this->MapperNormal); // Check to see if this is the first step, which means that this // is the boundary of the volume. If this is the case, use the // normal of the boundary. if (t == t1 && planeId >= 0 && xi[planeId/2] == extent[planeId]) { this->MapperNormal[0] = 0.0; this->MapperNormal[1] = 0.0; this->MapperNormal[2] = 0.0; this->MapperNormal[planeId/2] = 2.0*(planeId%2) - 1.0; if (spacing[planeId/2] < 0) { this->MapperNormal[planeId/2] = - this->MapperNormal[planeId/2]; } } else { // Set the normal from the direction of the gradient int *ci = this->CellIJK; double weights[8]; vtkVoxel::InterpolationFunctions(this->PCoords, weights); data->GetVoxelGradient(ci[0], ci[1], ci[2], scalars, this->Gradients); double v[3]; v[0] = v[1] = v[2] = 0.0; for (int k = 0; k < 8; k++) { double *pg = this->Gradients->GetTuple(k); v[0] += pg[0]*weights[k]; v[1] += pg[1]*weights[k]; v[2] += pg[2]*weights[k]; } double norm = vtkMath::Norm(v); if (norm > 0) { this->MapperNormal[0] = v[0]/norm; this->MapperNormal[1] = v[1]/norm; this->MapperNormal[2] = v[2]/norm; } } } // End of "if (opacity > opacityThreshold)" } // End of "if (t < tMin && t < this->GlobalTMin)" } // End of loop over volume components scalars->Delete(); return tMin; } //---------------------------------------------------------------------------- double vtkCellPicker::IntersectImageWithLine(const double p1[3], const double p2[3], double t1, double t2, vtkProp3D *prop, vtkImageMapper3D *imageMapper) { // Get the image information vtkImageData *data = imageMapper->GetInput(); double spacing[3], origin[3]; int extent[6]; data->GetSpacing(spacing); data->GetOrigin(origin); data->GetExtent(extent); // Get the plane equation for the slice double normal[4]; imageMapper->GetSlicePlaneInDataCoords(prop->GetMatrix(), normal); // Point the normal towards camera if (normal[0]*(p1[0] - p2[0]) + normal[1]*(p1[1] - p2[1]) + normal[2]*(p1[2] - p2[2]) < 0) { normal[0] = -normal[0]; normal[1] = -normal[1]; normal[2] = -normal[2]; normal[3] = -normal[3]; } // And finally convert normal to structured coords double xnormal[4]; xnormal[0] = normal[0]*spacing[0]; xnormal[1] = normal[1]*spacing[1]; xnormal[2] = normal[2]*spacing[2]; xnormal[3] = normal[3] + vtkMath::Dot(origin, normal); double l = vtkMath::Norm(xnormal); xnormal[0] /= l; xnormal[1] /= l; xnormal[2] /= l; xnormal[3] /= l; // Also convert ray to structured coords double x1[3], x2[3]; for (int i = 0; i < 3; i++) { x1[i] = (p1[i] - origin[i])/spacing[i]; x2[i] = (p2[i] - origin[i])/spacing[i]; } // Get the bounds to discover any cropping that has been applied double bounds[6]; imageMapper->GetBounds(bounds); // Convert bounds to structured coords for (int k = 0; k < 3; k++) { bounds[2*k] = (bounds[2*k] - origin[k])/spacing[k]; bounds[2*k+1] = (bounds[2*k+1] - origin[k])/spacing[k]; // It should be a multiple of 0.5, so round to closest multiple of 0.5 // (this reduces the impact of roundoff error from the above computation) bounds[2*k] = 0.5*vtkMath::Round(2.0*(bounds[2*k])); bounds[2*k+1] = 0.5*vtkMath::Round(2.0*(bounds[2*k+1])); } // Clip the ray with the extent int planeId, plane2Id; double tMin, tMax; if (!vtkBox::IntersectWithLine(bounds, x1, x2, tMin, tMax, 0, 0, planeId, plane2Id)) { return VTK_DOUBLE_MAX; } if (tMin != tMax) { // Intersect the ray with the slice plane double w1 = vtkMath::Dot(x1, xnormal) + xnormal[3]; double w2 = vtkMath::Dot(x2, xnormal) + xnormal[3]; if (w1*w2 > VTKCELLPICKER_VOXEL_TOL) { return VTK_DOUBLE_MAX; } if (w1*w2 < 0) { tMin = w1/(w1 - w2); } } // Make sure that intersection is within clipping planes if (tMin < t1 || tMin > t2) { return VTK_DOUBLE_MAX; } if (tMin < this->GlobalTMin) { // Compute the pick position in structured coords double x[3]; for (int j = 0; j < 3; j++) { x[j] = x1[j]*(1.0 - tMin) + x2[j]*tMin; // Do a bounds check. If beyond tolerance of bound, then // pick failed, but if within tolerance, clamp the coord // to the bound for robustness against roundoff errors. if (x[j] < bounds[2*j]) { if (x[j] < bounds[2*j] - VTKCELLPICKER_VOXEL_TOL) { return VTK_DOUBLE_MAX; } x[j] = bounds[2*j]; } else if (x[j] > bounds[2*j+1]) { if (x[j] > bounds[2*j+1] + VTKCELLPICKER_VOXEL_TOL) { return VTK_DOUBLE_MAX; } x[j] = bounds[2*j + 1]; } } this->ResetPickInfo(); this->Mapper = 0; this->DataSet = data; // Compute all the pick values this->SetImageDataPickInfo(x, extent); this->MapperPosition[0] = origin[0] + x[0]*spacing[0]; this->MapperPosition[1] = origin[1] + x[1]*spacing[1]; this->MapperPosition[2] = origin[2] + x[2]*spacing[2]; // Set the normal in mapper coordinates this->MapperNormal[0] = normal[0]; this->MapperNormal[1] = normal[1]; this->MapperNormal[2] = normal[2]; } return tMin; } //---------------------------------------------------------------------------- // This is a catch-all for Prop3D types that vtkCellPicker does not // recognize. It can be overridden in subclasses to provide support // for picking new Prop3D types. double vtkCellPicker::IntersectProp3DWithLine(const double *, const double *, double, double, double, vtkProp3D *, vtkAbstractMapper3D *) { return VTK_DOUBLE_MAX; } //---------------------------------------------------------------------------- // Clip a line with a collection of clipping planes, or return zero if // the line does not intersect the volume enclosed by the planes. // The result of the clipping is retured in t1 and t2, which will have // values between 0 and 1. The index of the frontmost intersected plane is // returned in planeId. int vtkCellPicker::ClipLineWithPlanes(vtkAbstractMapper3D *mapper, vtkMatrix4x4 *mat, const double p1[3], const double p2[3], double &t1, double &t2, int& planeId) { // The minPlaneId is the index of the plane that t1 lies on planeId = -1; t1 = 0.0; t2 = 1.0; double plane[4]; int numClipPlanes = mapper->GetNumberOfClippingPlanes(); for (int i = 0; i < numClipPlanes; i++) { mapper->GetClippingPlaneInDataCoords(mat, i, plane); double d1 = plane[0]*p1[0] + plane[1]*p1[1] + plane[2]*p1[2] + plane[3]; double d2 = plane[0]*p2[0] + plane[1]*p2[1] + plane[2]*p2[2] + plane[3]; // If both distances are negative, both points are outside if (d1 < 0 && d2 < 0) { return 0; } // If only one of the distances is negative, the line crosses the plane else if (d1 < 0 || d2 < 0) { // Compute fractional distance "t" of the crossing between p1 & p2 double t = 0.0; // The "if" here just avoids an expensive division when possible if (d1 != 0) { // We will never have d1==d2 since they have different signs t = d1/(d1 - d2); } // If point p1 was clipped, adjust t1 if (d1 < 0) { if (t >= t1) { t1 = t; planeId = i; } } // else point p2 was clipped, so adjust t2 else { if (t <= t2) { t2 = t; } } // If this happens, there's no line left if (t1 > t2) { return 0; } } } return 1; } //---------------------------------------------------------------------------- // Clip a line in structured coordinates with an extent. If the line // does not intersect the extent, the return value will be zero. // The fractional position of the new x1 with respect to the original line // is returned in tMin, and the index of the frontmost intersected plane // is returned in planeId. The planes are ordered as follows: // xmin, xmax, ymin, ymax, zmin, zmax. int vtkCellPicker::ClipLineWithExtent(const int extent[6], const double x1[3], const double x2[3], double &t1, double &t2, int &planeId) { double bounds[6]; bounds[0] = extent[0]; bounds[1] = extent[1]; bounds[2] = extent[2]; bounds[3] = extent[3]; bounds[4] = extent[4]; bounds[5] = extent[5]; int p2; return vtkBox::IntersectWithLine(bounds, x1, x2, t1, t2, 0, 0, planeId, p2); } //---------------------------------------------------------------------------- // Compute the cell normal either by interpolating the point normals, // or by computing the plane normal for 2D cells. int vtkCellPicker::ComputeSurfaceNormal(vtkDataSet *data, vtkCell *cell, const double *weights, double normal[3]) { vtkDataArray *normals = data->GetPointData()->GetNormals(); if (normals) { normal[0] = normal[1] = normal[2] = 0.0; double pointNormal[3]; vtkIdType numPoints = cell->GetNumberOfPoints(); for (vtkIdType k = 0; k < numPoints; k++) { normals->GetTuple(cell->PointIds->GetId(k), pointNormal); normal[0] += pointNormal[0]*weights[k]; normal[1] += pointNormal[1]*weights[k]; normal[2] += pointNormal[2]*weights[k]; } vtkMath::Normalize(normal); } else if (cell->GetCellDimension() == 2) { vtkPolygon::ComputeNormal(cell->Points, normal); } else { return 0; } return 1; } //---------------------------------------------------------------------------- // Use weights to compute the texture coordinates of a point on the cell. int vtkCellPicker::ComputeSurfaceTCoord(vtkDataSet *data, vtkCell *cell, const double *weights, double tcoord[3]) { vtkDataArray *tcoords = data->GetPointData()->GetTCoords(); if (tcoords) { tcoord[0] = tcoord[1] = tcoord[2] = 0.0; double pointTCoord[3]; int numComponents = tcoords->GetNumberOfComponents(); vtkIdType numPoints = cell->GetNumberOfPoints(); for (vtkIdType k = 0; k < numPoints; k++) { tcoords->GetTuple(cell->PointIds->GetId(k), pointTCoord); for (int i = 0; i < numComponents; i++) { tcoord[i] += pointTCoord[i]*weights[k]; } } return 1; } return 0; } //---------------------------------------------------------------------------- // Do an in-place replacement of a cell with a subcell of that cell void vtkCellPicker::SubCellFromCell(vtkGenericCell *cell, int subId) { switch(cell->GetCellType()) { case VTK_TRIANGLE_STRIP: { static int idx[2][3]={{0,1,2},{1,0,2}}; int *order = idx[subId & 1]; vtkIdType pointIds[3]; double points[3][3]; pointIds[0] = cell->PointIds->GetId(subId + order[0]); pointIds[1] = cell->PointIds->GetId(subId + order[1]); pointIds[2] = cell->PointIds->GetId(subId + order[2]); cell->Points->GetPoint(subId + order[0], points[0]); cell->Points->GetPoint(subId + order[1], points[1]); cell->Points->GetPoint(subId + order[2], points[2]); cell->SetCellTypeToTriangle(); cell->PointIds->SetId(0, pointIds[0]); cell->PointIds->SetId(1, pointIds[1]); cell->PointIds->SetId(2, pointIds[2]); cell->Points->SetPoint(0, points[0]); cell->Points->SetPoint(1, points[1]); cell->Points->SetPoint(2, points[2]); } break; case VTK_POLY_LINE: { vtkIdType pointIds[2]; double points[2][3]; pointIds[0] = cell->PointIds->GetId(subId); pointIds[1] = cell->PointIds->GetId(subId + 1); cell->Points->GetPoint(subId, points[0]); cell->Points->GetPoint(subId + 1, points[1]); cell->SetCellTypeToLine(); cell->PointIds->SetId(0, pointIds[0]); cell->PointIds->SetId(1, pointIds[1]); cell->Points->SetPoint(0, points[0]); cell->Points->SetPoint(1, points[1]); } break; case VTK_POLY_VERTEX: { double point[3]; vtkIdType pointId = cell->PointIds->GetId(subId); cell->Points->GetPoint(subId, point); cell->SetCellTypeToVertex(); cell->PointIds->SetId(0, pointId); cell->Points->SetPoint(0, point); } break; } } //---------------------------------------------------------------------------- int vtkCellPicker::HasSubCells(int cellType) { switch (cellType) { case VTK_TRIANGLE_STRIP: case VTK_POLY_LINE: case VTK_POLY_VERTEX: return 1; break; } return 0; } //---------------------------------------------------------------------------- // Extract a single subcell from a cell in a data set int vtkCellPicker::GetNumberOfSubCells(vtkIdList *pointIds, int cellType) { switch (cellType) { case VTK_TRIANGLE_STRIP: return pointIds->GetNumberOfIds() - 2; break; case VTK_POLY_LINE: return pointIds->GetNumberOfIds() - 1; break; case VTK_POLY_VERTEX: return pointIds->GetNumberOfIds(); break; } return 0; } //---------------------------------------------------------------------------- // Extract a single subcell from a cell in a data set. This method // requires a vtkIdList that contains the pointIds for the cell. void vtkCellPicker::GetSubCell(vtkDataSet *data, vtkIdList *ptIds, int subId, int cellType, vtkGenericCell *cell) { switch(cellType) { case VTK_TRIANGLE_STRIP: { static int idx[2][3]={{0,1,2},{1,0,2}}; int *order = idx[subId & 1]; vtkIdType pointIds[3]; double points[3][3]; pointIds[0] = ptIds->GetId(subId + order[0]); pointIds[1] = ptIds->GetId(subId + order[1]); pointIds[2] = ptIds->GetId(subId + order[2]); data->GetPoint(pointIds[0], points[0]); data->GetPoint(pointIds[1], points[1]); data->GetPoint(pointIds[2], points[2]); cell->SetCellTypeToTriangle(); cell->PointIds->SetId(0, pointIds[0]); cell->PointIds->SetId(1, pointIds[1]); cell->PointIds->SetId(2, pointIds[2]); cell->Points->SetPoint(0, points[0]); cell->Points->SetPoint(1, points[1]); cell->Points->SetPoint(2, points[2]); } break; case VTK_POLY_LINE: { vtkIdType pointIds[2]; double points[2][3]; pointIds[0] = ptIds->GetId(subId); pointIds[1] = ptIds->GetId(subId + 1); data->GetPoint(pointIds[0], points[0]); data->GetPoint(pointIds[1], points[1]); cell->SetCellTypeToLine(); cell->PointIds->SetId(0, pointIds[0]); cell->PointIds->SetId(1, pointIds[1]); cell->Points->SetPoint(0, points[0]); cell->Points->SetPoint(1, points[1]); } break; case VTK_POLY_VERTEX: { double point[3]; vtkIdType pointId = ptIds->GetId(subId); data->GetPoint(pointId, point); cell->SetCellTypeToVertex(); cell->PointIds->SetId(0, pointId); cell->Points->SetPoint(0, point); } break; } } //---------------------------------------------------------------------------- // Set all Cell and Point information, given a structured coordinate // and the extent of the data. void vtkCellPicker::SetImageDataPickInfo(const double x[3], const int extent[6]) { for (int j = 0; j < 3; j++) { double xj = x[j]; if (xj < extent[2*j]) { xj = extent[2*j]; } if (xj > extent[2*j+1]) { xj = extent[2*j+1]; } this->CellIJK[j] = vtkMath::Floor(xj); this->PCoords[j] = xj - this->CellIJK[j]; // Keep the cell in-bounds if it is on the edge if (this->CellIJK[j] == extent[2*j+1] && this->CellIJK[j] > extent[2*j]) { this->CellIJK[j] -= 1; this->PCoords[j] = 1.0; } this->PointIJK[j] = this->CellIJK[j] + (this->PCoords[j] >= 0.5); } // Stupid const/non-const, and I hate const_cast int ext[6]; ext[0] = extent[0]; ext[1] = extent[1]; ext[2] = extent[2]; ext[3] = extent[3]; ext[4] = extent[4]; ext[5] = extent[5]; this->PointId = vtkStructuredData::ComputePointIdForExtent(ext, this->PointIJK); this->CellId = vtkStructuredData::ComputeCellIdForExtent(ext, this->CellIJK); this->SubId = 0; } //---------------------------------------------------------------------------- // Given a structured position within the volume, and the point scalars, // compute the local opacity of the volume. double vtkCellPicker::ComputeVolumeOpacity( const int xi[3], const double pcoords[3], vtkImageData *data, vtkDataArray *scalars, vtkPiecewiseFunction *scalarOpacity, vtkPiecewiseFunction *gradientOpacity) { double opacity = 1.0; // Get interpolation weights from the pcoords double weights[8]; vtkVoxel::InterpolationFunctions(const_cast(pcoords), weights); // Get the volume extent to avoid out-of-bounds int extent[6]; data->GetExtent(extent); int scalarType = data->GetScalarType(); // Compute the increments for the three directions, checking the bounds vtkIdType xInc = 1; vtkIdType yInc = extent[1] - extent[0] + 1; vtkIdType zInc = yInc*(extent[3] - extent[2] + 1); if (xi[0] == extent[1]) { xInc = 0; } if (xi[1] == extent[3]) { yInc = 0; } if (xi[2] == extent[5]) { zInc = 0; } // Use the increments and weights to interpolate the data vtkIdType ptId = data->ComputePointId(const_cast(xi)); double val = 0.0; for (int j = 0; j < 8; j++) { vtkIdType ptInc = (j & 1)*xInc + ((j>>1) & 1)*yInc + ((j>>2) & 1)*zInc; val += weights[j]*scalars->GetComponent(ptId + ptInc, 0); } // Compute the ScalarOpacity if (scalarOpacity) { opacity *= scalarOpacity->GetValue(val); } else if (scalarType == VTK_FLOAT || scalarType == VTK_DOUBLE) { opacity *= val; } else { // Assume unsigned char opacity *= val/255.0; } // Compute gradient and GradientOpacity if (gradientOpacity) { data->GetVoxelGradient(xi[0], xi[1], xi[2], scalars, this->Gradients); double v[3]; v[0] = v[1] = v[2] = 0.0; for (int k = 0; k < 8; k++) { double *pg = this->Gradients->GetTuple(k); v[0] += pg[0]*weights[k]; v[1] += pg[1]*weights[k]; v[2] += pg[2]*weights[k]; } double grad = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); opacity *= gradientOpacity->GetValue(grad); } return opacity; }