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