/*=========================================================================
  This code provided as part of 
  Brooks, R., Arbel, T. "Improving the OrientedImage class for use in the
  Registration Framework", Insight Journal, October, 2007

  Author: Rupert Brooks
  Institution: Centre for Intelligent Machines, McGill University

=========================================================================*/
//#define PRINTFUNCS
#include "itkFastOrientedImage.h"
#include "itkImage.h"
#include "itkOrientedImage.h"
#include "time.h"
#include <fstream>
#include <ostream>

typedef unsigned char PixelType;
// for any sort of real test, make speed test num large, but
// that can make ctest time out. So here we haev it small
//#define SPEEDTESTNUM 50000000
#define SPEEDTESTNUM 50000
#define NUMTESTPOINTS 2
#define TOL 1e-12
#define TEST2D 1
#define TEST3D 1
#define TEST4D 1
#define SPEEDTEST 1
#define ALLCOMBINATIONS 1

std::ostream * outfile=NULL;


template<class TFastImage,class TImage,class TPlainImage>
void testfunction(typename TFastImage::Pointer& fastImage,
				  typename TImage::Pointer& regularImage,
				  typename TPlainImage::Pointer& plainImage,
				  unsigned int& any ) {
	double fastregratio, fastplainratio;
					  // now set all the directions, spacing and origin
					  typename TImage::PointType inPoint1,inPoint2,fastPoint,regularPoint,origin,
						  testPoints[NUMTESTPOINTS],fastPoints[NUMTESTPOINTS],
						  regularPoints[NUMTESTPOINTS],plainPoints[NUMTESTPOINTS];
					  typename TImage::SpacingType spacing;
					  typename TImage::DirectionType direction;
					  typename TImage::IndexType inIndex1,inIndex2,fastIndex,regularIndex,
						  testIndexes[NUMTESTPOINTS],fastIndexes[NUMTESTPOINTS],
						  regularIndexes[NUMTESTPOINTS],plainIndexes[NUMTESTPOINTS];
					  const unsigned int Dimension=TImage::ImageDimension;
					  //double temp;
					  itk::ContinuousIndex<double,Dimension> inContIndex1,inContIndex2,
						  fastContIndex,regularContIndex,
						  testContIndexes[NUMTESTPOINTS],fastContIndexes[NUMTESTPOINTS],
						  regularContIndexes[NUMTESTPOINTS],plainContIndexes[NUMTESTPOINTS];
					  bool fastResult,regularResult,equal;

					  // set up the test points /indices
					  for (unsigned int i=0;i<Dimension;i++) {
						  inIndex1[i]=i*5+3; //inside the image
						  inIndex2[i]=i*10+3; //outside the image
						  inPoint1[i]=i*2.1235; //inside the image
						  inPoint2[i]=i*25.23; // outside the image
						  inContIndex1[i]=i*5+2.23489238; //inside the image
						  inContIndex2[i]=i*15+2.93847927; //outside the image
						  for (unsigned int j=0;j<NUMTESTPOINTS;j++) {
							  testPoints[j][i]=(i*2.1235)*(j+1);
							  testIndexes[j][i]=(i*5+3)*(j+1);
							  testContIndexes[j][i]=(i*5+2.23489238)*(j+1);
						  }
					  }	// try the images BEFORE actually setting anything

					  //std::cout<<"Testing TCITPP inside"<<std::endl;
					  fastImage->TransformContinuousIndexToPhysicalPoint(inContIndex1,fastPoint);
					  regularImage->TransformContinuousIndexToPhysicalPoint(inContIndex1,regularPoint);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fabs(fastPoint[i]-regularPoint[i])>TOL) equal=false;
					  if(!equal) {
						  (*outfile)<<"TCITPP inside Failed test"<<fastPoint<<regularPoint<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						  //std::cout<<"OK"<<std::endl;
					  }

					  //std::cout<<"Testing TCITPP outside"<<std::endl;
					  fastImage->TransformContinuousIndexToPhysicalPoint(inContIndex2,fastPoint);
					  regularImage->TransformContinuousIndexToPhysicalPoint(inContIndex2,regularPoint);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fabs(fastPoint[i]-regularPoint[i])>TOL) equal=false;
					  if(!equal) {
						  (*outfile)<<"TCITPP outside Failed test"<<fastPoint<<regularPoint<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						  //std::cout<<"OK"<<std::endl;
					  }
					  //Speed test
#if SPEEDTEST
					  clock_t faststart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  fastImage->TransformContinuousIndexToPhysicalPoint(testContIndexes[i%NUMTESTPOINTS],fastPoints[(i*i)%NUMTESTPOINTS]);
					  }
					  clock_t fastend=clock();
					  clock_t regularstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  regularImage->TransformContinuousIndexToPhysicalPoint(testContIndexes[i%NUMTESTPOINTS],regularPoints[(i*i)%NUMTESTPOINTS]);
					  }
					  clock_t regularend=clock();
					  clock_t plainstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  plainImage->TransformContinuousIndexToPhysicalPoint(testContIndexes[i%NUMTESTPOINTS],plainPoints[(i*i)%NUMTESTPOINTS]);
					  }
					  clock_t plainend=clock();
					  (*outfile)<<"TCITPP ";
					  fastregratio=((double)(fastend-faststart))/((double)(regularend-regularstart));
					  fastplainratio=((double)(fastend-faststart))/((double)(plainend-plainstart));
					  (*outfile)<<"Fast: "<<((double)(fastend-faststart))/CLOCKS_PER_SEC
						  <<" Orient: "<<((double)(regularend-regularstart))/CLOCKS_PER_SEC
						  <<" ("<<fastregratio<<") "
						  <<" Image "<<((double)(plainend-plainstart))/CLOCKS_PER_SEC
						  <<" ("<<fastplainratio<<") "
						  <<std::endl;
#endif

					  //std::cout<<"Testing TITPP inside"<<std::endl;
					  fastImage->TransformIndexToPhysicalPoint(inIndex1,fastPoint);
					  regularImage->TransformIndexToPhysicalPoint(inIndex1,regularPoint);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fabs(fastPoint[i]-regularPoint[i])>TOL) equal=false;
					  if(!equal) {
						  (*outfile)<<"TCITPP inside Failed test"<<fastPoint<<regularPoint<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						  //std::cout<<"OK"<<std::endl;
					  }
					  //std::cout<<"Testing TITPP outside"<<std::endl;
					  fastImage->TransformIndexToPhysicalPoint(inIndex2,fastPoint);
					  regularImage->TransformIndexToPhysicalPoint(inIndex2,regularPoint);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fabs(fastPoint[i]-regularPoint[i])>TOL) equal=false;
					  if(!equal) {
						  (*outfile)<<"TITPP outside Failed test"<<fastPoint<<regularPoint<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						  //std::cout<<"OK"<<std::endl;
					  }

#if SPEEDTEST
					  //Speed test
					  faststart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  fastImage->TransformIndexToPhysicalPoint(testIndexes[i%NUMTESTPOINTS],fastPoints[(i*i)%NUMTESTPOINTS]);
					  }
					  fastend=clock();
					  regularstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  regularImage->TransformIndexToPhysicalPoint(testIndexes[i%NUMTESTPOINTS],regularPoints[(i*i)%NUMTESTPOINTS]);
					  }
					  regularend=clock();
					  plainstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  plainImage->TransformIndexToPhysicalPoint(testIndexes[i%NUMTESTPOINTS],plainPoints[(i*i)%NUMTESTPOINTS]);
					  }
					  plainend=clock();
					  //std::cout<<"Regular Start: "<<regularstart<<" End: "<<regularend<<std::endl;
					  (*outfile)<<"TITPP ";
					  fastregratio=((double)(fastend-faststart))/((double)(regularend-regularstart));
					  fastplainratio=((double)(fastend-faststart))/((double)(plainend-plainstart));
					  (*outfile)<<"Fast: "<<((double)(fastend-faststart))/CLOCKS_PER_SEC
						  <<" Orient: "<<((double)(regularend-regularstart))/CLOCKS_PER_SEC
						  <<" ("<<fastregratio<<") "
						  <<" Image "<<((double)(plainend-plainstart))/CLOCKS_PER_SEC
						  <<" ("<<fastplainratio<<") "
						  <<std::endl;
#endif

					  //(*outfile)<<"Testing TPPTI inside"<<std::endl;
					  fastResult=fastImage->TransformPhysicalPointToIndex(inPoint1,fastIndex);
					  regularResult=regularImage->TransformPhysicalPointToIndex(inPoint1,regularIndex);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fastIndex[i]!=regularIndex[i]) equal=false;
					  if(!equal ||(fastResult!=regularResult)) {
						  (*outfile)<<"TPPTI inside Failed test"<<fastIndex<<regularIndex<<fastResult<<regularResult<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						  //(*outfile)<<"OK"<<std::endl;
					  }

					  //(*outfile)<<"Testing TPPTI outside"<<std::endl;
					  fastResult=fastImage->TransformPhysicalPointToIndex(inPoint2,fastIndex);
					  regularResult=regularImage->TransformPhysicalPointToIndex(inPoint2,regularIndex);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fastIndex[i]!=regularIndex[i]) equal=false;
					  if(!equal ||(fastResult!=regularResult)) {
						  (*outfile)<<"TPPTI outside Failed test"<<fastIndex<<regularIndex<<fastResult<<regularResult<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						  //(*outfile)<<"OK"<<std::endl;
					  }

#if SPEEDTEST
					  //Speed test
					  faststart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  fastImage->TransformPhysicalPointToIndex(testPoints[i%NUMTESTPOINTS],fastIndexes[(i*i)%NUMTESTPOINTS]);
					  }
					  fastend=clock();
					  regularstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  regularImage->TransformPhysicalPointToIndex(testPoints[i%NUMTESTPOINTS],regularIndexes[(i*i)%NUMTESTPOINTS]);
					  }
					  regularend=clock();
					  plainstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  plainImage->TransformPhysicalPointToIndex(testPoints[i%NUMTESTPOINTS],plainIndexes[(i*i)%NUMTESTPOINTS]);
					  }
					  plainend=clock();
					  (*outfile)<<"TPPTI ";
					  fastregratio=((double)(fastend-faststart))/((double)(regularend-regularstart));
					  fastplainratio=((double)(fastend-faststart))/((double)(plainend-plainstart));
					  (*outfile)<<"Fast: "<<((double)(fastend-faststart))/CLOCKS_PER_SEC
						  <<" Orient: "<<((double)(regularend-regularstart))/CLOCKS_PER_SEC
						  <<" ("<<fastregratio<<") "
						  <<" Image "<<((double)(plainend-plainstart))/CLOCKS_PER_SEC
						  <<" ("<<fastplainratio<<") "
						  <<std::endl;
#endif

					  //(*outfile)<<"Testing TPPTCI inside"<<std::endl;
					  fastResult=fastImage->TransformPhysicalPointToContinuousIndex(inPoint1,fastContIndex);
					  regularResult=regularImage->TransformPhysicalPointToContinuousIndex(inPoint1,regularContIndex);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fabs(fastContIndex[i]-regularContIndex[i])>TOL) equal=false;
					  if(!equal ||(fastResult!=regularResult)) {
						  (*outfile)<<"TPPTCI inside Failed test"<<fastContIndex<<regularContIndex<<fastResult<<regularResult<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						  //(*outfile)<<"OK"<<std::endl;
					  }

					  //(*outfile)<<"Testing TPPTCI outside"<<std::endl;
					  fastResult=fastImage->TransformPhysicalPointToContinuousIndex(inPoint2,fastContIndex);
					  regularResult=regularImage->TransformPhysicalPointToContinuousIndex(inPoint2,regularContIndex);
					  equal=true;
					  for(unsigned int i=0;i<Dimension;i++) if(fabs(fastContIndex[i]-regularContIndex[i])>TOL) equal=false;
					  if(!equal ||(fastResult!=regularResult)) {
						  (*outfile)<<"TPPTCI outside Failed test"<<fastContIndex<<regularContIndex<<fastResult<<regularResult<<std::endl;
						  (*outfile)<<"   Origin:  "<<fastImage->GetOrigin()<<" "<<regularImage->GetOrigin()<<std::endl;
						  (*outfile)<<"   Spacing: "<<fastImage->GetSpacing()<<" "<<regularImage->GetSpacing()<<std::endl;
						  (*outfile)<<"   Direction: "<<std::endl<<fastImage->GetDirection()<<std::endl<<regularImage->GetDirection()<<std::endl;
						  any++;
					  } else {
						 // (*outfile)<<"OK"<<std::endl;
					  }

#if SPEEDTEST
					  //Speed test
					  faststart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  fastImage->TransformPhysicalPointToContinuousIndex(testPoints[i%NUMTESTPOINTS],fastContIndexes[(i*i)%NUMTESTPOINTS]);
					  }
					  fastend=clock();
					  regularstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  regularImage->TransformPhysicalPointToContinuousIndex(testPoints[i%NUMTESTPOINTS],regularContIndexes[(i*i)%NUMTESTPOINTS]);
					  }
					  regularend=clock();
					  plainstart=clock();
					  for (unsigned int i=0;i<SPEEDTESTNUM;i++) {
						  plainImage->TransformPhysicalPointToContinuousIndex(testPoints[i%NUMTESTPOINTS],plainContIndexes[(i*i)%NUMTESTPOINTS]);
					  }
					  plainend=clock();
					  (*outfile)<<"TPPTCI ";
					  fastregratio=((double)(fastend-faststart))/((double)(regularend-regularstart));
					  fastplainratio=((double)(fastend-faststart))/((double)(plainend-plainstart));
					  (*outfile)<<"Fast: "<<((double)(fastend-faststart))/CLOCKS_PER_SEC
						  <<" Orient: "<<((double)(regularend-regularstart))/CLOCKS_PER_SEC
						  <<" ("<<fastregratio<<") "
						  <<" Image "<<((double)(plainend-plainstart))/CLOCKS_PER_SEC
						  <<" ("<<fastplainratio<<") "
						  <<std::endl;
#endif

				  }


				  int main( int argc, char *argv[] )
				  {
					  unsigned int any=0;
					  //(*outfile).precision(15);
					  //for(unsigned int dim=2;dim<5;dim++)
					  std::ofstream logfile;
					  if(argc==2) {
						  logfile.open(argv[1]);
						  outfile=&logfile;
					  } else {
						  outfile=&std::cout;
					  }
					  
#if TEST2D
					  {
						  const unsigned int Dimension=2;//dim;
						  (*outfile)<<"************** TESTING DIMENSION "<<Dimension<<" **************"<<std::endl;
						  // typedefs
						  typedef itk::FastOrientedImage<PixelType,Dimension> FastImageType;
						  typedef itk::OrientedImage<PixelType,Dimension> RegularImageType;
						  typedef itk::Image<PixelType,Dimension> PlainImageType;
						  typedef RegularImageType::RegionType RegionType;
						  // first declare identical images of each type.
						  FastImageType::Pointer fastImage=FastImageType::New();
						  RegularImageType::Pointer regularImage=RegularImageType::New();
						  PlainImageType::Pointer plainImage=PlainImageType::New();
						  RegionType region;
						  RegionType::SizeType size;
						  RegionType::IndexType index;
						  region.SetSize(size);
						  region.SetIndex(index);
						  index.Fill(0);
						  for (unsigned int i=0;i<Dimension;i++) {
							  size[i]=(i+1)*10;  // so it will be 10x20x30x40... etc
						  }
						  fastImage->SetRegions(region);
						  regularImage->SetRegions(region);
						  plainImage->SetRegions(region);

#if ALLCOMBINATIONS
						  (*outfile)<<"Nothing has been set yet"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif

						  RegularImageType::PointType fastPoint,regularPoint,origin;
						  RegularImageType::SpacingType spacing;
						  RegularImageType::DirectionType direction;
						  //RegularImageType::IndexType inIndex1,fastIndex,regularIndex;

						  // first set the origin, direction and spacing to zero/unity
						  origin.Fill(0.0);
						  spacing.Fill(1.0);
						  direction.SetIdentity();

						  (*outfile)<<"    ============================================="<<std::endl;
						  fastImage->SetOrigin(origin);
						  fastImage->SetSpacing(spacing);
						  fastImage->SetDirection(direction);
						  regularImage->SetOrigin(origin);
						  regularImage->SetSpacing(spacing);
						  regularImage->SetDirection(direction);
						  plainImage->SetOrigin(origin);
						  plainImage->SetSpacing(spacing);
						  plainImage->SetDirection(direction);
						  (*outfile)<<"Unit/Zero Everything"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(1.1);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonUnitSpacing"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(1.1);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
						  (*outfile)<<"NonUnitSpacingNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(1.0);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  double c=sqrt(2.0);
						  direction[0][0]=c;
						  direction[0][1]=c;
						  direction[1][0]=-c;
						  direction[1][1]=c;
						  fastImage->SetDirection(direction);
						  regularImage->SetDirection(direction);
						  plainImage->SetDirection(direction);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagonalDirectionNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(0.0);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagDirection"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(0.9);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagDirectionNonunitspacing"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(2.0);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
						  (*outfile)<<"NonDiagDirectionNonUnitSpacingNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

					  }
#endif
#if TEST3D
					  {
						  const unsigned int Dimension=3;//dim;
						  (*outfile)<<"************** TESTING DIMENSION "<<Dimension<<" **************"<<std::endl;
						  // typedefs
						  typedef itk::FastOrientedImage<PixelType,Dimension> FastImageType;
						  typedef itk::OrientedImage<PixelType,Dimension> RegularImageType;
						  typedef itk::Image<PixelType,Dimension> PlainImageType;
						  typedef RegularImageType::RegionType RegionType;
						  // first declare identical images of each type.
						  FastImageType::Pointer fastImage=FastImageType::New();
						  RegularImageType::Pointer regularImage=RegularImageType::New();
						  PlainImageType::Pointer plainImage=PlainImageType::New();
						  RegionType region;
						  RegionType::SizeType size;
						  RegionType::IndexType index;
						  region.SetSize(size);
						  region.SetIndex(index);
						  index.Fill(0);
						  for (unsigned int i=0;i<Dimension;i++) {
							  size[i]=(i+1)*10;  // so it will be 10x20x30x40... etc
						  }
						  fastImage->SetRegions(region);
						  regularImage->SetRegions(region);
						  plainImage->SetRegions(region);
#if ALLCOMBINATIONS
						  (*outfile)<<"Nothing has been set yet"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  RegularImageType::PointType fastPoint,regularPoint,origin;
						  RegularImageType::SpacingType spacing;
						  RegularImageType::DirectionType direction;
						  //RegularImageType::IndexType inIndex1,fastIndex,regularIndex;

						  // first set the origin, direction and spacing to zero/unity
						  origin.Fill(0.0);
						  spacing.Fill(1.0);
						  direction.SetIdentity();

						  (*outfile)<<"    ============================================="<<std::endl;
						  fastImage->SetOrigin(origin);
						  fastImage->SetSpacing(spacing);
						  fastImage->SetDirection(direction);
						  regularImage->SetOrigin(origin);
						  regularImage->SetSpacing(spacing);
						  regularImage->SetDirection(direction);
						  plainImage->SetOrigin(origin);
						  plainImage->SetSpacing(spacing);
						  plainImage->SetDirection(direction);
						  (*outfile)<<"Unit/Zero Everything"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(1.1);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonUnitSpacing"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(1.1);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
						  (*outfile)<<"NonUnitSpacingNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(1.0);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  double c=sqrt(2.0);
						  direction[0][0]=c;
						  direction[0][1]=c;
						  direction[1][0]=-c;
						  direction[1][1]=c;
						  fastImage->SetDirection(direction);
						  regularImage->SetDirection(direction);
						  plainImage->SetDirection(direction);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagonalDirectionNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(0.0);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagDirection"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(0.9);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagDirectionNonunitspacing"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(2.0);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
						  (*outfile)<<"NonDiagDirectionNonUnitSpacingNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

					  }
#endif
#if TEST4D
					  {
						  const unsigned int Dimension=4;//dim;
						  (*outfile)<<"************** TESTING DIMENSION "<<Dimension<<" **************"<<std::endl;
						  // typedefs
						  typedef itk::FastOrientedImage<PixelType,Dimension> FastImageType;
						  typedef itk::OrientedImage<PixelType,Dimension> RegularImageType;
						  typedef itk::Image<PixelType,Dimension> PlainImageType;
						  typedef RegularImageType::RegionType RegionType;
						  // first declare identical images of each type.
						  FastImageType::Pointer fastImage=FastImageType::New();
						  RegularImageType::Pointer regularImage=RegularImageType::New();
						  PlainImageType::Pointer plainImage=PlainImageType::New();
						  RegionType region;
						  RegionType::SizeType size;
						  RegionType::IndexType index;
						  region.SetSize(size);
						  region.SetIndex(index);
						  index.Fill(0);
						  for (unsigned int i=0;i<Dimension;i++) {
							  size[i]=(i+1)*10;  // so it will be 10x20x30x40... etc
						  }
						  fastImage->SetRegions(region);
						  regularImage->SetRegions(region);
						  plainImage->SetRegions(region);

#if ALLCOMBINATIONS
						  (*outfile)<<"Nothing has been set yet"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  RegularImageType::PointType fastPoint,regularPoint,origin;
						  RegularImageType::SpacingType spacing;
						  RegularImageType::DirectionType direction;
						  //RegularImageType::IndexType inIndex1,fastIndex,regularIndex;

						  // first set the origin, direction and spacing to zero/unity
						  origin.Fill(0.0);
						  spacing.Fill(1.0);
						  direction.SetIdentity();

						  (*outfile)<<"    ============================================="<<std::endl;
						  fastImage->SetOrigin(origin);
						  fastImage->SetSpacing(spacing);
						  fastImage->SetDirection(direction);
						  regularImage->SetOrigin(origin);
						  regularImage->SetSpacing(spacing);
						  regularImage->SetDirection(direction);
						  plainImage->SetOrigin(origin);
						  plainImage->SetSpacing(spacing);
						  plainImage->SetDirection(direction);
						  (*outfile)<<"Unit/Zero Everything"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(1.1);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonUnitSpacing"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(1.1);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
						  (*outfile)<<"NonUnitSpacingNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(1.0);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  double c=sqrt(2.0);
						  direction[0][0]=c;
						  direction[0][1]=c;
						  direction[1][0]=-c;
						  direction[1][1]=c;
						  fastImage->SetDirection(direction);
						  regularImage->SetDirection(direction);
						  plainImage->SetDirection(direction);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagonalDirectionNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(0.0);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagDirection"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  spacing.Fill(0.9);
						  fastImage->SetSpacing(spacing);
						  regularImage->SetSpacing(spacing);
						  plainImage->SetSpacing(spacing);
#if ALLCOMBINATIONS
						  (*outfile)<<"NonDiagDirectionNonunitspacing"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);
#endif
						  (*outfile)<<"    ============================================="<<std::endl;
						  origin.Fill(2.0);
						  fastImage->SetOrigin(origin);
						  regularImage->SetOrigin(origin);
						  plainImage->SetOrigin(origin);
						  (*outfile)<<"NonDiagDirectionNonUnitSpacingNonZeroOrigin"<<std::endl;
						  testfunction<FastImageType,RegularImageType,PlainImageType>(fastImage,regularImage,plainImage,any);

					  }
#endif
					  (*outfile)<<any<<" errors detected"<<std::endl;
					  return any;
				  }