[vtkusers] x-ray effect in volume rendering

baliki2 at freemail.gr baliki2 at freemail.gr
Mon Oct 26 09:59:59 EDT 2009


It works just great! I'll try it in Java, too.

Thanks a lot Mr. Gasteiger!

> Hi "Baliki",
> 
> Below I have inserted a VTK-example program, which reads a raw volume
> dataset and performs normal volume rendering, X-ray and MIP projection. I
> use it for my students to illustrate these different approaches. The samples
> on each ray are interpolated with nearest neighbor (not the best, tri-linear
> is much better but expansive in computation time. It inFollowing steps you
> have to do:
> 1. Define "fname" (in main(...)) with the file name of you dataset
> 2  Specify your raw dataset for the vtkImageReader2 reader (In the example I
> use a volume version(!) of the foot dataset of the VTK distribution)
> 3. Turn on one of the render modes by uncomment one of the
> "currentRenderMode" line at the beginning of the example
> 
> I have included some comments which hopefully helps you to understand each
> step.
> 
> Best regards, Rocco
> 
> //--------------------------------------------------------------------------
> ------------
> //
> // This example reads a volume dataset and then displays it.
> //
> #include "vtkRenderer.h"
> #include "vtkRenderWindow.h"
> #include "vtkRenderWindowInteractor.h"
> #include "vtkCamera.h"
> #include "vtkInteractorStyleSwitch.h"
> #include "vtkImageReader2.h"
> #include "vtkActor.h"
> 
> #include "vtkColorTransferFunction.h"
> #include "vtkPiecewiseFunction.h"
> #include "vtkVolumeProperty.h"
> #include "vtkVolumeRayCastCompositeFunction.h"
> #include "vtkVolumeRayCastMapper.h"
> #include "vtkVolume.h"
> 
> 
> #define mode_MIP						0
> #define mode_XRay						1
> #define mode_Full						2
> 
> static int currentRenderMode = mode_XRay;
> 
> // Necessary for correct MIP and XRay scaling (0-1)
> static float globalMaximumScalarValue = 0;
> ////////////////////////////////////////////////////////////////
> // Definition and implementation of our own ray cast composite function
> ----------------------------------
> class MyVolumeRayCastCompositeFunction : public
> vtkVolumeRayCastCompositeFunction
> {
> 
> public:
> 	static MyVolumeRayCastCompositeFunction *New()
> 	{ return new MyVolumeRayCastCompositeFunction; }
> 
> protected:
> 
> 	////////////////////////////////////////////////////////////////
> 	// Do not change this code
> 	//
> 	// This is called from RenderAnImage (in vtkDepthPARCMapper.cxx)
> 	// It uses the integer data type flag that is passed in
> 	// to determine what type of ray needs to be cast (which is handled
> 	// by a template function.  It also uses the shading and
> 	// interpolation types to determine which template function
> 	// to call.
> 	void MyVolumeRayCastCompositeFunction::CastRay(
> vtkVolumeRayCastDynamicInfo *dynamicInfo,
> 		vtkVolumeRayCastStaticInfo *staticInfo )
> 	{
> 		void *data_ptr;
> 		data_ptr = staticInfo->ScalarDataPointer;
> 
> 		switch ( staticInfo->ScalarDataType ) {
> case VTK_UNSIGNED_CHAR:
> 	vtkCastRay_OwnWork( (unsigned char *)data_ptr, dynamicInfo,
> staticInfo );
> 	break;
> case VTK_UNSIGNED_SHORT:
> 	vtkCastRay_OwnWork( (unsigned short *)data_ptr, dynamicInfo,
> staticInfo );
> 	break;
> default:
> 	vtkWarningMacro ( << "Unsigned char and unsigned short are the only
> supported datatypes for rendering" );
> 	break;
> 		}
> 	}
> 
> 	// This is the template function that actually casts a ray and
> computes
> 	// the composite value.  This version uses nearest neighbor
> interpolation, and performs shading.
> 	template <class T>
> 	void vtkCastRay_OwnWork(T *data_ptr, vtkVolumeRayCastDynamicInfo
> *dynamicInfo, vtkVolumeRayCastStaticInfo *staticInfo )
> 	{
> 		unsigned char   *gmptr = NULL;
> 		float           accum_red_intensity, accum_green_intensity,
> accum_blue_intensity, accum_opacity;
> 		float           opacity;
> 		int             xinc, yinc, zinc;
> 		int             voxel[3];
> 		float           ray_position[3];
> 		T               *dptr;
> 		float           *ScalarToOpacityTF;
> 		float           *ScalarToColorTF;
> 		float           *red_d_shade, *green_d_shade, *blue_d_shade;
> 		float           *red_s_shade, *green_s_shade, *blue_s_shade;
> 		unsigned short  *encoded_normals;
> 		int             currentNormalIndex;
> 		float           red_shaded_value, green_shaded_value,
> blue_shaded_value;
> 		int             offset;
> 		int             steps_this_ray;
> 		int             scalar_value;
> 		float           r, g, b;
> 		int             num_steps;
> 		float           *ray_start, *ray_increment;
> 		int			    shading = staticInfo->Shading;
> 
> 
> 		// For each ray some information are fetched which are
> needed for color and shading
> 		// computation.
> 
> 		// Get ray parameters: Number of sample points, ray start
> position, and step size on the ray
> 		num_steps = dynamicInfo->NumberOfStepsToTake;
> 		ray_start = dynamicInfo->TransformedStart;
> 		ray_increment = dynamicInfo->TransformedIncrement;
> 
> 		// Get diffuse shading table pointers. This table is
> preprocessed by VTK and stores
> 		// for different normal directions its diffuse portion of
> the phong illumination term.
> 		red_d_shade   = staticInfo->RedDiffuseShadingTable;
> 		green_d_shade = staticInfo->GreenDiffuseShadingTable;
> 		blue_d_shade  = staticInfo->BlueDiffuseShadingTable;
> 
> 		// Get specular shading table pointers. This table is
> preprocessed by VTK and stores
> 		// for specular normal directions its diffuse portion of the
> phong illumination term.
> 		red_s_shade = staticInfo->RedSpecularShadingTable;
> 		green_s_shade = staticInfo->GreenSpecularShadingTable;
> 		blue_s_shade = staticInfo->BlueSpecularShadingTable;
> 
> 		// Get a pointer to the encoded normals of this volume. This
> table is preprocessed by VTK and stores
> 		// for each voxel its normal.
> 		encoded_normals = staticInfo->EncodedNormals;
> 
> 		// Get the scalar opacity transfer function which maps
> scalar input values to opacities
> 		ScalarToOpacityTF =
> staticInfo->Volume->GetCorrectedScalarOpacityArray();
> 
> 		// Get the color transfer function which maps scalar input
> values to RGB values
> 		ScalarToColorTF =  staticInfo->Volume->GetRGBArray();
> 
> 		// store increments (distance between sample points in x, y,
> and z direction) in local variables
> 		xinc = staticInfo->DataIncrement[0];
> 		yinc = staticInfo->DataIncrement[1];
> 		zinc = staticInfo->DataIncrement[2];
> 		//cerr<<xinc<<", "<<yinc<<", "<<zinc<<", "<<endl; //constant
> for a data set
> 
> 		// Initialize the ray position (in voxel space)
> 		ray_position[0] = ray_start[0];
> 		ray_position[1] = ray_start[1];
> 		ray_position[2] = ray_start[2];
> 
> 
> 		// Nearest neighbor interpolation - get nearest voxel (in
> voxel space) - you have to change this for tri-linear interpolation
> 		voxel[0] = vtkRoundFuncMacro( ray_position[0] );
> 		voxel[1] = vtkRoundFuncMacro( ray_position[1] );
> 		voxel[2] = vtkRoundFuncMacro( ray_position[2] );
> 
> 		// So far we haven't accumulated anything
> 		accum_red_intensity   = 0.0;
> 		accum_green_intensity = 0.0;
> 		accum_blue_intensity  = 0.0;
> 		// accumulated opacity
> 		accum_opacity         = 0.0;
> 
> 		// added for MIP mode rendering
> 		float maxOpacity      = 0.0;
> 		float maxValue        = 0.0;
> 
> 		currentNormalIndex    = 0;
> 
> 		// For each step along the ray do something (accumulate,
> find maximum opacity, ...)
> 		// "accum_opacity < 1.0" means early-ray-termination
> 		for ( steps_this_ray = 0; steps_this_ray < num_steps &&
> accum_opacity < 1.0; steps_this_ray++ ) {
> 
> 			// get offset position of current voxel in linear
> list
> 			offset = voxel[2] * zinc + voxel[1] * yinc +
> voxel[0];
> 			// get pointer on current voxel in linear list
> 			dptr = data_ptr + offset;
> 
> 			// get scalar value at current location
> 			scalar_value       = (int) *(dptr);
> 
> 			if (scalar_value > globalMaximumScalarValue)
> 			{
> 				globalMaximumScalarValue = scalar_value;
> 			}
> 
> 			// get opacity value for this scalar_value out of
> the ScalarOpacityTransferFunction
> 			opacity = ScalarToOpacityTF[scalar_value];
> 
> 			// If we have a valid opacity value, then compute
> the shading
> 			if ( opacity ) {
> 				switch (currentRenderMode)
> 				{
> 				case mode_MIP:
> 					// MIP - only maximal value or first
> local maximum above a threshold is visualized
> 					// check if current voxel opacity is
> above max opacity
> 					//if (opacity > maxOpacity)
> 					if (scalar_value > maxValue)
> 					{
> 						accum_red_intensity   =
> opacity * ScalarToColorTF[(scalar_value) * 3    ];
> 						accum_green_intensity =
> opacity * ScalarToColorTF[(scalar_value) * 3 + 1];
> 						accum_blue_intensity  =
> opacity * ScalarToColorTF[(scalar_value) * 3 + 2];
> 
> 						//maxOpacity    = opacity;
> 						maxValue    = scalar_value;
> 						accum_opacity = 1.0-opacity;
> 
> 					}
> 
> 					break;
> 				case mode_XRay:
> 					// X-ray projection mode:
> Accumulation of the scalar values (or opacity values)
> 					accum_red_intensity =
> accum_green_intensity = accum_blue_intensity = accum_red_intensity +
> scalar_value;
> 
> 					break;
> 				case mode_Full:
> 					// get r, g, and b value for current
> scalar value out of the ColorTransferFunction
> 					r = ScalarToColorTF[(scalar_value) *
> 3    ];
> 					g = ScalarToColorTF[(scalar_value) *
> 3 + 1];
> 					b = ScalarToColorTF[(scalar_value) *
> 3 + 2];
> 
> 					red_shaded_value   = opacity * r;
> 					green_shaded_value = opacity * g;
> 					blue_shaded_value  = opacity * b;
> 
> 
> 					// Accumulate color
> 					// (the color of the current value
> is weighted by the remaining opacity of the voxels
> 					// in front of the current one)
> 					accum_red_intensity   =
> red_shaded_value   * (1.0-accum_opacity) + accum_red_intensity;
> 					accum_green_intensity =
> green_shaded_value * (1.0-accum_opacity) + accum_green_intensity;
> 					accum_blue_intensity  =
> blue_shaded_value  * (1.0-accum_opacity) + accum_blue_intensity;
> 
> 					// Accumulate opacity
> 					accum_opacity  =
> opacity*(1.0-accum_opacity) + accum_opacity;
> 
> 					break;
> 				default:
> 					accum_red_intensity = 0.0;
> 					accum_green_intensity = 0.0;
> 					accum_blue_intensity = 0.0;
> 					accum_opacity = 1.0;
> 
> 					break;
> 				}
> 
> 			}
> 
> 			// Increment our position and compute new voxel
> location
> 			ray_position[0] += ray_increment[0];
> 			ray_position[1] += ray_increment[1];
> 			ray_position[2] += ray_increment[2];
> 			voxel[0] = vtkRoundFuncMacro( ray_position[0] );
> 			voxel[1] = vtkRoundFuncMacro( ray_position[1] );
> 			voxel[2] = vtkRoundFuncMacro( ray_position[2] );
> 		}// endFor each step along the ray
> ----------------------------------------------------------------------------
> -----------
> 
> 
> 		// Now we stored in accum_{red,green,blue}_intensity the
> accumulated RGB values along a ray
> 		// added for XRay mode rendering
> 		if( currentRenderMode == mode_XRay)
> 		{
> 			// divide accumulated scalar values by the number of
> steps times maximum scalar value to get the average
> 			accum_red_intensity = accum_green_intensity =
> accum_blue_intensity = accum_red_intensity /
> float(num_steps*globalMaximumScalarValue);
> 			accum_opacity = 1.0-accum_red_intensity;
> 
> 		}
> 
> 		// Cap the accumulated intensities at 1.0
> 		if ( accum_red_intensity > 1.0 )	{
> accum_red_intensity = 1.0; }
> 		if ( accum_green_intensity > 1.0 ){ accum_green_intensity =
> 1.0;}
> 		if ( accum_blue_intensity > 1.0 )	{
> accum_blue_intensity = 1.0;}
> 
> 		// Set the return pixel value.
> 		dynamicInfo->Color[0] = accum_red_intensity;
> 		dynamicInfo->Color[1] = accum_green_intensity;
> 		dynamicInfo->Color[2] = accum_blue_intensity;
> 		dynamicInfo->Color[3] = accum_opacity;
> 		dynamicInfo->NumberOfStepsTaken = steps_this_ray;
> 
> 
> 	}
> 
> };
> //
> ----------------------------------------------------------------------------
> -
> 
> 
> ////////////////////////////////////////////////////////////////
> //
> int main (int argc, char **argv)
> {
> 	//check for  arguments (e.g. a filename)
> 	const char* fname = "foot.raw";
> 
> 	// Create the renderer, the render window, and the interactor. The
> renderer
> 	// draws into the render window, the interactor enables mouse- and
> 	// keyboard-based interaction with the data within the render
> window.
> 	vtkRenderer *renderer = vtkRenderer::New();
> 	//renderer->SetBackground( 0.2, 0.2, 0.4 );
> 	renderer->SetBackground( 1.0, 1.0, 1.0 );
> 	vtkCamera *camera = renderer->GetActiveCamera();
> 	camera->ParallelProjectionOff();
> 	camera->SetViewUp (0, 0, -1);
> 	camera->SetPosition (-1, 2, -0.5);
> 	vtkRenderWindow *renderWindow = vtkRenderWindow::New();
> 	renderWindow->SetSize( 640, 480 );     // Set the size of the render
> window (in pixel).
> 	renderWindow->AddRenderer(renderer);
> 	vtkRenderWindowInteractor *interactionRenderer =
> vtkRenderWindowInteractor::New();
> 	interactionRenderer->SetRenderWindow(renderWindow);
> 
> reinterpret_cast<vtkInteractorStyleSwitch*>(interactionRenderer->GetInteract
> orStyle())->SetCurrentStyleToTrackballCamera();
> 
> 	vtkImageReader2 *volumeDataset = vtkImageReader2::New();
> 	volumeDataset->SetFileName(fname);
> 	volumeDataset->SetFileDimensionality(3);
> 	volumeDataset->SetDataExtent(0, 127, 0, 127, 0, 63);
> 	volumeDataset->SetDataSpacing(1.60938, 1.60938, 1.95313);
> 	volumeDataset->SetDataOrigin(0.0, 0.0, 0.0);
> 	volumeDataset->SetDataScalarTypeToUnsignedShort();
> 	volumeDataset->SetDataByteOrderToLittleEndian();
> 	volumeDataset->UpdateWholeExtent();
> 
> 
> 
> ////////////////////////////////////////////////////////////////////////////
> ///////////
> 	// Feel free to change this code
> 	// adjustments to the following code may be required and should be
> explored
> 
> 	// Create transfer function - mapping scalar values [0-...] to COLOR
> [0-1, 0-1, 0-1]
> 	vtkColorTransferFunction *colorTransferFunction =
> vtkColorTransferFunction::New();
> 	colorTransferFunction->AddRGBPoint(   0.0, 0.0,0.0,0.0);
> 	colorTransferFunction->AddRGBPoint( 500.0, 0.9,0.5,0.3);
> 	colorTransferFunction->AddRGBPoint(1100.0, 0.8,0.8,0.6);
> 	colorTransferFunction->AddRGBPoint(1200.0, 0.6,0.6,0.6);
> 
> 	vtkPiecewiseFunction *opacityTransferFunction =
> vtkPiecewiseFunction::New();
> 	opacityTransferFunction->AddPoint(    0, 0.0);
> 	//opacityTransferFunction->AddPoint(  980, 0.1);
> 	//opacityTransferFunction->AddPoint(  1055, 0.2);
> 	opacityTransferFunction->AddPoint(  1200, 0.05);
> 	opacityTransferFunction->AddPoint( 1300, 1.0);
> 
> 	// The property describes how the data will look
> 	vtkVolumeProperty *volumeProperty = vtkVolumeProperty::New();
> 	volumeProperty->SetColor(colorTransferFunction);
> 	volumeProperty->SetScalarOpacity(opacityTransferFunction);
> 	//volumeProperty->ShadeOn(); // request 3d shading (german:
> beleuchtungsberechnung anfordern) //has to be implemented in
> vtkCastRay_OwnWork
> 
> 	//vtkVolumeRayCastCompositeFunction *rayCompositeFunction =
> vtkVolumeRayCastCompositeFunction::New();
> 	MyVolumeRayCastCompositeFunction *rayCompositeFunction =
> MyVolumeRayCastCompositeFunction::New();
> 	// The mapper knows how to render the data
> 	vtkVolumeRayCastMapper *rayCastMapper =
> vtkVolumeRayCastMapper::New();
> 	rayCastMapper->SetInput( volumeDataset->GetOutput());
> 	rayCastMapper->SetVolumeRayCastFunction(rayCompositeFunction);
> 	//Feel free to change parameters for
> 	//rayCastMapper->AutoAdjustSampleDistancesOff();		//
> interactive framerates by lower quality (on/off)
> 	//rayCastMapper->SetSampleDistance(1.0);
> // adjust samples per ray rate if AutoAdjustSampleDistances is OFF
> 	//rayCastMapper->SetImageSampleDistance(1.0);			//
> adjust rays per pixel rate (e.g. 0.5 => 4 rays for each pixel, 2.0 => 1 ray
> for 4 pixel) if AutoAdjustSampleDistances is OFF
> 	//rayCastMapper->SetMaximumImageSampleDistance(6.0);	// max
> sample distance if AutoAdjustSampleDistances is ON
> 	//rayCastMapper->SetMinimumImageSampleDistance(1.0);	// min
> sample distance if AutoAdjustSampleDistances is ON
> 
> 	// The volume holds the mapper and the property and
> 	// can be used to position/orient the volume
> 	vtkVolume *volData= vtkVolume::New();
> 	volData->SetMapper(rayCastMapper);
> 	volData->SetProperty(volumeProperty);
> 
> 	// Actors are added to the renderer.
> 	renderer->AddVolume(volData);
> 
> 
> 	// adjust camera for good initial view
> 	renderer->ResetCamera();
> 	renderer->GetActiveCamera()->Dolly(1);
> 	renderer->ResetCameraClippingRange();
> 
> 	// Start the Interactor
> 	interactionRenderer->Initialize();
> 	interactionRenderer->Start();
> 
> 	// It is important to delete all objects created previously to
> prevent
> 	// memory leaks. In this case, since the program is on its way to
> 	// exiting, it is not so important. But in applications it is
> 	// essential.
> 	renderer->Delete();
> 	renderWindow->Delete();
> 	interactionRenderer->Delete();
> 
> 	return 0;
> }
> 
> 
> -----Original Message-----
> From: baliki2 at freemail.gr [mailto:baliki2 at freemail.gr]
> Sent: Friday, October 23, 2009 2:11 PM
> To: post at rocco-gasteiger.de
> Subject: Re: Re: [vtkusers] x-ray effect in volume rendering
> 
> 
> I don't really understand what you mean by "divide it by
> > numSteps*globalMaximumScalarValue at the end".
> 
> What i would like to get is something like this:
> http://www.onlinetelemedicine.com/HTML/product/sam_images/X-Ray.jpg
> If you look carefully, the bones edges are more bright than the "inside"
> bone area.
> I've achieved a similar effect, but not this exact one:
> http://img12.imageshack.us/i/sendj.png/
> Please download it and zoom in to see it better.
> 
> Thanks!
> 
> > Hello,
> >
> > Yes, you can! ;-)
> >
> > It is quite simply. In your composite function accumulate all sampled
> scalar
> > values (or opacity values) on your ray and divide it by
> > numSteps*globalMaximumScalarValue at the end (to get the average between 0
> > and 1). I've implement it and it should work.
> >
> > Best regards, Rocco
> >
> > --------------------------------------------------
> > Dipl.-Ing. Rocco Gasteiger
> > Otto-von-Guericke University
> > Faculty of Computer Science
> > Department of Simulation and Graphics
> > Universit=E4tsplatz 2, 39106 Magdeburg, Germany
> >
> > Office:  G29-223
> > Phone:   +49 391 67 127 59
> > Fax:=A0=A0   +49 391 67 111 64
> > Website: http://wwwisg.cs.uni-magdeburg.de/cvcms/  =
> >
> > --------------------------------------------------
> >
> >
> >
> > -----Original Message-----
> > From: vtkusers-bounces at vtk.org [mailto:vtkusers-bounces at vtk.org] On Behalf
> > Of baliki2 at freemail.gr
> > Sent: Friday, October 23, 2009 11:30 AM
> > To: vtkusers at vtk.org
> > Subject: [vtkusers] x-ray effect in volume rendering
> >
> >
> > Is there a way that i can achieve an x-ray effect in volume rendering of
> > DICOM images?
> > I want to see the constructed volume as a 2D xray-like "image".
> >
> > I've tied the composite function combined with some tranfer and color
> > functions, but i didn't take the exact x-ray effect. The images i use are
> of
> > short type.
> >
> >
> > _______________________________________________
> > Powered by www.kitware.com
> >
> > Visit other Kitware open-source projects at
> > http://www.kitware.com/opensource/opensource.html
> >
> > Please keep messages on-topic and check the VTK FAQ at:
> > http://www.vtk.org/Wiki/VTK_FAQ
> >
> > Follow this link to subscribe/unsubscribe:
> > http://www.vtk.org/mailman/listinfo/vtkusers
> >
> > _______________________________________________
> > Powered by www.kitware.com
> >
> > Visit other Kitware open-source projects at
> http://www.kitware.com/opensour=
> > ce/opensource.html
> >
> > Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/=
> > Wiki/VTK_FAQ
> >
> > Follow this link to subscribe/unsubscribe:
> > http://www.vtk.org/mailman/listinfo/vtkusers
> 
> _______________________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
> 
> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers





More information about the vtkusers mailing list