From 271df19d73efd8e91ac545bb18fa89ea1b8d30f2 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 21 May 2008 22:16:03 +0000 Subject: [PATCH] - First version of a VTK vtkImageData to vtkPolyData filter: mesh an iso-value in a 3D image. --- .gitignore | 1 + .../demo/Surface_mesher/VTK/CMakeLists.txt | 22 +++ .../VTK/mesh_a_VTK_3D_image.cpp | 142 ++++++++++++++++ .../CGAL/vtkSurfaceMesherContourFilter.h | 156 ++++++++++++++++++ 4 files changed, 321 insertions(+) create mode 100644 Surface_mesher/demo/Surface_mesher/VTK/mesh_a_VTK_3D_image.cpp create mode 100644 Surface_mesher/include/CGAL/vtkSurfaceMesherContourFilter.h diff --git a/.gitignore b/.gitignore index 77572386d32..19f852ed919 100644 --- a/.gitignore +++ b/.gitignore @@ -293,6 +293,7 @@ Surface_mesher/demo/Surface_mesher/.*.deps Surface_mesher/demo/Surface_mesher/VTK/*.exe Surface_mesher/demo/Surface_mesher/VTK/Makefile Surface_mesher/demo/Surface_mesher/VTK/mesh_a_3D_image +Surface_mesher/demo/Surface_mesher/VTK/mesh_a_VTK_3D_image Surface_mesher/demo/Surface_mesher/dump.ff Surface_mesher/demo/Surface_mesher/implicit_surface_mesher Surface_mesher/demo/Surface_mesher/my_makefile diff --git a/Surface_mesher/demo/Surface_mesher/VTK/CMakeLists.txt b/Surface_mesher/demo/Surface_mesher/VTK/CMakeLists.txt index 7d866bd186a..46f83eb11a6 100644 --- a/Surface_mesher/demo/Surface_mesher/VTK/CMakeLists.txt +++ b/Surface_mesher/demo/Surface_mesher/VTK/CMakeLists.txt @@ -45,6 +45,10 @@ SET (SRCS mesh_a_3D_image.cpp ) +SET (SRCS_VTK + mesh_a_VTK_3D_image.cpp +) + # Use the include path and library for Qt that is used by VTK. INCLUDE_DIRECTORIES( ${QT_INCLUDE_DIR} ${QT_QTGUI_INCLUDE_DIR} ${QT_QTCORE_INCLUDE_DIR}) @@ -62,4 +66,22 @@ TARGET_LINK_LIBRARIES( mesh_a_3D_image CGAL-ImageIO ) +SET (SRCS_VTK + mesh_a_VTK_3D_image.cpp +) + +ADD_EXECUTABLE( mesh_a_VTK_3D_image MACOSX_BUNDLE ${SRCS_VTK}) + +TARGET_LINK_LIBRARIES( mesh_a_VTK_3D_image + QVTK + ${QT_LIBRARIES} + vtkRendering + vtkGraphics + vtkIO + vtkCommon + CGAL + CGAL-ImageIO +) + + diff --git a/Surface_mesher/demo/Surface_mesher/VTK/mesh_a_VTK_3D_image.cpp b/Surface_mesher/demo/Surface_mesher/VTK/mesh_a_VTK_3D_image.cpp new file mode 100644 index 00000000000..c6f1868b3f1 --- /dev/null +++ b/Surface_mesher/demo/Surface_mesher/VTK/mesh_a_VTK_3D_image.cpp @@ -0,0 +1,142 @@ +/*========================================================================= + + Program: Visualization Toolkit + Module: $RCSfile: Medical1.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. + +=========================================================================*/ +// +// This example reads a volume dataset, extracts an isosurface that +// represents the skin and displays it. +// + +#include "vtkRenderer.h" +#include "vtkRenderWindow.h" +#include "vtkRenderWindowInteractor.h" +#include "vtkVolume16Reader.h" +#include "vtkPolyDataMapper.h" +#include "vtkActor.h" +#include "vtkOutlineFilter.h" +#include "vtkCamera.h" +#include "vtkProperty.h" +#include "vtkPolyDataNormals.h" +#include + +int main (int argc, char **argv) +{ + if (argc < 2) + { + cout << "Usage: " << argv[0] << " DATADIR/headsq/quarter" << endl; + return 1; + } + + // 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 *aRenderer = vtkRenderer::New(); + vtkRenderWindow *renWin = vtkRenderWindow::New(); + renWin->AddRenderer(aRenderer); + vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New(); + iren->SetRenderWindow(renWin); + + // The following reader is used to read a series of 2D slices (images) + // that compose the volume. The slice dimensions are set, and the + // pixel spacing. The data Endianness must also be specified. The reader + // usese the FilePrefix in combination with the slice number to construct + // filenames using the format FilePrefix.%d. (In this case the FilePrefix + // is the root name of the file: quarter.) + vtkVolume16Reader *v16 = vtkVolume16Reader::New(); + v16->SetDataDimensions (64,64); + v16->SetImageRange (1,93); + v16->SetDataByteOrderToLittleEndian(); + v16->SetFilePrefix (argv[1]); + v16->SetDataSpacing (3.2, 3.2, 1.5); + + // An isosurface, or contour value of 500 is known to correspond to the + // skin of the patient. Once generated, a vtkPolyDataNormals filter is + // is used to create normals for smooth surface shading during rendering. + vtkCGALSurfaceMesherContourFilter *skinExtractor = vtkCGALSurfaceMesherContourFilter::New(); + skinExtractor->SetInputConnection(v16->GetOutputPort()); + skinExtractor->SetValue(500); + vtkPolyDataNormals *skinNormals = vtkPolyDataNormals::New(); + skinNormals->SetInputConnection(skinExtractor->GetOutputPort()); + skinNormals->SetFeatureAngle(60.0); + vtkPolyDataMapper *skinMapper = vtkPolyDataMapper::New(); + skinMapper->SetInputConnection(skinNormals->GetOutputPort()); + skinMapper->ScalarVisibilityOff(); + vtkActor *skin = vtkActor::New(); + skin->SetMapper(skinMapper); + + // An outline provides context around the data. + // + vtkOutlineFilter *outlineData = vtkOutlineFilter::New(); + outlineData->SetInputConnection(v16->GetOutputPort()); + vtkPolyDataMapper *mapOutline = vtkPolyDataMapper::New(); + mapOutline->SetInputConnection(outlineData->GetOutputPort()); + vtkActor *outline = vtkActor::New(); + outline->SetMapper(mapOutline); + outline->GetProperty()->SetColor(0,0,0); + + // It is convenient to create an initial view of the data. The FocalPoint + // and Position form a vector direction. Later on (ResetCamera() method) + // this vector is used to position the camera to look at the data in + // this direction. + vtkCamera *aCamera = vtkCamera::New(); + aCamera->SetViewUp (0, 0, -1); + aCamera->SetPosition (0, 1, 0); + aCamera->SetFocalPoint (0, 0, 0); + aCamera->ComputeViewPlaneNormal(); + + // Actors are added to the renderer. An initial camera view is created. + // The Dolly() method moves the camera towards the FocalPoint, + // thereby enlarging the image. + aRenderer->AddActor(outline); + aRenderer->AddActor(skin); + aRenderer->SetActiveCamera(aCamera); + aRenderer->ResetCamera (); + aCamera->Dolly(1.5); + + // Set a background color for the renderer and set the size of the + // render window (expressed in pixels). + aRenderer->SetBackground(1,1,1); + renWin->SetSize(640, 480); + + // Note that when camera movement occurs (as it does in the Dolly() + // method), the clipping planes often need adjusting. Clipping planes + // consist of two planes: near and far along the view direction. The + // near plane clips out objects in front of the plane; the far plane + // clips out objects behind the plane. This way only what is drawn + // between the planes is actually rendered. + aRenderer->ResetCameraClippingRange (); + + // Initialize the event loop and then start it. + iren->Initialize(); + iren->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. + v16->Delete(); + skinExtractor->Delete(); + skinNormals->Delete(); + skinMapper->Delete(); + skin->Delete(); + outlineData->Delete(); + mapOutline->Delete(); + outline->Delete(); + aCamera->Delete(); + iren->Delete(); + renWin->Delete(); + aRenderer->Delete(); + + return 0; +} diff --git a/Surface_mesher/include/CGAL/vtkSurfaceMesherContourFilter.h b/Surface_mesher/include/CGAL/vtkSurfaceMesherContourFilter.h new file mode 100644 index 00000000000..f41f3080d21 --- /dev/null +++ b/Surface_mesher/include/CGAL/vtkSurfaceMesherContourFilter.h @@ -0,0 +1,156 @@ +// Copyright (c) 2008 GeometryFactory, Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Laurent Rineau + +#ifndef CGAL_VTK_SURFACE_MESHER_CONTOUR_FILTER_H +#define CGAL_VTK_SURFACE_MESHER_CONTOUR_FILTER_H + +#include + +class vtkCGALSurfaceMesherContourFilter : public vtkPolyDataAlgorithm +{ +public: + static vtkCGALSurfaceMesherContourFilter *New(); + vtkTypeMacro(vtkCGALSurfaceMesherContourFilter,vtkPolyDataAlgorithm); + void PrintSelf(ostream& os, vtkIndent indent); + + // Methods to set contour values + vtkSetMacro(Value,double); + vtkGetMacro(Value,double); + +protected: + vtkCGALSurfaceMesherContourFilter(); + ~vtkCGALSurfaceMesherContourFilter(); + + virtual int FillInputPortInformation(int port, vtkInformation *info); + + virtual int RequestData(vtkInformation *, vtkInformationVector **, vtkInformationVector *); + + double Value; + +private: + vtkCGALSurfaceMesherContourFilter(const vtkCGALSurfaceMesherContourFilter&); // Not implemented. + void operator=(const vtkCGALSurfaceMesherContourFilter&); // Not implemented. +}; + +// IMPLEMENTATION + +#include "vtkCellArray.h" +#include "vtkImageData.h" +#include "vtkInformation.h" +#include "vtkInformationVector.h" +#include "vtkObjectFactory.h" +#include "vtkPointData.h" +#include "vtkPolyData.h" + +#ifndef CGAL_USE_VTK +# define CGAL_USE_VTK +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkCGALSurfaceMesherContourFilter); + +vtkCGALSurfaceMesherContourFilter::vtkCGALSurfaceMesherContourFilter() +{ + Value = 0.; +} + +vtkCGALSurfaceMesherContourFilter::~vtkCGALSurfaceMesherContourFilter() +{ +} + +int +vtkCGALSurfaceMesherContourFilter:: +FillInputPortInformation(int, vtkInformation *info) +{ + info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData"); + return 1; +} + +void vtkCGALSurfaceMesherContourFilter::PrintSelf(ostream& os, + vtkIndent indent) +{ + this->Superclass::PrintSelf(os,indent); + + os << indent << "Value: " << this->Value << "\n"; +} + +int vtkCGALSurfaceMesherContourFilter::RequestData( + vtkInformation *vtkNotUsed(request), + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) +{ + // get the info objects + vtkInformation *inInfo = inputVector[0]->GetInformationObject(0); + vtkInformation *outInfo = outputVector->GetInformationObject(0); + + // get the input and ouptut + vtkImageData *inData = vtkImageData::SafeDownCast( + inInfo->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData *output = vtkPolyData::SafeDownCast( + outInfo->Get(vtkDataObject::DATA_OBJECT())); + + typedef CGAL::Surface_mesh_default_triangulation_3 Tr; + + // c2t3 + typedef CGAL::Complex_2_in_triangulation_3 C2t3; + + typedef Tr::Geom_traits GT; + typedef CGAL::Gray_level_image_3 Gray_level_image; + typedef CGAL::Implicit_surface_3 Surface_3; + + Tr tr; // 3D-Delaunay triangulation + C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation + + CGAL::Image_3 image; + if(!image.read_vtk_image_data(inData)) + return 0; + Gray_level_image gray_level_image(image, Value); + + GT::FT radius = std::max(image.xdim() * image.vx(), + std::max(image.ydim() * image.vy(), + image.zdim() * image.vz()) + ); + GT::Sphere_3 bounding_sphere(GT::Point_3(image.xdim() * image.vx()/2., + image.ydim() * image.vy()/2., + image.zdim() * image.vz()/2.), + radius*radius); + // definition of the surface, with 10^-2 as relative precision + Surface_3 surface(gray_level_image, bounding_sphere, 1e-5); + CGAL::Surface_mesh_default_criteria_3 criteria(30., + radius/100., + radius/500.); + // meshing surface, with the "manifold without boundary" algorithm + CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Manifold_tag()); + + CGAL::output_c2t3_to_vtk_polydata(c2t3, output); + output->Squeeze(); + + return 1; +} + +#endif // CGAL_VTK_SURFACE_MESHER_CONTOUR_FILTER_H