A sliced 3D volume writer, with contours

slice_contours.cpp

// --------------------------------------------------------------------------------------
// File:          slice_contours.cxx
// Date:          Nov 18, 2014
// Author:        code@oscaresteban.es (Oscar Esteban)
// Version:       1.0 beta
// License:       GPLv3 - 29 June 2007
// Short Summary:
// --------------------------------------------------------------------------------------
//
// Copyright (c) 2014, code@oscaresteban.es (Oscar Esteban)
// with Signal Processing Lab 5, EPFL (LTS5-EPFL)
// and Biomedical Image Technology, UPM (BIT-UPM)
// All rights reserved.
//
// This file is part of ACWE-Reg
//
// ACWE-Reg is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ACWE-Reg is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with ACWE-Reg.  If not, see <http://www.gnu.org/licenses/>.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//

#ifndef SLICE_CONTOURS_H_
#define SLICE_CONTOURS_H_


#include <boost/program_options.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/optional.hpp>
#include <boost/filesystem.hpp>

#include <vector>

#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkImageActor.h>
#include <vtkImageProperty.h>
#include <vtkImageMapper3D.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkImageData.h>
#include <vtkSphereSource.h>
#include <vtkMetaImageWriter.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkImageStencil.h>
#include <vtkPointData.h>
#include <vtkCutter.h>
#include <vtkPlane.h>
#include <vtkStripper.h>
#include <vtkLinearExtrusionFilter.h>
#include <vtkImageMapper.h>
#include <vtkImageSliceMapper.h>
#include <vtkImageSlice.h>
#include <vtkPolyDataReader.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkPNGWriter.h>
#include <vtkCamera.h>
#include <vtkCornerAnnotation.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkRescaleIntensityImageFilter.h>
#include <itkImageToVTKImageFilter.h>


#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkInteractorStyleImage.h>
#include <vtkRenderer.h>
#include <vtkWindowToImageFilter.h>

#if VTK_MAJOR_VERSION == 6
// These 2 lines are needed due to:
// http://www.vtk.org/Wiki/VTK/VTK_6_Migration/Factories_now_require_defines
// Taken from: https://gist.github.com/certik/5687727
#include <vtkAutoInit.h>
#define vtkRenderingCore_AUTOINIT 4(vtkInteractionStyle,vtkRenderingFreeType,vtkRenderingFreeTypeOpenGL,vtkRenderingOpenGL)
#define vtkRenderingVolume_AUTOINIT 1(vtkRenderingVolumeOpenGL)
#else
#include <vtkGraphicsFactory.h>
#include <vtkImagingFactory.h>
#endif

namespace bpo = boost::program_options;
namespace bfs = boost::filesystem;

static const unsigned int DIMENSION = 3;

typedef itk::Image<float, DIMENSION>                         ImageType;
typedef typename ImageType::Pointer                          ImagePointer;
typedef itk::ImageFileReader<ImageType>                      ReaderType;
typedef typename ReaderType::Pointer                         ReaderPointer;
typedef itk::ImageToVTKImageFilter<ImageType>                ConnectorType;
typedef itk::RescaleIntensityImageFilter< ImageType, ImageType > RescaleFilterType;

#endif /* SLICE_CONTOURS_H_ */



/**
 * This example generates a sphere, cuts it with a plane and, therefore, generates a circlular contour (vtkPolyData).
 * Subsequently a binary image representation (vtkImageData) is extracted from it. Internally vtkPolyDataToImageStencil and
 * vtkLinearExtrusionFilter are utilized. Both the circular poly data (circle.vtp) and the resultant image (labelImage.mhd)
 * are saved to disk.
 */
int main(int argc, char *argv[]) {

#if VTK_MAJOR_VERSION < 6
	// Setup offscreen rendering
	vtkSmartPointer<vtkGraphicsFactory> graphics_factory = vtkSmartPointer<vtkGraphicsFactory>::New();
	graphics_factory->SetOffScreenOnlyMode( 1);
	graphics_factory->SetUseMesaClasses( 1 );

	vtkSmartPointer<vtkImagingFactory> imaging_factory =  vtkSmartPointer<vtkImagingFactory>::New();
	imaging_factory->SetUseMesaClasses( 1 ); 
#endif

	std::string image;
	std::vector<std::string> surfaces, rsurfaces;
	int slice_num, nimages;
	std::string axisname = "axial";
	std::vector<int> axislist;

	bpo::options_description all_desc("Usage");
	all_desc.add_options()
		("help,h", "show help message")
		("input-image,i", bpo::value<std::string>(&image)->required(), "reference image file")
		("surfaces,S", bpo::value<std::vector<std::string>	>(&surfaces)->multitoken(), "surfaces (vtk files)")
		("ref-surfaces,R", bpo::value<std::vector<std::string>	>(&rsurfaces)->multitoken(), "reference surfaces (vtk files)")
		("num-slices,n", bpo::value < int >(&nimages)->default_value(14), "number of slices")
		("slice,s", bpo::value < int >(&slice_num)->default_value(-1), "slice number")
		("axis,a", bpo::value < std::vector<int> >(&axislist), "axes to be extracted")
		("all-axis,A", bpo::bool_switch(), "export all axes");

	bpo::variables_map vm;

	try {
		bpo::store(	bpo::parse_command_line( argc, argv, all_desc ), vm);

		if (vm.count("help")) {
			std::cout << all_desc << std::endl;
			return 1;
		}

		bpo::notify(vm);
	} catch ( bpo::error& e ) {
		std::cerr << "Error: " << e.what() << std::endl << std::endl;
		std::cerr << all_desc << std::endl;
		return 1;
	}

	ReaderPointer r = ReaderType::New();
	r->SetFileName(image);
	r->Update();

	RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
	rescaleFilter->SetInput(r->GetOutput());
	rescaleFilter->SetOutputMinimum(0);
	rescaleFilter->SetOutputMaximum(255);
	rescaleFilter->Update();

	ImagePointer im = rescaleFilter->GetOutput();
	ImageType::DirectionType dir;
	dir.SetIdentity();
	ImageType::PointType zero; zero.Fill(0.0);
	ImageType::PointType orig = im->GetOrigin();
	orig = im->GetDirection() * (orig - zero);
	im->SetOrigin(orig);
	im->SetDirection(dir);

	ConnectorType::Pointer connector = ConnectorType::New();
	connector->SetInput(im);
	connector->Update();
	vtkSmartPointer<vtkImageData> vtkim = connector->GetOutput();


	ImageType::PointType c;
	ImageType::SizeType s = im->GetLargestPossibleRegion().GetSize();
	ImageType::IndexType idx;

	for( size_t i = 0; i < DIMENSION; i++) {
		idx[i] = int(0.5*(s[i]-1));
	}

	if (axislist.size() == 0) {
		axislist.push_back(2);
	}

	if (vm.count("all-axis") && vm["all-axis"].as<bool>()) {
		axislist.empty();

		for (size_t aa = 0; aa < 3; aa++) axislist.push_back(aa);
	}

	size_t nsurf = surfaces.size();
	size_t nrsurf = rsurfaces.size();

	std::vector< vtkSmartPointer<vtkPolyData> > vp;
	for(size_t surf = 0; surf < surfaces.size(); surf++) {
		vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
		reader->SetFileName(surfaces[surf].c_str());
		reader->Update();
		vp.push_back(reader->GetOutput());
	}

	std::vector< vtkSmartPointer<vtkPolyData> > vpr;
	for(size_t surf = 0; surf < rsurfaces.size(); surf++) {
		vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
		reader->SetFileName(rsurfaces[surf].c_str());
		reader->Update();
		vpr.push_back(reader->GetOutput());
	}

	double totalims = 1.0 * (nimages+1);

	for(size_t ax = 0; ax < axislist.size(); ax++) {
		int axis = axislist[ax];
		if (axis==1) axisname = "coronal";
		if (axis==0) axisname = "sagittal";

		double normal[DIMENSION];
		for (size_t i = 0; i < DIMENSION; i++) {
			normal[i] = (i==axis)?1.0:0.0;
		}

		// Image slice: http://www.cmake.org/Wiki/VTK/Examples/Cxx/Images/ImageSliceMapper
		vtkSmartPointer<vtkImageSliceMapper> imageSliceMapper = vtkSmartPointer<vtkImageSliceMapper>::New();
		imageSliceMapper->SliceFacesCameraOn();
		imageSliceMapper->SetOrientation(axis);
#if VTK_MAJOR_VERSION < 6
		imageSliceMapper->SetInput(connector->GetOutput());
#else
		imageSliceMapper->SetInputData(connector->GetOutput());

#endif

		for (size_t sl = 0; sl < nimages; sl++) {
			double factor = ((sl+1) / totalims);
			idx[axis] = int( (s[axis]-1)* factor);
			im->TransformIndexToPhysicalPoint(idx, c);

			imageSliceMapper->SetSliceNumber(idx[axis]);

			// Setup renderers
			vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();

			for(size_t surf = 0; surf < rsurfaces.size(); surf++) {
				vtkSmartPointer<vtkCutter> cutter = vtkSmartPointer<vtkCutter>::New();
#if VTK_MAJOR_VERSION < 6
				cutter->SetInput(vpr[surf]);
#else
				cutter->SetInputData(vpr[surf]);
#endif
				vtkSmartPointer<vtkPlane> cutPlane = vtkSmartPointer<vtkPlane>::New();
				cutPlane->SetOrigin(c[0], c[1], c[2]);
				cutPlane->SetNormal(normal);
				cutter->SetCutFunction(cutPlane);
				vtkSmartPointer<vtkStripper> stripper = vtkSmartPointer<vtkStripper>::New();
				stripper->SetInputConnection(cutter->GetOutputPort()); // valid circle
				stripper->Update();

				vtkSmartPointer<vtkPolyDataMapper> contourmapper = vtkSmartPointer<vtkPolyDataMapper>::New();
				contourmapper->SetInputConnection(stripper->GetOutputPort());
				vtkSmartPointer<vtkActor> outlineActor = vtkSmartPointer<vtkActor>::New();
				outlineActor->SetMapper(contourmapper);
				outlineActor->GetProperty()->SetColor(1.0,1.0,0.2);
				outlineActor->GetProperty()->SetLineWidth(2);
				renderer->AddActor(outlineActor);
			}

			for(size_t surf = 0; surf < surfaces.size(); surf++) {
				vtkSmartPointer<vtkCutter> cutter = vtkSmartPointer<vtkCutter>::New();
#if VTK_MAJOR_VERSION < 6
				cutter->SetInput(vp[surf]);
#else
				cutter->SetInputData(vp[surf]);
#endif
				vtkSmartPointer<vtkPlane> cutPlane = vtkSmartPointer<vtkPlane>::New();
				cutPlane->SetOrigin(c[0], c[1], c[2]);
				cutPlane->SetNormal(normal);
				cutter->SetCutFunction(cutPlane);
				vtkSmartPointer<vtkStripper> stripper = vtkSmartPointer<vtkStripper>::New();
				stripper->SetInputConnection(cutter->GetOutputPort()); // valid circle
				stripper->Update();

				vtkSmartPointer<vtkPolyDataMapper> contourmapper = vtkSmartPointer<vtkPolyDataMapper>::New();
				contourmapper->SetInputConnection(stripper->GetOutputPort());
				vtkSmartPointer<vtkActor> outlineActor = vtkSmartPointer<vtkActor>::New();
				outlineActor->SetMapper(contourmapper);
				outlineActor->GetProperty()->SetColor(0.5,0.6,1.0);
				outlineActor->GetProperty()->SetLineWidth(2);
				renderer->AddActor(outlineActor);
			}

			vtkSmartPointer<vtkImageSlice> imageSlice = vtkSmartPointer<vtkImageSlice>::New();
			imageSlice->GetProperty()->SetInterpolationTypeToNearest();
			imageSlice->SetMapper(imageSliceMapper);
			// imageSlice->SetDisplayLocationToBackground();
			renderer->AddViewProp(imageSlice);

			// Setup render window
			int viewExtent[2];
			switch(axis) {
			case 0:
				viewExtent[0] = s[1]*5;
				viewExtent[1] = s[2]*5;
				break;
			case 1:
				viewExtent[0] = s[0]*5;
				viewExtent[1] = s[2]*5;
				break;
			default:
				viewExtent[0] = s[0]*5;
				viewExtent[1] = s[1]*5;
				break;
			}

			std::stringstream slicenumtext;
			slicenumtext << idx[axis];
			vtkSmartPointer<vtkTextActor> textActor = vtkSmartPointer<vtkTextActor>::New();
			textActor->GetTextProperty()->SetFontSize(40);
			textActor->SetPosition( int(0.25*viewExtent[0]), int(0.25*viewExtent[1]));
			textActor->SetInput( slicenumtext.str().c_str() );
			renderer->AddActor2D( textActor );

	//		vtkSmartPointer<vtkCornerAnnotation> cornerAnnotation = vtkSmartPointer<vtkCornerAnnotation>::New();
	//		cornerAnnotation->SetLinearFontScaleFactor(2);
	//		cornerAnnotation->SetNonlinearFontScaleFactor(1);
	//		cornerAnnotation->SetMaximumFontSize(30);
	//		cornerAnnotation->SetText( 0, slicenumtext.str().c_str() );
	//		renderer->AddViewProp(cornerAnnotation);

			vtkSmartPointer<vtkCamera> cam = renderer->GetActiveCamera();
			cam->SetPosition(normal);
			renderer->ResetCamera();

			vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
			renderWindow->SetOffScreenRendering(1);
			renderWindow->AddRenderer(renderer);


			// Setup render window interactor
		//	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
		//	vtkSmartPointer<vtkInteractorStyleImage> style = vtkSmartPointer<vtkInteractorStyleImage>::New();
		//	renderWindowInteractor->SetInteractorStyle(style);
		//	// Render and start interaction
		//	renderWindowInteractor->SetRenderWindow(renderWindow);
		//	renderWindowInteractor->Initialize();
		//
		//	renderWindowInteractor->Start();
		//
			vtkSmartPointer<vtkWindowToImageFilter> windowToImageFilter =
					vtkSmartPointer<vtkWindowToImageFilter>::New();
			windowToImageFilter->SetInput(renderWindow);
			windowToImageFilter->Update();

			vtkSmartPointer<vtkPNGWriter> writer = vtkSmartPointer<vtkPNGWriter>::New();
			writer->SetInputConnection(windowToImageFilter->GetOutputPort());

			std::stringstream ss;
			ss << axisname << std::setw(4) << std::setfill('0') << idx[axis] << ".png";
			writer->SetFileName(ss.str().c_str());
			writer->Write();
		}
	}

	return EXIT_SUCCESS;
}