/*========================================================================= Program: Visualization Toolkit Module: vtkVolumeRayCastMapper.cxx Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen All rights reserved. See Copyright.txt or http://www.kitware.com/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notice for more information. =========================================================================*/ #include "vtkVolumeRayCastMapper.h" #include "vtkCamera.h" #include "vtkDataArray.h" #include "vtkEncodedGradientEstimator.h" #include "vtkEncodedGradientShader.h" #include "vtkFiniteDifferenceGradientEstimator.h" #include "vtkGarbageCollector.h" #include "vtkGraphicsFactory.h" #include "vtkImageData.h" #include "vtkMath.h" #include "vtkMultiThreader.h" #include "vtkObjectFactory.h" #include "vtkPlaneCollection.h" #include "vtkPointData.h" #include "vtkRenderWindow.h" #include "vtkRenderer.h" #include "vtkTimerLog.h" #include "vtkTransform.h" #include "vtkVolumeProperty.h" #include "vtkVolumeRayCastFunction.h" #include "vtkRayCastImageDisplayHelper.h" #include vtkStandardNewMacro(vtkVolumeRayCastMapper); vtkCxxSetObjectMacro(vtkVolumeRayCastMapper,VolumeRayCastFunction, vtkVolumeRayCastFunction ); #define vtkVRCMultiplyPointMacro( 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 vtkVRCMultiplyViewPointMacro( 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[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]; } #define vtkVRCMultiplyNormalMacro( 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] // Construct a new vtkVolumeRayCastMapper with default values vtkVolumeRayCastMapper::vtkVolumeRayCastMapper() { this->SampleDistance = 1.0; this->ImageSampleDistance = 1.0; this->MinimumImageSampleDistance = 1.0; this->MaximumImageSampleDistance = 10.0; this->AutoAdjustSampleDistances = 1; this->VolumeRayCastFunction = NULL; this->GradientEstimator = vtkFiniteDifferenceGradientEstimator::New(); this->GradientShader = vtkEncodedGradientShader::New(); 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->VolumeMatrix = vtkMatrix4x4::New(); this->PerspectiveTransform = vtkTransform::New(); this->VoxelsTransform = vtkTransform::New(); this->VoxelsToViewTransform = vtkTransform::New(); this->ImageMemorySize[0] = 0; this->ImageMemorySize[1] = 0; this->Threader = vtkMultiThreader::New(); this->Image = NULL; this->RowBounds = NULL; this->OldRowBounds = NULL; this->RenderTimeTable = NULL; this->RenderVolumeTable = NULL; this->RenderRendererTable = NULL; this->RenderTableSize = 0; this->RenderTableEntries = 0; this->ZBuffer = NULL; this->ZBufferSize[0] = 0; this->ZBufferSize[1] = 0; this->ZBufferOrigin[0] = 0; this->ZBufferOrigin[1] = 0; this->ImageDisplayHelper = vtkRayCastImageDisplayHelper::New(); this->IntermixIntersectingGeometry = 1; } // Destruct a vtkVolumeRayCastMapper - clean up any memory used vtkVolumeRayCastMapper::~vtkVolumeRayCastMapper() { if ( this->GradientEstimator ) { this->GradientEstimator->UnRegister(this); this->GradientEstimator = NULL; } this->GradientShader->Delete(); this->SetVolumeRayCastFunction(NULL); this->PerspectiveMatrix->Delete(); this->ViewToWorldMatrix->Delete(); this->ViewToVoxelsMatrix->Delete(); this->VoxelsToViewMatrix->Delete(); this->WorldToVoxelsMatrix->Delete(); this->VoxelsToWorldMatrix->Delete(); this->VolumeMatrix->Delete(); this->VoxelsTransform->Delete(); this->VoxelsToViewTransform->Delete(); this->PerspectiveTransform->Delete(); this->ImageDisplayHelper->Delete(); this->Threader->Delete(); if ( this->Image ) { delete [] this->Image; } if ( this->RenderTableSize ) { delete [] this->RenderTimeTable; delete [] this->RenderVolumeTable; delete [] this->RenderRendererTable; } if ( this->RowBounds ) { delete [] this->RowBounds; delete [] this->OldRowBounds; } } float vtkVolumeRayCastMapper::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 vtkVolumeRayCastMapper::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 vtkVolumeRayCastMapper::SetNumberOfThreads( int num ) { this->Threader->SetNumberOfThreads( num ); } int vtkVolumeRayCastMapper::GetNumberOfThreads() { if (this->Threader) { return this->Threader->GetNumberOfThreads(); } return 0; } void vtkVolumeRayCastMapper::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(); } float vtkVolumeRayCastMapper::GetGradientMagnitudeScale() { if ( !this->GradientEstimator ) { vtkErrorMacro( "You must have a gradient estimator set to get the scale" ); return 1.0; } return this->GradientEstimator->GetGradientMagnitudeScale(); } float vtkVolumeRayCastMapper::GetGradientMagnitudeBias() { if ( !this->GradientEstimator ) { vtkErrorMacro( "You must have a gradient estimator set to get the bias" ); return 1.0; } return this->GradientEstimator->GetGradientMagnitudeBias(); } void vtkVolumeRayCastMapper::ReleaseGraphicsResources(vtkWindow *) { } void vtkVolumeRayCastMapper::Render( vtkRenderer *ren, vtkVolume *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(); } int scalarType = this->GetInput()->GetPointData()->GetScalars()->GetDataType(); if (scalarType != VTK_UNSIGNED_SHORT && scalarType != VTK_UNSIGNED_CHAR) { vtkErrorMacro ("Cannot volume render data of type " << vtkImageScalarTypeNameMacro(scalarType) << ", only unsigned char or unsigned short."); return; } // Start timing now. We didn't want to capture the update of the // input data in the times this->Timer->StartTimer(); this->ConvertCroppingRegionPlanesToVoxels(); this->UpdateShadingTables( ren, vol ); // This is the input of this mapper vtkImageData *input = this->GetInput(); // Get the camera from the renderer vtkCamera *cam = ren->GetActiveCamera(); // Get the aspect ratio from the renderer. This is needed for the // computation of the perspective matrix ren->ComputeAspect(); double *aspect = ren->GetAspect(); // 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->Identity(); this->PerspectiveTransform->Concatenate( cam->GetProjectionTransformMatrix(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->ComputeMatrices( input, vol ); // 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 double 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 *= sqrt(oldTime / newTime); this->ImageSampleDistance = (this->ImageSampleDistance>this->MaximumImageSampleDistance)? (this->MaximumImageSampleDistance):(this->ImageSampleDistance); this->ImageSampleDistance = (this->ImageSampleDistanceMinimumImageSampleDistance)? (this->MinimumImageSampleDistance):(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 int width, height; ren->GetTiledSize(&width, &height); this->ImageViewportSize[0] = static_cast(width/this->ImageSampleDistance); this->ImageViewportSize[1] = static_cast(height/this->ImageSampleDistance); // Compute row bounds. This will also compute the size of the image to // render, allocate the space if necessary, and clear the image where // required if ( this->ComputeRowBounds( vol, ren ) ) { vtkVolumeRayCastStaticInfo *staticInfo = new vtkVolumeRayCastStaticInfo; staticInfo->ClippingPlane = NULL; staticInfo->Volume = vol; staticInfo->Renderer = ren; staticInfo->ScalarDataPointer = this->GetInput()->GetPointData()->GetScalars()->GetVoidPointer(0); staticInfo->ScalarDataType = this->GetInput()->GetPointData()->GetScalars()->GetDataType(); // Do we need to capture the z buffer to intermix intersecting // geometry? If so, do it here if ( this->IntermixIntersectingGeometry && ren->GetNumberOfPropsRendered() ) { int x1, x2, y1, y2; // turn this->ImageOrigin into (x1,y1) in window (not viewport!) // coordinates. x1 = static_cast ( viewport[0] * static_cast(renWinSize[0]) + static_cast(this->ImageOrigin[0]) * this->ImageSampleDistance ); y1 = static_cast ( viewport[1] * static_cast(renWinSize[1]) + static_cast(this->ImageOrigin[1]) * this->ImageSampleDistance); // compute z buffer size this->ZBufferSize[0] = static_cast( static_cast(this->ImageInUseSize[0]) * this->ImageSampleDistance); this->ZBufferSize[1] = static_cast( static_cast(this->ImageInUseSize[1]) * this->ImageSampleDistance); // Use the size to compute (x2,y2) in window coordinates x2 = x1 + this->ZBufferSize[0] - 1; y2 = y1 + this->ZBufferSize[1] - 1; // This is the z buffer origin (in viewport coordinates) this->ZBufferOrigin[0] = static_cast( static_cast(this->ImageOrigin[0]) * this->ImageSampleDistance); this->ZBufferOrigin[1] = static_cast( static_cast(this->ImageOrigin[1]) * this->ImageSampleDistance); // Capture the z buffer this->ZBuffer = ren->GetRenderWindow()->GetZbufferData(x1,y1,x2,y2); } // This must be done before FunctionInitialize since FunctionInitialize // depends on the gradient opacity constant (computed in here) to // determine whether to save the gradient magnitudes vol->UpdateTransferFunctions( ren ); // Requires UpdateTransferFunctions to have been called first this->VolumeRayCastFunction->FunctionInitialize( ren, vol, staticInfo ); vol->UpdateScalarOpacityforSampleSize( ren, this->SampleDistance ); staticInfo->CameraThickness = static_cast(ren->GetActiveCamera()->GetThickness()); // Copy the viewToVoxels matrix to 16 floats int i, j; for ( j = 0; j < 4; j++ ) { for ( i = 0; i < 4; i++ ) { staticInfo->ViewToVoxelsMatrix[j*4+i] = static_cast(this->ViewToVoxelsMatrix->GetElement(j,i)); } } // Copy the worldToVoxels matrix to 16 floats for ( j = 0; j < 4; j++ ) { for ( i = 0; i < 4; i++ ) { staticInfo->WorldToVoxelsMatrix[j*4+i] = static_cast(this->WorldToVoxelsMatrix->GetElement(j,i)); } } // Copy the voxelsToWorld matrix to 16 floats for ( j = 0; j < 4; j++ ) { for ( i = 0; i < 4; i++ ) { staticInfo->VoxelsToWorldMatrix[j*4+i] = static_cast(this->VoxelsToWorldMatrix->GetElement(j,i)); } } if ( this->ClippingPlanes ) { this->InitializeClippingPlanes( staticInfo, this->ClippingPlanes ); } else { staticInfo->NumberOfClippingPlanes = 0; } // Copy in the image info staticInfo->ImageInUseSize[0] = this->ImageInUseSize[0]; staticInfo->ImageInUseSize[1] = this->ImageInUseSize[1]; staticInfo->ImageMemorySize[0] = this->ImageMemorySize[0]; staticInfo->ImageMemorySize[1] = this->ImageMemorySize[1]; staticInfo->ImageViewportSize[0] = this->ImageViewportSize[0]; staticInfo->ImageViewportSize[1] = this->ImageViewportSize[1]; staticInfo->ImageOrigin[0] = this->ImageOrigin[0]; staticInfo->ImageOrigin[1] = this->ImageOrigin[1]; staticInfo->Image = this->Image; staticInfo->RowBounds = this->RowBounds; // Set the number of threads to use for ray casting, // then set the execution method and do it. this->Threader->SetSingleMethod( VolumeRayCastMapper_CastRays, static_cast(staticInfo)); this->Threader->SingleMethodExecute(); if ( !ren->GetRenderWindow()->GetAbortRender() ) { float depth; if ( this->IntermixIntersectingGeometry ) { depth = this->MinimumViewDistance; } else { depth = -1; } this->ImageDisplayHelper-> RenderTexture( vol, ren, this->ImageMemorySize, this->ImageViewportSize, this->ImageInUseSize, this->ImageOrigin, depth, this->Image ); 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 ( staticInfo->ClippingPlane ) { delete [] staticInfo->ClippingPlane; } delete staticInfo; if ( this->ZBuffer ) { delete [] this->ZBuffer; this->ZBuffer = NULL; } } } VTK_THREAD_RETURN_TYPE VolumeRayCastMapper_CastRays( void *arg ) { // Get the info out of the input structure int threadID = ((vtkMultiThreader::ThreadInfo *)(arg))->ThreadID; int threadCount = ((vtkMultiThreader::ThreadInfo *)(arg))->NumberOfThreads; vtkVolumeRayCastStaticInfo *staticInfo = (vtkVolumeRayCastStaticInfo *)((vtkMultiThreader::ThreadInfo *)arg)->UserData; int i, j, k; unsigned char *ucptr; vtkVolumeRayCastMapper *me = vtkVolumeRayCastMapper::SafeDownCast( staticInfo->Volume->GetMapper() ); if ( !me ) { vtkGenericWarningMacro("The volume does not have a ray cast mapper!"); return VTK_THREAD_RETURN_VALUE; } vtkVolumeRayCastDynamicInfo *dynamicInfo = new vtkVolumeRayCastDynamicInfo; // Initialize this to avoid purify problems dynamicInfo->ScalarValue = 0; float *rayStart = dynamicInfo->TransformedStart; float *rayEnd = dynamicInfo->TransformedEnd; float *rayDirection = dynamicInfo->TransformedDirection; float *rayStep = dynamicInfo->TransformedIncrement; float norm; float viewRay[3]; float rayCenter[3]; float absStep[3]; float voxelPoint[4]; // We need to know what the center ray is (in voxel coordinates) to // compute sampling distance later on. Save it in an instance variable. // This is tbe near end of the center ray in view coordinates // Convert it to voxel coordinates viewRay[0] = viewRay[1] = viewRay[2] = 0.0; vtkVRCMultiplyPointMacro( viewRay, rayStart, staticInfo->ViewToVoxelsMatrix ); // This is the far end of the center ray in view coordinates // Convert it to voxel coordiantes viewRay[2] = 1.0; vtkVRCMultiplyPointMacro( viewRay, voxelPoint, staticInfo->ViewToVoxelsMatrix ); // Turn these two points into a vector rayCenter[0] = voxelPoint[0] - rayStart[0]; rayCenter[1] = voxelPoint[1] - rayStart[1]; rayCenter[2] = voxelPoint[2] - rayStart[2]; // normalize the vector based on world coordinate distance // This way we can scale by sample distance and it will work // out even though we are in voxel coordinates rayCenter[0] /= staticInfo->CameraThickness; rayCenter[1] /= staticInfo->CameraThickness; rayCenter[2] /= staticInfo->CameraThickness; float centerScale = sqrt( (rayCenter[0] * rayCenter[0]) + (rayCenter[1] * rayCenter[1]) + (rayCenter[2] * rayCenter[2]) ); rayCenter[0] /= centerScale; rayCenter[1] /= centerScale; rayCenter[2] /= centerScale; float bounds[6]; int dim[3]; me->GetInput()->GetDimensions(dim); bounds[0] = bounds[2] = bounds[4] = 0.0; bounds[1] = dim[0]-1; bounds[3] = dim[1]-1; bounds[5] = dim[2]-1; // If we have a simple crop box then we can tighten the bounds if ( me->Cropping && me->CroppingRegionFlags == 0x2000 ) { bounds[0] = me->VoxelCroppingRegionPlanes[0]; bounds[1] = me->VoxelCroppingRegionPlanes[1]; bounds[2] = me->VoxelCroppingRegionPlanes[2]; bounds[3] = me->VoxelCroppingRegionPlanes[3]; bounds[4] = me->VoxelCroppingRegionPlanes[4]; bounds[5] = me->VoxelCroppingRegionPlanes[5]; } bounds[0] = (bounds[0] < 0)?(0):(bounds[0]); bounds[0] = (bounds[0] > dim[0]-1)?(dim[0]-1):(bounds[0]); bounds[1] = (bounds[1] < 0)?(0):(bounds[1]); bounds[1] = (bounds[1] > dim[0]-1)?(dim[0]-1):(bounds[1]); bounds[2] = (bounds[2] < 0)?(0):(bounds[2]); bounds[2] = (bounds[2] > dim[1]-1)?(dim[1]-1):(bounds[2]); bounds[3] = (bounds[3] < 0)?(0):(bounds[3]); bounds[3] = (bounds[3] > dim[1]-1)?(dim[1]-1):(bounds[3]); bounds[4] = (bounds[4] < 0)?(0):(bounds[4]); bounds[4] = (bounds[4] > dim[2]-1)?(dim[2]-1):(bounds[4]); bounds[5] = (bounds[5] < 0)?(0):(bounds[5]); bounds[5] = (bounds[5] > dim[2]-1)?(dim[2]-1):(bounds[5]); // The rounding tie-breaker is used to protect against a very small rounding error // introduced in the QuickFloor function, which is used by the vtkFloorFuncMacro macro. bounds[1] -= vtkFastNumericConversion::RoundingTieBreaker(); bounds[3] -= vtkFastNumericConversion::RoundingTieBreaker(); bounds[5] -= vtkFastNumericConversion::RoundingTieBreaker(); int *imageInUseSize = staticInfo->ImageInUseSize; int *imageMemorySize = staticInfo->ImageMemorySize; int *imageViewportSize = staticInfo->ImageViewportSize; int *imageOrigin = staticInfo->ImageOrigin; int *rowBounds = staticInfo->RowBounds; unsigned char *imagePtr = staticInfo->Image; float sampleDistance = me->GetSampleDistance(); float val; vtkRenderWindow *renWin = staticInfo->Renderer->GetRenderWindow(); // Compute the offset valuex for viewing rays - this is the 1 / fullSize // value to add to the computed location so that they falls between // -1 + 1/fullSize and 1 - 1/fullSize and are each 2/fullSize apart. // fullSize is the viewport size along the corresponding direction (in // pixels) float offsetX = 1.0 / static_cast(imageViewportSize[0]); float offsetY = 1.0 / static_cast(imageViewportSize[1]); // Some variables needed for non-subvolume cropping float fullRayStart[3]; float fullRayEnd[3]; float fullRayDirection[3]; int bitLoop, bitFlag; float rgbaArray[40], distanceArray[10], scalarArray[10]; float tmp, tmpArray[4]; int arrayCount; for ( j = 0; j < imageInUseSize[1]; j++ ) { if ( j%threadCount != threadID ) { continue; } if ( !threadID ) { if ( renWin->CheckAbortStatus() ) { break; } } else if ( renWin->GetAbortRender() ) { break; } ucptr = imagePtr + 4*(j*imageMemorySize[0] + rowBounds[j*2]); // compute the view point y value for this row. Do this by // taking our pixel position, adding the image origin then dividing // by the full image size to get a number from 0 to 1-1/fullSize. Then, // multiply by two and subtract one to get a number from // -1 to 1 - 2/fullSize. Then ass offsetX (which is 1/fullSize) to // center it. viewRay[1] = ((static_cast(j) + static_cast(imageOrigin[1])) / imageViewportSize[1]) * 2.0 - 1.0 + offsetY; for ( i = rowBounds[j*2]; i <= rowBounds[j*2+1]; i++ ) { // Initialize for the cases where the ray doesn't intersect anything ucptr[0] = 0; ucptr[1] = 0; ucptr[2] = 0; ucptr[3] = 0; // Convert the view coordinates for the start and end of the // ray into voxel coordinates, then compute the origin and direction. // compute the view point x value for this pixel. Do this by // taking our pixel position, adding the image origin then dividing // by the full image size to get a number from 0 to 1-1/fullSize. Then, // multiply by two and subtract one to get a number from // -1 to 1 - 2/fullSize. Then ass offsetX (which is 1/fullSize) to // center it. viewRay[0] = ((static_cast(i) + static_cast(imageOrigin[0])) / imageViewportSize[0]) * 2.0 - 1.0 + offsetX; // Now transform this point with a z value of 0 for the ray start, and // a z value of 1 for the ray end. This corresponds to the near and far // plane locations. If IntermixIntersectingGeometry is on, then use // the zbuffer value instead of 1.0 viewRay[2] = 0.0; vtkVRCMultiplyPointMacro( viewRay, rayStart, staticInfo->ViewToVoxelsMatrix ); viewRay[2] = (me->ZBuffer)?(me->GetZBufferValue(i,j)):(1.0); vtkVRCMultiplyPointMacro( viewRay, rayEnd, staticInfo->ViewToVoxelsMatrix ); rayDirection[0] = rayEnd[0] - rayStart[0]; rayDirection[1] = rayEnd[1] - rayStart[1]; rayDirection[2] = rayEnd[2] - rayStart[2]; // If cropping is off, or we are just doing a subvolume, we can // do the easy thing here if ( !me->Cropping || me->CroppingRegionFlags == 0x2000 ) { if ( me->ClipRayAgainstVolume( dynamicInfo, bounds ) && ( staticInfo->NumberOfClippingPlanes == 0 || me->ClipRayAgainstClippingPlanes( dynamicInfo, staticInfo ) ) ) { rayDirection[0] = rayEnd[0] - rayStart[0]; rayDirection[1] = rayEnd[1] - rayStart[1]; rayDirection[2] = rayEnd[2] - rayStart[2]; // Find the length of the input ray. It is not normalized norm = sqrt( rayDirection[0] * rayDirection[0] + rayDirection[1] * rayDirection[1] + rayDirection[2] * rayDirection[2] ); // Normalize this ray into rayStep rayStep[0] = rayDirection[0] / norm; rayStep[1] = rayDirection[1] / norm; rayStep[2] = rayDirection[2] / norm; // Correct for perspective in the sample distance. 1.0 over the // dot product between this ray and the center ray is the // correction factor to allow samples to be taken on parallel // planes rather than concentric hemispheres. This factor will // compute to be 1.0 for parallel. val = ( rayStep[0] * rayCenter[0] + rayStep[1] * rayCenter[1] + rayStep[2] * rayCenter[2] ); norm = (val != 0)?(1.0/val):(1.0); // Now multiple the normalized step by the sample distance and this // correction factor to find the actual step rayStep[0] *= norm * sampleDistance * centerScale; rayStep[1] *= norm * sampleDistance * centerScale; rayStep[2] *= norm * sampleDistance * centerScale; absStep[0] = ( rayStep[0] < 0.0 )?(-rayStep[0]):(rayStep[0]); absStep[1] = ( rayStep[1] < 0.0 )?(-rayStep[1]):(rayStep[1]); absStep[2] = ( rayStep[2] < 0.0 )?(-rayStep[2]):(rayStep[2]); if ( absStep[0] >= absStep[1] && absStep[0] >= absStep[2] ) { dynamicInfo->NumberOfStepsToTake = static_cast ((rayEnd[0]-rayStart[0]) / rayStep[0]); } else if ( absStep[1] >= absStep[2] && absStep[1] >= absStep[0] ) { dynamicInfo->NumberOfStepsToTake = static_cast ((rayEnd[1]-rayStart[1]) / rayStep[1]); } else { dynamicInfo->NumberOfStepsToTake = static_cast ((rayEnd[2]-rayStart[2]) / rayStep[2]); } me->VolumeRayCastFunction->CastRay( dynamicInfo, staticInfo ); if ( dynamicInfo->Color[3] > 0.0 ) { val = dynamicInfo->Color[0]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[0] = static_cast(val); val = dynamicInfo->Color[1]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[1] = static_cast(val); val = dynamicInfo->Color[2]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[2] = static_cast(val); val = dynamicInfo->Color[3]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[3] = static_cast(val); } } } // Otherwise, cropping is on and we don't have a simple subvolume. // We'll have to cast a ray for each of the 27 regions that is on // and composite the results. else { // We'll keep an array of regions that we intersect, the arrayCount // variable will count how many of them we have arrayCount = 0; // Save the ray start, end, and direction. We will modify this // during each iteration of the loop for the current cropping region. fullRayStart[0] = rayStart[0]; fullRayStart[1] = rayStart[1]; fullRayStart[2] = rayStart[2]; fullRayEnd[0] = rayEnd[0]; fullRayEnd[1] = rayEnd[1]; fullRayEnd[2] = rayEnd[2]; fullRayDirection[0] = rayDirection[0]; fullRayDirection[1] = rayDirection[1]; fullRayDirection[2] = rayDirection[2]; // Loop through the twenty seven cropping regions bitFlag = 1; for ( bitLoop = 0; bitLoop < 27; bitLoop++ ) { // Check if this cropping region is on if ( !(me->CroppingRegionFlags & bitFlag) ) { bitFlag = bitFlag << 1; continue; } // Restore the ray information rayStart[0] = fullRayStart[0]; rayStart[1] = fullRayStart[1]; rayStart[2] = fullRayStart[2]; rayEnd[0] = fullRayEnd[0]; rayEnd[1] = fullRayEnd[1]; rayEnd[2] = fullRayEnd[2]; rayDirection[0] = fullRayDirection[0]; rayDirection[1] = fullRayDirection[1]; rayDirection[2] = fullRayDirection[2]; // Figure out the bounds of the cropping region // along the X axis switch ( bitLoop % 3 ) { case 0: bounds[0] = 0; bounds[1] = me->VoxelCroppingRegionPlanes[0]; break; case 1: bounds[0] = me->VoxelCroppingRegionPlanes[0]; bounds[1] = me->VoxelCroppingRegionPlanes[1]; break; case 2: bounds[0] = me->VoxelCroppingRegionPlanes[1]; bounds[1] = (staticInfo->DataSize[0] - 1) - vtkFastNumericConversion::RoundingTieBreaker(); break; } // Figure out the bounds of the cropping region // along the Y axis switch ( (bitLoop % 9) / 3 ) { case 0: bounds[2] = 0; bounds[3] = me->VoxelCroppingRegionPlanes[2]; break; case 1: bounds[2] = me->VoxelCroppingRegionPlanes[2]; bounds[3] = me->VoxelCroppingRegionPlanes[3]; break; case 2: bounds[2] = me->VoxelCroppingRegionPlanes[3]; bounds[3] = (staticInfo->DataSize[1] - 1) - vtkFastNumericConversion::RoundingTieBreaker(); break; } // Figure out the bounds of the cropping region // along the Z axis switch ( bitLoop / 9 ) { case 0: bounds[4] = 0; bounds[5] = me->VoxelCroppingRegionPlanes[4]; break; case 1: bounds[4] = me->VoxelCroppingRegionPlanes[4]; bounds[5] = me->VoxelCroppingRegionPlanes[5]; break; case 2: bounds[4] = me->VoxelCroppingRegionPlanes[5]; bounds[5] = (staticInfo->DataSize[2] - 1) - vtkFastNumericConversion::RoundingTieBreaker(); break; } // Check against the bounds of the volume for ( k = 0; k < 3; k++ ) { if ( bounds[2*k] < 0 ) { bounds[2*k] = 0; } if ( bounds[2*k + 1] > ((staticInfo->DataSize[k]-1) - vtkFastNumericConversion::RoundingTieBreaker()) ) { bounds[2*k + 1] = ((staticInfo->DataSize[k] - 1) - vtkFastNumericConversion::RoundingTieBreaker()); } } // Clip against the volume and the clipping planes if ( me->ClipRayAgainstVolume( dynamicInfo, bounds ) && ( staticInfo->NumberOfClippingPlanes == 0 || me->ClipRayAgainstClippingPlanes( dynamicInfo, staticInfo ) ) ) { // The ray start and end may have changed - recompute // the direction rayDirection[0] = rayEnd[0] - rayStart[0]; rayDirection[1] = rayEnd[1] - rayStart[1]; rayDirection[2] = rayEnd[2] - rayStart[2]; // Find the length of the ray. It is not normalized yet norm = sqrt( rayDirection[0] * rayDirection[0] + rayDirection[1] * rayDirection[1] + rayDirection[2] * rayDirection[2] ); // Normalize this ray into rayStep rayStep[0] = rayDirection[0] / norm; rayStep[1] = rayDirection[1] / norm; rayStep[2] = rayDirection[2] / norm; // Correct for perspective in the sample distance. 1.0 over the // dot product between this ray and the center ray is the // correction factor to allow samples to be taken on parallel // planes rather than concentric hemispheres. This factor will // compute to be 1.0 for parallel. val = ( rayStep[0] * rayCenter[0] + rayStep[1] * rayCenter[1] + rayStep[2] * rayCenter[2] ); norm = (val != 0)?(1.0/val):(1.0); // Now multiple the normalized step by the sample distance and this // correction factor to find the actual step rayStep[0] *= norm * sampleDistance * centerScale; rayStep[1] *= norm * sampleDistance * centerScale; rayStep[2] *= norm * sampleDistance * centerScale; // Find the major direction to determine the number of // steps to take absStep[0] = ( rayStep[0] < 0.0 )?(-rayStep[0]):(rayStep[0]); absStep[1] = ( rayStep[1] < 0.0 )?(-rayStep[1]):(rayStep[1]); absStep[2] = ( rayStep[2] < 0.0 )?(-rayStep[2]):(rayStep[2]); if ( absStep[0] >= absStep[1] && absStep[0] >= absStep[2] ) { dynamicInfo->NumberOfStepsToTake = static_cast ((rayEnd[0]-rayStart[0]) / rayStep[0]); } else if ( absStep[1] >= absStep[2] && absStep[1] >= absStep[0] ) { dynamicInfo->NumberOfStepsToTake = static_cast ((rayEnd[1]-rayStart[1]) / rayStep[1]); } else { dynamicInfo->NumberOfStepsToTake = static_cast ((rayEnd[2]-rayStart[2]) / rayStep[2]); } // Cast the ray me->VolumeRayCastFunction->CastRay( dynamicInfo, staticInfo ); // If the ray returns a non-transparent color, store this // in our arrays of distances and colors if ( dynamicInfo->Color[3] > 0.0 ) { // Figure out the distance from this ray start to the full ray // start and use this to sort the ray segments for ( k = 0; k < 3; k++ ) { if ( absStep[k] >= absStep[(k+1)%3] && absStep[k] >= absStep[(k+2)%3] ) { distanceArray[arrayCount] = (rayStart[k] - fullRayStart[k]) / rayStep[k]; break; } } // Store the ray color rgbaArray[4*arrayCount ] = dynamicInfo->Color[0]; rgbaArray[4*arrayCount+1] = dynamicInfo->Color[1]; rgbaArray[4*arrayCount+2] = dynamicInfo->Color[2]; rgbaArray[4*arrayCount+3] = dynamicInfo->Color[3]; scalarArray[arrayCount] = dynamicInfo->ScalarValue; if ( !staticInfo->MIPFunction ) { // Do a sort pass (one iteration of bubble sort each time an // element is added. The array stores element from farthest // to closest for ( k = arrayCount; k > 0 && distanceArray[k] > distanceArray[k-1]; k-- ) { tmp = distanceArray[k]; distanceArray[k] = distanceArray[k-1]; distanceArray[k-1] = tmp; tmpArray[0] = rgbaArray[4*k ]; tmpArray[1] = rgbaArray[4*k+1]; tmpArray[2] = rgbaArray[4*k+2]; tmpArray[3] = rgbaArray[4*k+3]; rgbaArray[4*k ] = rgbaArray[4*(k-1) ]; rgbaArray[4*k+1] = rgbaArray[4*(k-1)+1]; rgbaArray[4*k+2] = rgbaArray[4*(k-1)+2]; rgbaArray[4*k+3] = rgbaArray[4*(k-1)+3]; rgbaArray[4*(k-1) ] = tmpArray[0]; rgbaArray[4*(k-1)+1] = tmpArray[1]; rgbaArray[4*(k-1)+2] = tmpArray[2]; rgbaArray[4*(k-1)+3] = tmpArray[3]; } } arrayCount++; } } // Move the bit over by one bitFlag = bitFlag << 1; } // We have encountered something in at least one crop region - // merge all results into one RGBA value if ( arrayCount ) { // Is this MIP compositing? We need to treat this differently if ( staticInfo->MIPFunction ) { dynamicInfo->Color[0] = 0.0; dynamicInfo->Color[1] = 0.0; dynamicInfo->Color[2] = 0.0; dynamicInfo->Color[3] = 0.0; dynamicInfo->ScalarValue = 0.0; // If we are maximizing the opacity, find the max Color[3] if ( staticInfo->MaximizeOpacity ) { for ( k = 0; k < arrayCount; k++ ) { if ( rgbaArray[k*4+3] > dynamicInfo->Color[3] ) { dynamicInfo->Color[0] = rgbaArray[k*4 ]; dynamicInfo->Color[1] = rgbaArray[k*4+1]; dynamicInfo->Color[2] = rgbaArray[k*4+2]; dynamicInfo->Color[3] = rgbaArray[k*4+3]; } } } // Otherwise we are maximizing scalar value else { for ( k = 0; k < arrayCount; k++ ) { if ( scalarArray[k] > dynamicInfo->ScalarValue ) { dynamicInfo->Color[0] = rgbaArray[k*4 ]; dynamicInfo->Color[1] = rgbaArray[k*4+1]; dynamicInfo->Color[2] = rgbaArray[k*4+2]; dynamicInfo->Color[3] = rgbaArray[k*4+3]; dynamicInfo->ScalarValue = scalarArray[k]; } } } } else { // Now we have the sorted distances / colors, put them together // in a back-to-front order. First, initialize the color to black // and the remaining opacity (color[3]) to 1.0 dynamicInfo->Color[0] = 0.0; dynamicInfo->Color[1] = 0.0; dynamicInfo->Color[2] = 0.0; dynamicInfo->Color[3] = 1.0; // Now do alpha blending, keeping remaining opacity in color[3] for ( k = 0; k < arrayCount; k++ ) { dynamicInfo->Color[0] = dynamicInfo->Color[0] * (1.0 - rgbaArray[k*4 + 3]) + rgbaArray[k*4 + 0]; dynamicInfo->Color[1] = dynamicInfo->Color[1] * (1.0 - rgbaArray[k*4 + 3]) + rgbaArray[k*4 + 1]; dynamicInfo->Color[2] = dynamicInfo->Color[2] * (1.0 - rgbaArray[k*4 + 3]) + rgbaArray[k*4 + 2]; dynamicInfo->Color[3] *= 1.0 - rgbaArray[k*4 + 3]; } // Take 1.0 - color[3] to convert from remaining opacity to alpha. dynamicInfo->Color[3] = 1.0 - dynamicInfo->Color[3]; } val = dynamicInfo->Color[0]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[0] = static_cast(val); val = dynamicInfo->Color[1]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[1] = static_cast(val); val = dynamicInfo->Color[2]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[2] = static_cast(val); val = dynamicInfo->Color[3]*255.0; val = (val > 255.0)?(255.0):(val); val = (val < 0.0)?( 0.0):(val); ucptr[3] = static_cast(val); } } // Increment the image pointer ucptr+=4; } } delete dynamicInfo; return VTK_THREAD_RETURN_VALUE; } double vtkVolumeRayCastMapper::GetZBufferValue(int x, int y) { int xPos, yPos; xPos = static_cast(static_cast(x) * this->ImageSampleDistance); yPos = static_cast(static_cast(y) * this->ImageSampleDistance); xPos = (xPos >= this->ZBufferSize[0])?(this->ZBufferSize[0]-1):(xPos); yPos = (yPos >= this->ZBufferSize[1])?(this->ZBufferSize[1]-1):(yPos); return *(this->ZBuffer + yPos*this->ZBufferSize[0] + xPos); } int vtkVolumeRayCastMapper::ComputeRowBounds(vtkVolume *vol, vtkRenderer *ren) { 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; float bounds[6]; int dim[3]; this->GetInput()->GetDimensions(dim); bounds[0] = bounds[2] = bounds[4] = 0.0; // The rounding tie-breaker is used to protect against a very small rounding error // introduced in the QuickFloor function, which is used by the vtkFloorFuncMacro macro. bounds[1] = static_cast(dim[0]-1) - vtkFastNumericConversion::RoundingTieBreaker(); bounds[3] = static_cast(dim[1]-1) - vtkFastNumericConversion::RoundingTieBreaker(); bounds[5] = static_cast(dim[2]-1) - vtkFastNumericConversion::RoundingTieBreaker(); 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; } // If we have a simple crop box then we can tighten the bounds // See prior explanation of RoundingTieBreaker if ( this->Cropping && this->CroppingRegionFlags == 0x2000 ) { bounds[0] = this->VoxelCroppingRegionPlanes[0]; bounds[1] = this->VoxelCroppingRegionPlanes[1] - vtkFastNumericConversion::RoundingTieBreaker(); bounds[2] = this->VoxelCroppingRegionPlanes[2]; bounds[3] = this->VoxelCroppingRegionPlanes[3] - vtkFastNumericConversion::RoundingTieBreaker(); bounds[4] = this->VoxelCroppingRegionPlanes[4]; bounds[5] = this->VoxelCroppingRegionPlanes[5] - vtkFastNumericConversion::RoundingTieBreaker(); } // 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]; vtkVRCMultiplyPointMacro( 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)); // We have min/max values from -1.0 to 1.0 now - we want to convert // these to pixel locations. Give a couple of pixels of breathing room // on each side if possible minX = ( minX + 1.0 ) * 0.5 * this->ImageViewportSize[0] - 2; minY = ( minY + 1.0 ) * 0.5 * this->ImageViewportSize[1] - 2; maxX = ( maxX + 1.0 ) * 0.5 * this->ImageViewportSize[0] + 2; maxY = ( maxY + 1.0 ) * 0.5 * this->ImageViewportSize[1] + 2; // If we are outside the view frustum return 0 - there is no need // to render anything if ( ( minX < 0 && maxX < 0 ) || ( minY < 0 && maxY < 0 ) || ( minX > this->ImageViewportSize[0]-1 && maxX > this->ImageViewportSize[0]-1 ) || ( minY > this->ImageViewportSize[1]-1 && maxY > this->ImageViewportSize[1]-1 ) ) { return 0; } int oldImageMemorySize[2]; oldImageMemorySize[0] = this->ImageMemorySize[0]; oldImageMemorySize[1] = this->ImageMemorySize[1]; // Swap the row bounds int *tmpptr; tmpptr = this->RowBounds; this->RowBounds = this->OldRowBounds; this->OldRowBounds = tmpptr; // Check the bounds - the volume might project outside of the // viewing box / frustum so clip it if necessary minX = (minX<0)?(0):(minX); minY = (minY<0)?(0):(minY); maxX = (maxX>this->ImageViewportSize[0]-1)? (this->ImageViewportSize[0]-1):(maxX); maxY = (maxY>this->ImageViewportSize[1]-1)? (this->ImageViewportSize[1]-1):(maxY); // Create the new image, and set its size and position this->ImageInUseSize[0] = static_cast(maxX - minX + 1.0); this->ImageInUseSize[1] = static_cast(maxY - minY + 1.0); // What is a power of 2 size big enough to fit this image? this->ImageMemorySize[0] = 32; this->ImageMemorySize[1] = 32; while ( this->ImageMemorySize[0] < this->ImageInUseSize[0] ) { this->ImageMemorySize[0] *= 2; } while ( this->ImageMemorySize[1] < this->ImageInUseSize[1] ) { this->ImageMemorySize[1] *= 2; } this->ImageOrigin[0] = static_cast(minX); this->ImageOrigin[1] = static_cast(minY); // If the old image size is much too big (more than twice in // either direction) then set the old width to 0 which will // cause the image to be recreated if ( oldImageMemorySize[0] > 2*this->ImageMemorySize[0] || oldImageMemorySize[1] > 2*this->ImageMemorySize[1] ) { oldImageMemorySize[0] = 0; } // If the old image is big enough (but not too big - we handled // that above) then we'll bump up our required size to the // previous one. This will keep us from thrashing. if ( oldImageMemorySize[0] >= this->ImageMemorySize[0] && oldImageMemorySize[1] >= this->ImageMemorySize[1] ) { this->ImageMemorySize[0] = oldImageMemorySize[0]; this->ImageMemorySize[1] = oldImageMemorySize[1]; } // Do we already have a texture big enough? If not, create a new one and // clear it. if ( !this->Image || this->ImageMemorySize[0] > oldImageMemorySize[0] || this->ImageMemorySize[1] > oldImageMemorySize[1] ) { // If there is an image there must be row bounds if ( this->Image ) { delete [] this->Image; delete [] this->RowBounds; delete [] this->OldRowBounds; } this->Image = new unsigned char[(this->ImageMemorySize[0] * this->ImageMemorySize[1] * 4)]; // Create the row bounds array. This will store the start / stop pixel // for each row. This helps eleminate work in areas outside the bounding // hexahedron since a bounding box is not very tight. We keep the old ones // too to help with only clearing where required this->RowBounds = new int [2*this->ImageMemorySize[1]]; this->OldRowBounds = new int [2*this->ImageMemorySize[1]]; for ( i = 0; i < this->ImageMemorySize[1]; i++ ) { this->RowBounds[i*2] = this->ImageMemorySize[0]; this->RowBounds[i*2+1] = -1; this->OldRowBounds[i*2] = this->ImageMemorySize[0]; this->OldRowBounds[i*2+1] = -1; } ucptr = this->Image; for ( i = 0; i < this->ImageMemorySize[0]*this->ImageMemorySize[1]; i++ ) { *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; } } // If we are inside the volume our row bounds indicate every ray must be // cast - we don't need to intersect with the 12 lines if ( insideFlag ) { for ( j = 0; j < this->ImageInUseSize[1]; j++ ) { this->RowBounds[j*2] = 0; this->RowBounds[j*2+1] = this->ImageInUseSize[0] - 1; } } else { // create an array of lines where the y value of the first vertex is less // than or equal to the y value of the second vertex. There are 12 lines, // each containing x1, y1, x2, y2 values. float lines[12][4]; float x1, y1, x2, y2; int xlow, xhigh; int lineIndex[12][2] = {{0,1}, {2,3}, {4,5}, {6,7}, {0,2}, {1,3} ,{4,6}, {5,7}, {0,4}, {1,5}, {2,6}, {3,7}}; for ( i = 0; i < 12; i++ ) { x1 = (viewPoint[lineIndex[i][0]][0]+1.0) * 0.5*this->ImageViewportSize[0] - this->ImageOrigin[0]; y1 = (viewPoint[lineIndex[i][0]][1]+1.0) * 0.5*this->ImageViewportSize[1] - this->ImageOrigin[1]; x2 = (viewPoint[lineIndex[i][1]][0]+1.0) * 0.5*this->ImageViewportSize[0] - this->ImageOrigin[0]; y2 = (viewPoint[lineIndex[i][1]][1]+1.0) * 0.5*this->ImageViewportSize[1] - this->ImageOrigin[1]; if ( y1 < y2 ) { lines[i][0] = x1; lines[i][1] = y1; lines[i][2] = x2; lines[i][3] = y2; } else { lines[i][0] = x2; lines[i][1] = y2; lines[i][2] = x1; lines[i][3] = y1; } } // Now for each row in the image, find out the start / stop pixel // If min > max, then no intersection occurred for ( j = 0; j < this->ImageInUseSize[1]; j++ ) { this->RowBounds[j*2] = this->ImageMemorySize[0]; this->RowBounds[j*2+1] = -1; for ( i = 0; i < 12; i++ ) { if ( j >= lines[i][1] && j <= lines[i][3] && ( lines[i][1] != lines[i][3] ) ) { x1 = lines[i][0] + (static_cast(j) - lines[i][1])/(lines[i][3] - lines[i][1]) * (lines[i][2] - lines[i][0] ); xlow = static_cast(x1 + 1.5); xhigh = static_cast(x1 - 1.0); xlow = (xlow<0)?(0):(xlow); xlow = (xlow>this->ImageInUseSize[0]-1)? (this->ImageInUseSize[0]-1):(xlow); xhigh = (xhigh<0)?(0):(xhigh); xhigh = (xhigh>this->ImageInUseSize[0]-1)?( this->ImageInUseSize[0]-1):(xhigh); if ( xlow < this->RowBounds[j*2] ) { this->RowBounds[j*2] = xlow; } if ( xhigh > this->RowBounds[j*2+1] ) { this->RowBounds[j*2+1] = xhigh; } } } // If they are the same this is either a point on the cube or // all lines were out of bounds (all on one side or the other) // It is safe to ignore the point (since the ray isn't likely // to travel through it enough to actually take a sample) and it // must be ignored in the case where all lines are out of range if ( this->RowBounds[j*2] == this->RowBounds[j*2+1] ) { this->RowBounds[j*2] = this->ImageMemorySize[0]; this->RowBounds[j*2+1] = -1; } } } for ( j = this->ImageInUseSize[1]; j < this->ImageMemorySize[1]; j++ ) { this->RowBounds[j*2] = this->ImageMemorySize[0]; this->RowBounds[j*2+1] = -1; } for ( j = 0; j < this->ImageMemorySize[1]; j++ ) { // New bounds are not overlapping with old bounds - clear between // old bounds only if ( this->RowBounds[j*2+1] < this->OldRowBounds[j*2] || this->RowBounds[j*2] > this->OldRowBounds[j*2+1] ) { ucptr = this->Image + 4*( j*this->ImageMemorySize[0] + this->OldRowBounds[j*2] ); for ( i = 0; i <= (this->OldRowBounds[j*2+1] - this->OldRowBounds[j*2]); i++ ) { *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; } } // New bounds do overlap with old bounds else { // Clear from old min to new min ucptr = this->Image + 4*( j*this->ImageMemorySize[0] + this->OldRowBounds[j*2] ); for ( i = 0; i < (this->RowBounds[j*2] - this->OldRowBounds[j*2]); i++ ) { *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; } // Clear from new max to old max ucptr = this->Image + 4*( j*this->ImageMemorySize[0] + this->RowBounds[j*2+1]+1 ); for ( i = 0; i < (this->OldRowBounds[j*2+1] - this->RowBounds[j*2+1]); i++ ) { *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; *(ucptr++) = 0; } } } return 1; } void vtkVolumeRayCastMapper::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. float volumeOrigin[3]; double *bds = data->GetBounds(); volumeOrigin[0] = bds[0]; volumeOrigin[1] = bds[2]; volumeOrigin[2] = bds[4]; // Get the dimensions of the data. int volumeDimensions[3]; data->GetDimensions( volumeDimensions ); vtkTransform *voxelsTransform = this->VoxelsTransform; vtkTransform *voxelsToViewTransform = this->VoxelsToViewTransform; // 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. this->VolumeMatrix->DeepCopy( vol->GetMatrix() ); voxelsToViewTransform->SetMatrix( VolumeMatrix ); // 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(); } void vtkVolumeRayCastMapper::InitializeClippingPlanes( vtkVolumeRayCastStaticInfo *staticInfo, vtkPlaneCollection *planes ) { vtkPlane *onePlane; double worldNormal[3], worldOrigin[3]; double volumeOrigin[4]; int i; float *worldToVoxelsMatrix; float *voxelsToWorldMatrix; int count; float *clippingPlane; float t; count = planes->GetNumberOfItems(); staticInfo->NumberOfClippingPlanes = count; if ( !count ) { return; } worldToVoxelsMatrix = staticInfo->WorldToVoxelsMatrix; voxelsToWorldMatrix = staticInfo->VoxelsToWorldMatrix; staticInfo->ClippingPlane = new float [4*count]; // loop through all the clipping planes for ( i = 0; i < count; i++ ) { onePlane = static_cast(planes->GetItemAsObject(i)); onePlane->GetNormal(worldNormal); onePlane->GetOrigin(worldOrigin); clippingPlane = staticInfo->ClippingPlane + 4*i; vtkVRCMultiplyNormalMacro( worldNormal, clippingPlane, voxelsToWorldMatrix ); vtkVRCMultiplyPointMacro( worldOrigin, volumeOrigin, worldToVoxelsMatrix ); t = sqrt( clippingPlane[0]*clippingPlane[0] + clippingPlane[1]*clippingPlane[1] + clippingPlane[2]*clippingPlane[2] ); if ( t ) { clippingPlane[0] /= t; clippingPlane[1] /= t; clippingPlane[2] /= t; } clippingPlane[3] = -(clippingPlane[0]*volumeOrigin[0] + clippingPlane[1]*volumeOrigin[1] + clippingPlane[2]*volumeOrigin[2]); } } int vtkVolumeRayCastMapper::ClipRayAgainstClippingPlanes( vtkVolumeRayCastDynamicInfo *dynamicInfo, vtkVolumeRayCastStaticInfo *staticInfo ) { float *clippingPlane; int i; float rayDir[3]; float t, point[3], dp; float *rayStart, *rayEnd; rayStart = dynamicInfo->TransformedStart; rayEnd = dynamicInfo->TransformedEnd; rayDir[0] = rayEnd[0] - rayStart[0]; rayDir[1] = rayEnd[1] - rayStart[1]; rayDir[2] = rayEnd[2] - rayStart[2]; // loop through all the clipping planes for ( i = 0; i < staticInfo->NumberOfClippingPlanes; i++ ) { clippingPlane = staticInfo->ClippingPlane + 4*i; dp = clippingPlane[0]*rayDir[0] + clippingPlane[1]*rayDir[1] + clippingPlane[2]*rayDir[2]; if ( dp != 0.0 ) { t = -( clippingPlane[0]*rayStart[0] + clippingPlane[1]*rayStart[1] + clippingPlane[2]*rayStart[2] + clippingPlane[3]) / dp; if ( t > 0.0 && t < 1.0 ) { point[0] = rayStart[0] + t*rayDir[0]; point[1] = rayStart[1] + t*rayDir[1]; point[2] = rayStart[2] + t*rayDir[2]; if ( dp > 0.0 ) { rayStart[0] = point[0]; rayStart[1] = point[1]; rayStart[2] = point[2]; } else { rayEnd[0] = point[0]; rayEnd[1] = point[1]; rayEnd[2] = point[2]; } rayDir[0] = rayEnd[0] - rayStart[0]; rayDir[1] = rayEnd[1] - rayStart[1]; rayDir[2] = rayEnd[2] - rayStart[2]; } // If the clipping plane is outside the ray segment, then // figure out if that means the ray segment goes to zero (if so // return 0) or doesn't affect it (if so do nothing) else { if ( dp >= 0.0 && t >= 1.0 ) { return 0; } if ( dp <= 0.0 && t <= 0.0 ) { return 0; } } } } return 1; } int vtkVolumeRayCastMapper::ClipRayAgainstVolume( vtkVolumeRayCastDynamicInfo *dynamicInfo, float bounds[6] ) { int loop; float diff; float t; float *rayStart, *rayEnd, *rayDirection; rayStart = dynamicInfo->TransformedStart; rayEnd = dynamicInfo->TransformedEnd; rayDirection = dynamicInfo->TransformedDirection; if ( rayStart[0] >= bounds[1] || rayStart[1] >= bounds[3] || rayStart[2] >= bounds[5] || rayStart[0] < bounds[0] || rayStart[1] < bounds[2] || rayStart[2] < bounds[4] ) { for ( loop = 0; loop < 3; loop++ ) { diff = 0; if ( rayStart[loop] < (bounds[2*loop]+0.01) ) { diff = (bounds[2*loop]+0.01) - rayStart[loop]; } else if ( rayStart[loop] > (bounds[2*loop+1]-0.01) ) { diff = (bounds[2*loop+1]-0.01) - rayStart[loop]; } if ( diff ) { if ( rayDirection[loop] != 0.0 ) { t = diff / rayDirection[loop]; } else { t = -1.0; } if ( t > 0.0 ) { rayStart[0] += rayDirection[0] * t; rayStart[1] += rayDirection[1] * t; rayStart[2] += rayDirection[2] * t; } } } } // If the voxel still isn't inside the volume, then this ray // doesn't really intersect the volume if ( rayStart[0] >= bounds[1] || rayStart[1] >= bounds[3] || rayStart[2] >= bounds[5] || rayStart[0] < bounds[0] || rayStart[1] < bounds[2] || rayStart[2] < bounds[4] ) { return 0; } // The ray does intersect the volume, and we have a starting // position that is inside the volume if ( rayEnd[0] >= bounds[1] || rayEnd[1] >= bounds[3] || rayEnd[2] >= bounds[5] || rayEnd[0] < bounds[0] || rayEnd[1] < bounds[2] || rayEnd[2] < bounds[4] ) { for ( loop = 0; loop < 3; loop++ ) { diff = 0; if ( rayEnd[loop] < (bounds[2*loop]+0.01) ) { diff = (bounds[2*loop]+0.01) - rayEnd[loop]; } else if ( rayEnd[loop] > (bounds[2*loop+1]-0.01) ) { diff = (bounds[2*loop+1]-0.01) - rayEnd[loop]; } if ( diff ) { if ( rayDirection[loop] != 0.0 ) { t = diff / rayDirection[loop]; } else { t = 1.0; } if ( t < 0.0 ) { rayEnd[0] += rayDirection[0] * t; rayEnd[1] += rayDirection[1] * t; rayEnd[2] += rayDirection[2] * t; } } } } // To be absolutely certain our ray remains inside the volume, // recompute the ray direction (since it has changed - it is not // normalized and therefore changes when start/end change) and move // the start/end points in by 1/1000th of the distance. float offset; offset = (rayEnd[0] - rayStart[0])*0.001; rayStart[0] += offset; rayEnd[0] -= offset; offset = (rayEnd[1] - rayStart[1])*0.001; rayStart[1] += offset; rayEnd[1] -= offset; offset = (rayEnd[2] - rayStart[2])*0.001; rayStart[2] += offset; rayEnd[2] -= offset; if ( rayEnd[0] >= bounds[1] || rayEnd[1] >= bounds[3] || rayEnd[2] >= bounds[5] || rayEnd[0] < bounds[0] || rayEnd[1] < bounds[2] || rayEnd[2] < bounds[4] ) { return 0; } return 1; } void vtkVolumeRayCastMapper::UpdateShadingTables( vtkRenderer *ren, vtkVolume *vol ) { int shading; vtkVolumeProperty *volume_property; volume_property = vol->GetProperty(); shading = volume_property->GetShade(); this->GradientEstimator->SetInput( this->GetInput() ); if ( shading ) { this->GradientShader->UpdateShadingTable( ren, vol, this->GradientEstimator ); } } float vtkVolumeRayCastMapper::GetZeroOpacityThreshold( vtkVolume *vol ) { return( this->VolumeRayCastFunction->GetZeroOpacityThreshold( vol ) ); } // Print method for vtkVolumeRayCastMapper void vtkVolumeRayCastMapper::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os,indent); os << indent << "Sample Distance: " << this->SampleDistance << "\n"; os << indent << "Image Sample Distance: " << this->ImageSampleDistance << "\n"; os << indent << "Minimum Image Sample Distance: " << this->MinimumImageSampleDistance << "\n"; os << indent << "Maximum Image Sample Distance: " << this->MaximumImageSampleDistance << "\n"; os << indent << "Auto Adjust Sample Distances: " << this->AutoAdjustSampleDistances << "\n"; os << indent << "Intermix Intersecting Geometry: " << (this->IntermixIntersectingGeometry ? "On\n" : "Off\n"); if ( this->VolumeRayCastFunction ) { os << indent << "Ray Cast Function: " << this->VolumeRayCastFunction<<"\n"; } else { os << indent << "Ray Cast Function: (none)\n"; } if ( this->GradientEstimator ) { os << indent << "Gradient Estimator: " << (this->GradientEstimator) << endl; } else { os << indent << "Gradient Estimator: (none)" << endl; } if ( this->GradientShader ) { os << indent << "Gradient Shader: " << (this->GradientShader) << endl; } else { os << indent << "Gradient Shader: (none)" << endl; } } //---------------------------------------------------------------------------- void vtkVolumeRayCastMapper::ReportReferences(vtkGarbageCollector* collector) { this->Superclass::ReportReferences(collector); // These filters share our input and are therefore involved in a // reference loop. vtkGarbageCollectorReport(collector, this->GradientEstimator, "GradientEstimator"); }