/*=========================================================================

  Program:   FusionViewer
  Module:    $RCSfile: ItkImageProvider.cpp,v $
  Language:  C++
  Date:      $Date: 2008/01/11 20:29:35 $
  Version:   $Revision: 1.2 $
  
  Copyright (c) Insightful Corporation.  All rights reserved.
  See Copyright.txt 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.

=========================================================================*/

/**\class ImageProvider
 *\brief Adapter for itk::Image so that any 3D itk::Image instantiation can be
 *	displayed on a 32-bit color screen.
 */
 
//
//	Type Hierarchy
//
//	             ImageProvider
//	                  |
//	          ImageProviderPrivate
//	                  |
//	ImageProviderHelperImpl<PixelType, ScaleType>
//
//	ImageProvider is the public interface.  Users call factory methods in
//	ItkImageProvider for construction.  ImageProviderPrivate adds some virtual
//	methods we would like to access from the factory methods.  We do not want
//	these methods in the public interface since they have ITK types as
//	parameters.  ImageProviderHelperImpl is templated on each pixel type we
//	would like to use and does the work of converting the native pixel type
//	to 8 or 12 bit pixels.

#include "stdafx.h"

#include <stdio.h>
#include <string>
#include <sstream>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkImageIOBase.h>
#include <itkImageSeriesReader.h>
#include <itkImageSliceConstIteratorWithIndex.h>
#include <itkMinimumMaximumImageCalculator.h>
#include "ImageProvider.h"
#include "ItkImageProvider.h"
#include "ItkIOException.h"

#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"


#ifdef I2TK_IO
	#include <itkExtendedImageIOFactory.h>
#endif

#define DIMENSION	3


// --------------------------------------------------------------------------
//	Private classes
// --------------------------------------------------------------------------

/** \class ImageProviderPrivate
 * \brief Adds virtual methods to ImageProvider that we do not want exposed.
 */

class ImageProviderPrivate : public ImageProvider
{
public:
	virtual void	Import(int *size, float *pixelSpacing, void *pixels) = 0;
	virtual void	OpenFile(itk::ImageIOBase::Pointer ioBase, const char *filename) = 0;
	virtual void	OpenFile(itk::ImageIOBase::Pointer ioBase, const std::vector< std::string >& filenames) = 0;
	virtual void	SaveImage(const char *filename) = 0;
};

/** \class ImageProviderHelperImpl
 * \brief Implementation of the ImageProvider interface templated on the image pixel
 *	type (T).  This is a private subclass of ImageProvider.
 */

template <class T, class ScaleType>
class ImageProviderHelperImpl : public ImageProviderPrivate
{
public:
	ImageProviderHelperImpl();
	virtual ~ImageProviderHelperImpl();
	
	virtual void LoadSlice(int firstDirection, int secondDirection, int slice,
						   IntermediatePixelType *outBuffer);
	
	virtual void LoadSlice(int firstDirection, int secondDirection, int slice,
						   int startX, int startY, int width, int height,
						   float offset, float slope, GrayPixelType *colormap,
						   bool invert, GrayPixelType *outBuffer);

	virtual void MeasureROI(int startX, int startY, int width, int height,
							int firstDirection, int secondDirection, int slice,
							double& mean, double& stdDev, int& numPixels);
	
	virtual int MeasureLineProfile(int x1, int y1, int z1, int x2, int y2, int z2,
								   double *elems, int maxNumElems);
	
	virtual double GetPixel(int x, int y, int z);
	virtual void *GetPixelBuffer()	{ return m_image->GetBufferPointer(); }
	virtual void PixelUpdate()		{ CalculateScaling(); }

	virtual const std::type_info& GetPixelType() { return typeid(T); }
	virtual double GetMinValue()	{ return m_minValue; }
	virtual double GetMaxValue()	{ return m_maxValue; }
	
	virtual int GetSize(int i)		{ return m_size[i]; }
	virtual float GetPixelSpacing(int i) { return m_spacing[i]; }
	
	virtual bool IsByteType()		{ return m_fByteType; }

	void SetNoScale() { m_fByteType = true; }

	virtual void Import(int *size, float *pixelSpacing, void *pixels);
	virtual void OpenFile(itk::ImageIOBase::Pointer ioBase, const char *filename);
	virtual void OpenFile(itk::ImageIOBase::Pointer ioBase, const std::vector< std::string >& filenames);
	virtual void SaveImage(const char *filename);

protected:
	typedef itk::Image< T, DIMENSION >  ImageType3D;
	typedef itk::ImageSliceConstIteratorWithIndex< ImageType3D >	SliceIteratorType;
	
	typename ImageType3D::Pointer	m_image;				// ITK image
	ScaleType   m_slope;				// pixel value scaling (i.e., window and level)
	T			m_offset;
	bool		m_useFixedPoint;		// true if fixed point scaling should be used
	T			m_minValue, m_maxValue;	// minimum and maximum pixel values
	int			m_width, m_height, m_numSlices;
	int			m_size[DIMENSION];
	float		m_spacing[DIMENSION];
	bool		m_fByteType;
	
	void	CalculateScaling();
	
	void	CopySlice(SliceIteratorType inputIt, IntermediatePixelType *outBuffer);
	void	CopySlice(SliceIteratorType inputIt, float offset, float slope,
					  GrayPixelType *colormap, bool invert, GrayPixelType *outBuffer);
	
	void	SetImageIndex(typename ImageType3D::IndexType &idx,
						  const typename ImageType3D::SizeType &size,
						  int i, int newValue);
};

// Helper functions
ImageProviderPrivate *AllocateImage(const std::type_info &pixelType);
long PinValue(long value, long lowerBound, long upperBound);


// --------------------------------------------------------------------------
//	ItkImageProvider namespace
// --------------------------------------------------------------------------

// --------------------------------------------------------------------------
//		OpenFile
// --------------------------------------------------------------------------

/** \class OpenFile
 * \brief Open an image file.
 */
ImageProvider *ItkImageProvider::OpenFile(const char *filename, const std::type_info *pixelType)
{
	ImageProviderPrivate *pImage = NULL;

	try
	{
		#ifdef I2TK_IO
			itk::ExtendedImageIOFactory::RegisterBuiltInFactories();
		#endif

		// Determine the type of data in the image file and create the appropriate
		// reader implementation
		itk::ImageIOBase::Pointer	ioBase = itk::ImageIOFactory::CreateImageIO(filename, itk::ImageIOFactory::ReadMode);
		if (!ioBase)
		{
			std::stringstream errStr;
			errStr << "Could not open file " << filename << std::ends;
			throw ItkIOException(errStr.str().c_str());
		}

		ioBase->SetFileName(filename);
		ioBase->ReadImageInformation();
		
		if (!pixelType)
			pImage = AllocateImage(ioBase->GetComponentTypeInfo());
		else
			pImage = AllocateImage(*pixelType);

		pImage->OpenFile(ioBase, filename);
	}
	catch (const itk::ExceptionObject &err)
	{
		std::stringstream errStr;
		errStr << "Could not open file " << filename << std::endl << err.GetDescription() << std::ends;
		throw ItkIOException(errStr.str().c_str());
	}

	return pImage;
}

/** \class OpenFile
 * \brief Open an image volume.
 */

ImageProvider *ItkImageProvider::OpenFile(const std::vector< std::string >& filenames,
										  const std::type_info *pixelType)
{
	ImageProviderPrivate *pImage = NULL;

	try
	{
		#ifdef I2TK_IO
			itk::ExtendedImageIOFactory::RegisterBuiltInFactories();
		#endif
		
		// Determine the type of data in the image file and create the appropriate
		// reader implementation
		itk::ImageIOBase::Pointer	ioBase = itk::ImageIOFactory::CreateImageIO(filenames.front().c_str(), itk::ImageIOFactory::ReadMode);
		if (!ioBase)
		{
			std::stringstream errStr;
			errStr << "Could not open file sequence starting with " << filenames[0] << std::ends;
			throw ItkIOException(errStr.str().c_str());
		}
		
		ioBase->SetFileName(filenames.front().c_str());
		ioBase->ReadImageInformation();
		
		if (!pixelType)
			pImage = AllocateImage(ioBase->GetComponentTypeInfo());
		else
			pImage = AllocateImage(*pixelType);
		pImage->OpenFile(ioBase, filenames);
	}
	catch (const itk::ExceptionObject &err)
	{
		std::stringstream errStr;
		errStr << "Could not open file sequence starting with " << filenames[0] << std::endl << err.GetDescription() << std::ends;
		throw ItkIOException(errStr.str().c_str());
	}
	return pImage;
}

// --------------------------------------------------------------------------
//		CreateImage
// --------------------------------------------------------------------------

ImageProvider *ItkImageProvider::CreateImage(int *size, float *pixelSpacing,
											 const std::type_info& pixelType,
											 void *pixels)
{
	ImageProviderPrivate *pImage = AllocateImage(pixelType);
	pImage->Import(size, pixelSpacing, pixels);
	return pImage;
}


// --------------------------------------------------------------------------
//		SaveImage
// --------------------------------------------------------------------------

void ItkImageProvider::SaveImage(const char *filename, ImageProvider *pImage)
{
	ImageProviderPrivate *pPrivateImage = dynamic_cast<ImageProviderPrivate*>(pImage);
	if (pPrivateImage)
	{
		try
		{
			pPrivateImage->SaveImage(filename);
		}
		catch (const itk::ExceptionObject &err)
		{
			std::stringstream errStr;
			errStr << "Could not save image to file " << filename << std::endl << err.GetDescription() << std::ends;
			throw ItkIOException(errStr.str().c_str());
		}
	}
	else
	{
		throw ItkIOException("The image was not created by ItkImageProvider.");
	}
}


// --------------------------------------------------------------------------
//	File private helper functions
// --------------------------------------------------------------------------

// --------------------------------------------------------------------------
//		AllocateImage
// --------------------------------------------------------------------------
//	Instantiate the appropriate helper object based on the pixel type.

ImageProviderPrivate *AllocateImage(const std::type_info &pixelType)
{
	ImageProviderPrivate *pHelper = NULL;
	
	if (pixelType == typeid(float))
		pHelper = new ImageProviderHelperImpl<float, float>();
	else if ((pixelType == typeid(long)) || (pixelType == typeid(int)))
		pHelper = new ImageProviderHelperImpl<long, float>();
	else if ((pixelType == typeid(unsigned long)) || (pixelType == typeid(unsigned int)))
		pHelper = new ImageProviderHelperImpl<unsigned long, float>();
	else if (pixelType == typeid(short))
		pHelper = new ImageProviderHelperImpl<short, float>();
	else if (pixelType == typeid(unsigned short))
		pHelper = new ImageProviderHelperImpl<unsigned short, float>();
	else if (pixelType == typeid(char))
		pHelper = new ImageProviderHelperImpl<short, float>();
	else if (pixelType == typeid(unsigned char))
	{
		ImageProviderHelperImpl<unsigned char, float>
				*helper = new ImageProviderHelperImpl<unsigned char, float>();
		pHelper = helper;
		helper->SetNoScale();
	}
	else
	{
		throw ItkIOException((std::string("Unsupported pixel type: ") +
							  std::string(pixelType.name())).c_str());
	}

	return pHelper;
}


// --------------------------------------------------------------------------
//		PinValue
// --------------------------------------------------------------------------
//	Pin value to the lower and upper bounds.

long PinValue(long value, long lowerBound, long upperBound)
{
	if (value < lowerBound)
		value = lowerBound;
	if (value > upperBound)
		value = upperBound;
	
	return value;
}


// --------------------------------------------------------------------------
//	ImageProviderHelperImpl
// --------------------------------------------------------------------------

// --------------------------------------------------------------------------
//	ImageProviderHelperImpl constructor
// --------------------------------------------------------------------------

template <class T, class ScaleType>
ImageProviderHelperImpl<T, ScaleType>::ImageProviderHelperImpl()
	: m_fByteType(false)
{
}


// --------------------------------------------------------------------------
//		~ImageProviderHelperImpl
// --------------------------------------------------------------------------

template <class T, class ScaleType>
ImageProviderHelperImpl<T, ScaleType>::~ImageProviderHelperImpl()
{
}


// --------------------------------------------------------------------------
//		OpenFile
// --------------------------------------------------------------------------
//	Open the file at filename using ioBase.

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::OpenFile(itk::ImageIOBase::Pointer ioBase,
													 const char *filename)
{
	// Create an ITK file reading data graph and get its output
	typedef itk::ImageFileReader< ImageType3D > ReaderType;
	typename ReaderType::Pointer	fileReader = ReaderType::New();
	// ** BUG with ITK CVS version 01/12/2005: BMP reader does not read the image correctly if given an ImageIOBase
	 fileReader->SetImageIO(ioBase);
	fileReader->SetFileName(filename);
	fileReader->Update();
	
	m_image = fileReader->GetOutput();
	
	// Retrieve image dimensions and pixel spacing
	m_size[0] = ioBase->GetDimensions(0);
	m_size[1] = ioBase->GetDimensions(1);
			
	// Check for square pixels. Ignore decimal places smaller than 10e-5.
	if (static_cast<int>(ioBase->GetSpacing(0) * 10000.0f + 0.5f) != static_cast<int>(ioBase->GetSpacing(1) * 10000.0f + 0.5f))
		throw std::invalid_argument("Image pixels must be square.");
	m_spacing[0] = ioBase->GetSpacing(0);
	m_spacing[1] = ioBase->GetSpacing(1);
	
	if (ioBase->GetNumberOfDimensions() > 2)
	{
		m_spacing[2] = ioBase->GetSpacing(2);
		m_size[2] = ioBase->GetDimensions(2);
	}
	else
	{
		m_spacing[2] = 1.0f;
		m_size[2] = 1;
	}
	
	CalculateScaling();
}

// --------------------------------------------------------------------------
//		OpenFile
// --------------------------------------------------------------------------
//	Open the file at filename using ioBase.

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::OpenFile(itk::ImageIOBase::Pointer ioBase,
												    const std::vector< std::string >& filenames)
{
	// Create an ITK file reading data graph and get its output
	typedef itk::ImageSeriesReader< ImageType3D > ReaderType;
	typename ReaderType::Pointer	fileReader = ReaderType::New();
	// fileReader->SetImageIO(ioBase);	// Setting ioBase does not work for TIFF series
	fileReader->SetFileNames(filenames);
	fileReader->Update();
	
	m_image = fileReader->GetOutput();
	
	// Retrieve image dimensions and pixel spacing
	m_size[0] = ioBase->GetDimensions(0);
	m_size[1] = ioBase->GetDimensions(1);
	m_size[2] = static_cast<int>(filenames.size());
			
	// Check for square pixels. Ignore decimal places smaller than 10e-5.
	if (static_cast<int>(ioBase->GetSpacing(0) * 10000.0f + 0.5f) != static_cast<int>(ioBase->GetSpacing(1) * 10000.0f + 0.5f))
		throw std::invalid_argument("Image pixels must be square.");
	m_spacing[0] = ioBase->GetSpacing(0);
	m_spacing[1] = ioBase->GetSpacing(1);
	
	if (m_size[2] > 1)
	{
		m_spacing[2] = ioBase->GetSpacing(2);
	}
	else
	{
		m_spacing[2] = 1.0f;
	}
	
	CalculateScaling();
}

// --------------------------------------------------------------------------
//		SaveImage
// --------------------------------------------------------------------------

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::SaveImage(const char *filename)
{
	typedef itk::ImageFileWriter<ImageType3D> WriterType;
	typename WriterType::Pointer	fileWriter = WriterType::New();
	fileWriter->SetFileName(filename);
	fileWriter->SetInput(m_image);
	fileWriter->Update();
}


// --------------------------------------------------------------------------
//		Import
// --------------------------------------------------------------------------

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::Import(int *size, float *pixelSpacing, void *pixels)
{
	typename ImageType3D::SizeType		imageSize;
	typename ImageType3D::SpacingType	spacing;
	
	for (unsigned int i = 0; i < ImageType3D::ImageDimension; i++) {
		imageSize[i] = m_size[i] = size[i];
		spacing[i] = m_spacing[i] = pixelSpacing[i];
	}

	// Create the image
	typename ImageType3D::IndexType start;
	start.Fill(0);

	typename ImageType3D::RegionType region;
	region.SetSize(imageSize);
	region.SetIndex(start);

	m_image = ImageType3D::New();
	m_image->SetRegions(region);
	m_image->Allocate();

	m_image->SetSpacing(spacing);
	
	// Calculate size
	int		totalBytes = imageSize[0];
	for (unsigned int i = 1; i < ImageType3D::ImageDimension; i++)
		totalBytes *= imageSize[i];
	totalBytes = totalBytes * sizeof(typename ImageType3D::PixelType);
	
	if (pixels)
	{
		// Copy image data
		::memcpy(m_image->GetBufferPointer(), pixels, totalBytes);

		CalculateScaling();
	}
	else
	{
		::memset(m_image->GetBufferPointer(), 0, totalBytes);

		m_minValue = m_maxValue = m_offset = static_cast<T>(0.0);
		m_slope = static_cast<ScaleType>(0.0);
	}
}


// --------------------------------------------------------------------------
//		CalculateScaling
// --------------------------------------------------------------------------
//	Find the scaling factors needed to linearly scale the native pixel values
//	into 12 bits.

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::CalculateScaling()
{
	typedef itk::MinimumMaximumImageCalculator< ImageType3D >	MinMaxCalculatorType;
	typename MinMaxCalculatorType::Pointer	minMaxCalculator = MinMaxCalculatorType::New();
	minMaxCalculator->SetImage(m_image);
	minMaxCalculator->Compute();
	
	m_minValue = minMaxCalculator->GetMinimum();
	m_maxValue = minMaxCalculator->GetMaximum();
	
	if ((m_minValue >= static_cast<ScaleType>(0.0)) && (m_maxValue <= static_cast<ScaleType>(255.0)))
		m_fByteType = true;

	if (m_fByteType)
	{
		m_slope = static_cast<ScaleType>(1.0);
		m_offset = 0;
	}
	else if (m_maxValue == m_minValue)
	{
		// Trap here to avoid a slope of zero
		m_slope = static_cast<ScaleType>(1.0);
		m_offset = m_minValue;
	}
	else
	{
		m_slope = static_cast<ScaleType>(4095.0) / static_cast<ScaleType>((m_maxValue - m_minValue));
		m_offset = m_minValue;
	}
	
	// If the pixels are any integer type with a maximum value in the 12-bit range,
	// set a flag for performing the scaling using fixed point arithmetic
	if ((typeid(T) == typeid(float)) || (typeid(T) == typeid(double)))
	{
		m_useFixedPoint = false;
	}
	else
	{
		if (m_maxValue - m_minValue <= 4095)
			m_useFixedPoint = true;
		else
			m_useFixedPoint = false;
	}	
}


// --------------------------------------------------------------------------
//		CopySlice
// --------------------------------------------------------------------------

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::CopySlice(SliceIteratorType inputIt,
													  IntermediatePixelType *outBuffer)
{
	if ((m_offset == static_cast<ScaleType>(0.0)) && (m_slope == static_cast<ScaleType>(1.0)))
	{
		while (!inputIt.IsAtEndOfSlice())
		{
			while (!inputIt.IsAtEndOfLine())
			{
				*outBuffer++ = static_cast<IntermediatePixelType>(inputIt.Get());
				++inputIt;
			}
			
			inputIt.NextLine();
		}
	}
	else if (m_slope == static_cast<ScaleType>(1.0))
	{
		while (!inputIt.IsAtEndOfSlice())
		{
			while (!inputIt.IsAtEndOfLine())
			{
				*outBuffer++ = static_cast<IntermediatePixelType>(inputIt.Get() - m_offset);
				++inputIt;
			}
			
			inputIt.NextLine();
		}
	}
	else
	{
		while (!inputIt.IsAtEndOfSlice())
		{
			while (!inputIt.IsAtEndOfLine())
			{
				*outBuffer++ = static_cast<IntermediatePixelType>((static_cast<ScaleType>(inputIt.Get()) - m_offset) * m_slope +
																   static_cast<ScaleType>(0.5));
				++inputIt;
			}
			
			inputIt.NextLine();
		}
	}
}


// --------------------------------------------------------------------------
//		CopySlice
// --------------------------------------------------------------------------

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::CopySlice(SliceIteratorType inputIt,
													  float offset, float slope,
													  GrayPixelType *colormap,
													  bool invert,
													  GrayPixelType *outBuffer)
{
	const int kNumComponents = 3;		// number of components in the colormap
	GrayPixelType outVal;
	GrayPixelType *rgb;
	
	if (m_useFixedPoint)
	{
		// Calculate fixed point representation for the slope
		int fixSlope = static_cast<int>(slope * 65536.0f + 0.5f);
		
		while (!inputIt.IsAtEndOfSlice())
		{
			while (!inputIt.IsAtEndOfLine())
			{
				int		in_pixel = static_cast<int>(inputIt.Get() - offset);
				
				if (in_pixel <= 0)
				{
					outVal = 0;
				}
				else
				{
					in_pixel = (in_pixel * fixSlope) >> 16;
					
					if (in_pixel > 0xff)
						outVal = 0xff;
					else
						outVal = static_cast<GrayPixelType>(in_pixel);
				}
				
				if (invert)
					outVal = 255 - outVal;
				
				rgb = colormap + (outVal * kNumComponents) + 2;
				*outBuffer++ = *rgb--;
				*outBuffer++ = *rgb--;
				*outBuffer++ = *rgb;				
				*outBuffer++ = 255;
				++inputIt;
			}
			
			inputIt.NextLine();
		}
	}
	else
	{
		while (!inputIt.IsAtEndOfSlice())
		{
			while (!inputIt.IsAtEndOfLine())
			{
				ScaleType   in_pixel = static_cast<ScaleType>(inputIt.Get()) - offset;
				
				if (in_pixel <= 0)
					outVal = 0;
				else
				{
					in_pixel = in_pixel * slope + static_cast<ScaleType>(0.5);
					if (in_pixel > 255)
						outVal = 255;
					else
						outVal = static_cast<GrayPixelType>(in_pixel);
				}
				
				if (invert)
					outVal = 255 - outVal;
				
				rgb = colormap + (outVal * kNumComponents) + 2;				
				*outBuffer++ = *rgb--;
				*outBuffer++ = *rgb--;
				*outBuffer++ = *rgb;				
				*outBuffer++ = 255;
				++inputIt;
			}
			
			inputIt.NextLine();
		}
	}
}


// --------------------------------------------------------------------------
//		LoadSlice
// --------------------------------------------------------------------------
//	Set up an ITK image iterator to copy a slice.

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::LoadSlice(int firstDirection, int secondDirection, int slice,
													  IntermediatePixelType *outBuffer)
{
	try
	{
		if ((firstDirection < 0) || (secondDirection < 0) || (firstDirection > 2) ||
				(secondDirection > 2) || (firstDirection == secondDirection))
			throw std::invalid_argument("Invalid slice directions.");

		int	zDirection = 3 - (firstDirection + secondDirection);
		
		// Initialize region for getting a slice
		typename ImageType3D::RegionType	inputRegion = m_image->GetLargestPossibleRegion();
		typename ImageType3D::SizeType		size = inputRegion.GetSize();
		
		// Set the slice number to extract
		typename ImageType3D::IndexType	start = inputRegion.GetIndex();
		SetImageIndex(start, size, zDirection, slice);
		
		size[zDirection] = 1;
		
		// Set the extraction region
		typename ImageType3D::RegionType	desiredRegion;
		desiredRegion.SetSize(size);
		desiredRegion.SetIndex(start);
		
		// Setup slice iterator
		SliceIteratorType		inputIt(m_image, desiredRegion);
		inputIt.SetFirstDirection(firstDirection);
		inputIt.SetSecondDirection(secondDirection);
		
		// Use the iterator to copy the slice
		CopySlice(inputIt, outBuffer);
	}
	catch (const itk::ExceptionObject &err)
	{
		throw std::runtime_error(err.GetDescription());
	}
}


// --------------------------------------------------------------------------
//		LoadSlice
// --------------------------------------------------------------------------
//	Set up an ITK image iterator to copy a slice.

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::LoadSlice(int firstDirection, int secondDirection, int slice,
													  int startX, int startY, int width, int height,
													  float offset, float slope, GrayPixelType*colormap,
													  bool invert, GrayPixelType *outBuffer)
{
	try
	{
		if ((firstDirection < 0) || (secondDirection < 0) || (firstDirection > 2) ||
			(secondDirection > 2) || (firstDirection == secondDirection))
			throw std::invalid_argument("Invalid slice directions.");
		
		int	zDirection = 3 - (firstDirection + secondDirection);
		
		// Initialize region for getting a slice
		typename ImageType3D::RegionType	inputRegion = m_image->GetLargestPossibleRegion();
		typename ImageType3D::SizeType		size = inputRegion.GetSize();
		
		// Set the slice number to extract
		typename ImageType3D::IndexType	start = inputRegion.GetIndex();
		SetImageIndex(start, size, zDirection, slice);
		SetImageIndex(start, size, firstDirection, startX);
		SetImageIndex(start, size, secondDirection, startY);
		
		size[zDirection] = 1;
		if (width < static_cast<int>(size[firstDirection]))
			size[firstDirection] = width;
		if (height < static_cast<int>(size[secondDirection]))
			size[secondDirection] = height;
		
		// Set the extraction region
		typename ImageType3D::RegionType	desiredRegion;
		desiredRegion.SetSize(size);
		desiredRegion.SetIndex(start);
		
		// Setup slice iterator
		SliceIteratorType		inputIt(m_image, desiredRegion);
		inputIt.SetFirstDirection(firstDirection);
		inputIt.SetSecondDirection(secondDirection);
		
		// Use the iterator to copy the slice
		CopySlice(inputIt, offset, slope, colormap, invert, outBuffer);
	}
	catch (const itk::ExceptionObject &err)
	{
		throw std::runtime_error(err.GetDescription());
	}
}


// --------------------------------------------------------------------------
//		MeasureROI
// --------------------------------------------------------------------------

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::MeasureROI(int startX, int startY, int width, int height,
													   int firstDirection, int secondDirection, int slice,
													   double& mean, double& stdDev, int& numPixels)
{
	try
	{
		if ((firstDirection < 0) || (secondDirection < 0) || (firstDirection > 2) ||
			(secondDirection > 2) || (firstDirection == secondDirection))
			throw std::invalid_argument("Invalid slice directions.");
		
		int	zDirection = 3 - (firstDirection + secondDirection);
		
		// Initialize region for getting a slice
		typename ImageType3D::RegionType	inputRegion = m_image->GetLargestPossibleRegion();
		typename ImageType3D::SizeType		size = inputRegion.GetSize();
		
		// Set the slice number to extract
		typename ImageType3D::IndexType	start = inputRegion.GetIndex();
		SetImageIndex(start, size, zDirection, slice);
		SetImageIndex(start, size, firstDirection, startX);
		SetImageIndex(start, size, secondDirection, startY);
		
		size[zDirection] = 1;
		if (width < static_cast<int>(size[firstDirection]))
			size[firstDirection] = width;
		if (height < static_cast<int>(size[secondDirection]))
			size[secondDirection] = height;
		
		// Set the extraction region
		typename ImageType3D::RegionType	desiredRegion;
		desiredRegion.SetSize(size);
		desiredRegion.SetIndex(start);
		
		// Setup slice iterator
		SliceIteratorType		inputIt(m_image, desiredRegion);
		inputIt.SetFirstDirection(firstDirection);
		inputIt.SetSecondDirection(secondDirection);

		// Use the iterator to measure pixel values
		double totalValue = 0.0;
		double totalSquaredValues = 0.0;
		int	pixelCount = 0;
		while (!inputIt.IsAtEndOfSlice())
		{
			while (!inputIt.IsAtEndOfLine())
			{
				double val = inputIt.Get();
				totalValue += val;
				totalSquaredValues += val * val;

				// ::fprintf(f, "%g\n", val);

				++inputIt;
				pixelCount++;
			}
			
			inputIt.NextLine();
		}
		
		mean = totalValue / static_cast<double>(pixelCount);
		stdDev = ::sqrt((pixelCount * totalSquaredValues - totalValue  * totalValue) / (pixelCount * pixelCount));
		numPixels = pixelCount;
	}
	catch (const itk::ExceptionObject &err)
	{
		throw std::runtime_error(err.GetDescription());
	}
}


/* find maximum of a and b */
#define MAX(a,b) (((a)>(b))?(a):(b))

/* absolute value of a */
#define ABS(a) (((a)<0) ? -(a) : (a))

/* take sign of a, either -1, 0, or 1 */
#define ZSGN(a) (((a)<0) ? -1 : (a)>0 ? 1 : 0)


// --------------------------------------------------------------------------
//		MeasureLineProfile
// --------------------------------------------------------------------------
//	3d Bresenham line drawing algorithm, from ftp://ftp.uu.net/usenet/comp.sources.unix/volume26/line3d.

template <class T, class ScaleType>
int ImageProviderHelperImpl<T, ScaleType>::MeasureLineProfile(int x1, int y1, int z1, int x2, int y2, int z2,
															  double *elems, int maxNumElems)
{
	int i = 0;
	
	int xd, yd, zd;
	int x, y, z;
	int ax, ay, az;
	int sx, sy, sz;
	int dx, dy, dz;
	
	dx = x2 - x1;
	dy = y2 - y1;
	dz = z2 - z1;
	
	ax = ABS(dx) << 1;
	ay = ABS(dy) << 1;
	az = ABS(dz) << 1;
	
	sx = ZSGN(dx);
	sy = ZSGN(dy);
	sz = ZSGN(dz);
	
	x = x1;
	y = y1;
	z = z1;
	
	if (ax >= MAX(ay, az))            /* x dominant */
	{
		yd = ay - (ax >> 1);
		zd = az - (ax >> 1);
		for (;;)
		{
			if (i > maxNumElems)
				break;
			
			elems[i++] = GetPixel(x, y, z);
			
			if (x == x2)
				break;
			
			if (yd >= 0)
			{
				y += sy;
				yd -= ax;
			}
			
			if (zd >= 0)
			{
				z += sz;
				zd -= ax;
			}
			
			x += sx;
			yd += ay;
			zd += az;
		}
	}
	else if (ay >= MAX(ax, az))            /* y dominant */
	{
		xd = ax - (ay >> 1);
		zd = az - (ay >> 1);
		for (;;)
		{
			if (i > maxNumElems)
				break;
			
			elems[i++] = GetPixel(x, y, z);
			
			if (y == y2)
				break;
			
			if (xd >= 0)
			{
				x += sx;
				xd -= ay;
			}
			
			if (zd >= 0)
			{
				z += sz;
				zd -= ay;
			}
			
			y += sy;
			xd += ax;
			zd += az;
		}
	}
	else if (az >= MAX(ax, ay))            /* z dominant */
	{
		xd = ax - (az >> 1);
		yd = ay - (az >> 1);
		for (;;)
		{
			if (i > maxNumElems)
				break;
			
			elems[i++] = GetPixel(x, y, z);
			
			if (z == z2)
				break;
			
			if (xd >= 0)
			{
				x += sx;
				xd -= az;
			}
			
			if (yd >= 0)
			{
				y += sy;
				yd -= az;
			}
			
			z += sz;
			xd += ax;
			yd += ay;
		}
	}
	
	return i;
}


// --------------------------------------------------------------------------
//		SetImageIndex
// --------------------------------------------------------------------------

template <class T, class ScaleType>
void ImageProviderHelperImpl<T, ScaleType>::SetImageIndex(typename ImageType3D::IndexType &idx,
														  const typename ImageType3D::SizeType &size,
														  int i, int newValue)
{
	if (newValue < 0)
		idx[i] = 0;
	else if (newValue >= static_cast<int>(size[i]))
		idx[i] = size[i] - 1;
	else
		idx[i] = newValue;
}


// --------------------------------------------------------------------------
//		GetPixel
// --------------------------------------------------------------------------
//	Return the pixel value at x, y, z.

template <class T, class ScaleType>
double ImageProviderHelperImpl<T, ScaleType>::GetPixel(int x, int y, int z)
{
	// Ensure input coordinates are inside the image
	typename ImageType3D::IndexType	pixelIndex;
	typename ImageType3D::SizeType	size = m_image->GetLargestPossibleRegion().GetSize();
	pixelIndex[0] = PinValue(x, 0, size[0] - 1);
	pixelIndex[1] = PinValue(y, 0, size[1] - 1);
	pixelIndex[2] = PinValue(z, 0, size[2] - 1);
	
	return static_cast<double>(m_image->GetPixel(pixelIndex));
}