#include "ImageProcessorTest.h"
#include "ImageProcessor.h"
#include "Log.h"
#include "TimeUtil.h"
#include <set>
#include <string.h>

CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(ImageProcessorTest, "XTCAVUnitTest");
CPPUNIT_REGISTRY_ADD_TO_DEFAULT("XTCAVUnitTest");

void ImageProcessorTest::setUp() {
  _singleLineTiff = TIFFOpen("TestImage-FirstLine.tif", "r");
  _doubleLineTiff = TIFFOpen("TestImage-DoubleLine.tif", "r");
  _fourLinesTiff = TIFFOpen("TestImage-FourLines.tif", "r");
  _roi1Tiff = TIFFOpen("TestImage-ROI-600-300x750-500.tif", "r");
  _roiBoxTiff = TIFFOpen("TestImage-Box-30x200-at-250x300.tif", "r");
  _trexTiff1 = TIFFOpen("TREX-Sample-2014-01-29-143735.tif","r");
  _trexTiff2 = TIFFOpen("TREX-Sample-2014-01-29-140544.tif","r");
  _trexTiff3 = TIFFOpen("TREX-Sample-2014-01-29-144155.tif","r");
  _backgroundTiff = TIFFOpen("TestImage-AverageBackground.tif","r");
  _roiStripesTiff = TIFFOpen("TestImage-ROIStripes.tif","r");
  _xtcavBackgroundTiff = TIFFOpen("XTCAV-background.tif", "r");
  _xtcavOnTiff = TIFFOpen("XTCAV-on.tif", "r");
  _testEnergyTiff = TIFFOpen("TestImage-Energy.tif", "r");

  if (_singleLineTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-FirstLine.tif file", false);
  }

  if (_doubleLineTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-DoubleLine.tif file", false);
  }

  if (_fourLinesTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-FourLines.tif file", false);
  }

  if (_roi1Tiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-ROI-600-300x750-500.tif file", false);
  }

  if (_roiBoxTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-Box-30x200-at-250x300.tif file", false);
  }

  if (_trexTiff1 == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TREX-Sample-2014-01-29-143735.tif file", false);
  }

  if (_trexTiff2 == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TREX-Sample-2014-01-29-140544.tif file", false);
  }

  if (_trexTiff3 == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TREX-Sample-2014-01-29-144155.tif file", false);
  }

  if (_backgroundTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-AverageBackground.tif file", false);
  }

  if (_roiStripesTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-ROIStripes.tif file", false);
  }

  if (_xtcavBackgroundTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open XTCAV-background.tif file", false);
  }

  if (_xtcavOnTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open XTCAV-on.tif file", false);
  }

  if (_testEnergyTiff == NULL) {
    CPPUNIT_ASSERT_MESSAGE("Failed to open TestImage-Energy.tif file", false);
  }

}

void ImageProcessorTest::tearDown() {
  TIFFClose(_singleLineTiff);
}

bool ImageProcessorTest::loadImage(TIFF *_tiffImage) {
  unsigned int w;
  unsigned int h;

  TIFFGetField(_tiffImage, TIFFTAG_IMAGEWIDTH, &w);
  TIFFGetField(_tiffImage, TIFFTAG_IMAGELENGTH, &h);

  Log::getInstance() << Log::showtime << Log::flagGlobal << Log::dpInfo
		     << "W=" << w << " H=" << h << Log::dp;

  _image.ndims = 2;

  _image.dims[0].size = w;
  _image.dims[0].offset = 0;
  _image.dims[0].binning = 0;
  _image.dims[0].reverse = 0;

  _image.dims[1].size = h;
  _image.dims[1].offset = 0;
  _image.dims[1].binning = 0;
  _image.dims[1].reverse = 0;

  _image.dataType = NDUInt16;
  _image.dataSize = w * h * 4;
  _image.pData = realloc(_image.pData, _image.dataSize);

  tsize_t scanline = TIFFScanlineSize(_tiffImage);
  Log::getInstance() << Log::showtime << Log::flagGlobal << Log::dpInfo
		     << "scanline=" << scanline << Log::dp;
  
  tdata_t buf;
  buf = _TIFFmalloc(scanline);
  uint16 *buf16;
  uint16 *pData16 = (uint16 *) _image.pData;

  _sum0 = 0;
  _sum1 = 0;
  _sum2 = 0;

  for (uint32 row = 0; row < h; ++row) {
    CPPUNIT_ASSERT(TIFFReadScanline(_tiffImage, buf, row) > 0);
    buf16 = (uint16 *) buf;
    for (unsigned int i = 0; i < w; ++i) {
      pData16[i] = buf16[i];
    } 

    _sum0 += pData16[0];
    _sum1 += pData16[0] * (row + 1);
    _sum2 += pData16[0]*pData16[0] * (row+1)*(row+1);

    if (row < 3 || row > h - 4 ||
	(h / 2 -1 <= row && row <= h / 2 + 1)) {
      Log::getInstance() << Log::showtime << Log::flagGlobal << Log::dpInfo
			 << "[row " << (int) row << "]: "
			 << (int) pData16[0] << " "
			 << (int) pData16[1] << " "
			 << (int) pData16[2] << " ... "
			 << (int) pData16[w/2-1] << " "
			 << (int) pData16[w/2] << " "
			 << (int) pData16[w/2+1] << " ... "	
			 << (int) pData16[w-3] << " "
			 << (int) pData16[w-2] << " "
			 << (int) pData16[w-1] << " "
			 << Log::dp;
    }
    pData16 += w;
  }

  _TIFFfree(buf);

  return true;
}


/**
 * @author L.Piccoli
 */
void ImageProcessorTest::testProcessReferenceSingleLine() {
  ImageProcessor p;
  
  CPPUNIT_ASSERT(loadImage(_singleLineTiff));
  CPPUNIT_ASSERT(p.processReference(&_image));
  /*
  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum1 (test) = " << _sum1
		     << "; Sum1 = " << p._referenceSum1[0] 
		     << Log::dp;

  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum2 (test) = " << _sum2
		     << "; Sum2 = " << p._referenceSum2[0] 
		     << Log::dp;


  CPPUNIT_ASSERT(p._referenceSum0[0] == _sum0);
  CPPUNIT_ASSERT(p._referenceSum1[0] == _sum1);
  CPPUNIT_ASSERT(p._referenceSum2[0] == _sum2);
  */
}

void ImageProcessorTest::testProcessReferenceDoubleLine() {
  ImageProcessor p;
  
  CPPUNIT_ASSERT(loadImage(_doubleLineTiff));
  CPPUNIT_ASSERT(p.processReference(&_image));
  /*
  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum1 (test) = " << _sum1
		     << "; Sum1 = " << p._referenceSum1[0] 
		     << Log::dp;

  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum2 (test) = " << _sum2
		     << "; Sum2 = " << p._referenceSum2[0] 
		     << Log::dp;

  CPPUNIT_ASSERT(p._referenceSum0[0] == _sum0);
  CPPUNIT_ASSERT(p._referenceSum1[0] == _sum1);
  CPPUNIT_ASSERT(p._referenceSum2[0] == _sum2);

  */
}

void ImageProcessorTest::testProcessReferenceFourLines() {
  ImageProcessor p;
  
  CPPUNIT_ASSERT(loadImage(_fourLinesTiff));
  CPPUNIT_ASSERT(p.processReference(&_image));
  /*
  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum1 (test) = " << _sum1
		     << "; Sum1 = " << p._referenceSum1[0] 
		     << Log::dp;

  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum2 (test) = " << _sum2
		     << "; Sum2 = " << p._referenceSum2[0] 
		     << Log::dp;
  CPPUNIT_ASSERT(p._referenceSum0[0] == _sum0);
  CPPUNIT_ASSERT(p._referenceSum1[0] == _sum1);
  CPPUNIT_ASSERT(p._referenceSum2[0] == _sum2);
  */
}

void ImageProcessorTest::testProcessDoubleLine() {
  ImageProcessor p;

  CPPUNIT_ASSERT(loadImage(_singleLineTiff));
  CPPUNIT_ASSERT(p.processReference(&_image));
  /*
  double _referenceSum1 = _sum1;
  double _referenceSum2 = _sum2;

  CPPUNIT_ASSERT(loadImage(_doubleLineTiff));
  CPPUNIT_ASSERT(p.process(&_image));

  CPPUNIT_ASSERT(p._sum0[0] == _sum0);
  CPPUNIT_ASSERT(p._sum1[0] == _sum1);
  CPPUNIT_ASSERT(p._sum2[0] == _sum2);

  double _sum3 = _referenceSum1 - _sum1;
  double _sum4 = _referenceSum2 - _sum2;
  */
}

void ImageProcessorTest::testLoadReference() {
  ImageProcessor p;

  CPPUNIT_ASSERT(p.loadReference("InvalidFile.tif") == false);
  CPPUNIT_ASSERT(p.loadReference("TestImage-FourLines.tif"));
}

/**
 * The reference file has two big gray boxes side-by-side
 * with pixel intensity of 14417. The boxes are 700w 200h
 * and 424w 500h. All other pixels are 0 intensity. The
 * average should be 14417 * 200 / 1024
 */
void ImageProcessorTest::testFindAverageIntensity() {
  ImageProcessor p;

  CPPUNIT_ASSERT(p.loadReference("TestImage-AverageBackground.tif"));
  //CPPUNIT_ASSERT(p.loadReference("TestImage-PrintedBeam.tif")); // average=8427

  unsigned long expectedAverage = 14417 * 200 / 1024;
  CPPUNIT_ASSERT(p.findAverageIntensity() == expectedAverage);
}

void ImageProcessorTest::testProcessLoadedReference() {
  ImageProcessor p;

  CPPUNIT_ASSERT(p.processReference() == false);
  CPPUNIT_ASSERT(p.loadReference("TestImage-FourLines.tif"));

  ImageProcessor p2;
  
  CPPUNIT_ASSERT(p.processReference());
  CPPUNIT_ASSERT(loadImage(_fourLinesTiff));
  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum1 (test) = " << _sum1
		     << "; Sum1 = " << p._referenceSum1[0] 
		     << Log::dp;

  Log::getInstance() << Log::showtime << Log::flagReference << Log::dpInfo
		     << "Sum2 (test) = " << _sum2
		     << "; Sum2 = " << p._referenceSum2[0] 
		     << Log::dp;
  CPPUNIT_ASSERT(p._referenceSum0[0] == _sum0);
  CPPUNIT_ASSERT(p._referenceSum1[0] == _sum1);
  CPPUNIT_ASSERT(p._referenceSum2[0] == _sum2);
}

void ImageProcessorTest::testSumDoubleLine() {
  ImageProcessor p;

  TimeAverage avg(1, "X/Y Sum Time");
  
  CPPUNIT_ASSERT(loadImage(_doubleLineTiff));
  avg.start();
  CPPUNIT_ASSERT(p.calculateSums(&_image));
  avg.end();
  avg.show();

  for (unsigned int i = 0; i < p._sumX.size(); i++) {
    CPPUNIT_ASSERT(p._sumX[i] == 130);
  }
  for (unsigned int i = 0; i < p._sumX.size(); i++) {
    if (p._sumY[i] != 0) {
      CPPUNIT_ASSERT(p._sumY[i] == 65 * 1024);
    }
  }
}

void ImageProcessorTest::testFindROIPointsX1() {
  ImageProcessor p;
  unsigned long bgThreshold = 10;
  unsigned expectedROI = 30; // this is the gray box size
  Point before;
  Point after;
  unsigned int scanInterval = expectedROI / 2;

  TimeAverage avg(1, "Find X ROI points Time");
  
  unsigned long avgFound;
  CPPUNIT_ASSERT(loadImage(_roiBoxTiff));
  avg.start();
  CPPUNIT_ASSERT(p.findROIPoints(&_image, scanInterval, bgThreshold, before, after, true, &avgFound));
  avg.end();

  int expectedBefore = 247;
  int expectedAfter = 282;

  std::cout << "left=" << before.x << " expected=" << expectedBefore << std::endl;
  std::cout << "right=" << after.x << " expected=" << expectedAfter << std::endl;

  CPPUNIT_ASSERT(before.x == expectedBefore);
  CPPUNIT_ASSERT(after.x == expectedAfter);
  avg.show();
}

// double the expected ROI size 
void ImageProcessorTest::testFindROIPointsX2() {
  ImageProcessor p;
  unsigned long bgThreshold = 10;
  unsigned expectedROI = 30; // this is the gray box size
  Point before;
  Point after;
  unsigned int scanInterval = expectedROI / 2;

  unsigned long avgFound;
  TimeAverage avg(1, "Find X ROI points Time");
  CPPUNIT_ASSERT(loadImage(_roiBoxTiff));

  avg.clear();
  expectedROI = 60; 
  scanInterval = expectedROI / 2;
  avg.start();
  CPPUNIT_ASSERT(p.findROIPoints(&_image, scanInterval, bgThreshold, before, after, true, &avgFound));
  avg.end();

  int expectedBefore = 247;
  int expectedAfter = 287;

  std::cout << "left=" << before.x << " expected=" << expectedBefore << std::endl;
  std::cout << "right=" << after.x << " expected=" << expectedAfter << std::endl;

  CPPUNIT_ASSERT(before.x == expectedBefore);
  CPPUNIT_ASSERT(after.x == expectedAfter);
  avg.show();
}

// Try in the other direction - box is 200px tall
void ImageProcessorTest::testFindROIPointsY1() {
  ImageProcessor p;
  unsigned long bgThreshold = 10;
  unsigned expectedROI = 200; // this is the gray box size
  Point before;
  Point after;
  unsigned int scanInterval = expectedROI / 2;

  unsigned long avgFound;
  TimeAverage avg(1, "Find Y ROI points Time");
  CPPUNIT_ASSERT(loadImage(_roiBoxTiff));

  avg.clear();
  expectedROI = 200; 
  scanInterval = expectedROI / 2;
  avg.start();
  CPPUNIT_ASSERT(p.findROIPoints(&_image, scanInterval, bgThreshold, before, after, false, &avgFound));
  avg.end();

  int expectedBefore = 262;
  int expectedAfter = 512;

  std::cout << "left=" << before.y << " expected=" << expectedBefore << std::endl;
  std::cout << "right=" << after.y << " expected=" << expectedAfter << std::endl;

  CPPUNIT_ASSERT(before.y == expectedBefore);
  CPPUNIT_ASSERT(after.y == expectedAfter);
  avg.show();
}

void ImageProcessorTest::testFindROIPointsY2() {
  ImageProcessor p;
  unsigned long bgThreshold = 10;
  unsigned expectedROI = 200; // this is the gray box size
  Point before;
  Point after;
  unsigned int scanInterval = expectedROI / 2;

  unsigned long avgFound;
  TimeAverage avg(1, "Find Y ROI points Time");
  CPPUNIT_ASSERT(loadImage(_roiBoxTiff));

  avg.clear();
  expectedROI = 400; 
  scanInterval = expectedROI / 2;
  avg.start();
  CPPUNIT_ASSERT(p.findROIPoints(&_image, scanInterval, bgThreshold, before, after, false, &avgFound));
  avg.end();

  int expectedBefore = 212;
  int expectedAfter = 512;

  std::cout << "left=" << before.y << " expected=" << expectedBefore << std::endl;
  std::cout << "right=" << after.y << " expected=" << expectedAfter << std::endl;

  CPPUNIT_ASSERT(before.y == expectedBefore);
  CPPUNIT_ASSERT(after.y == expectedAfter);
  avg.show();
}

void ImageProcessorTest::testROI1() {
  ImageProcessor p;
  unsigned long threshold = 150; // bg pixel intensity is less than 250
  int expectedROIX = 30;
  int expectedROIY = 100;

  TimeAverage avg(1, "TREX img ROI Time");
  
  CPPUNIT_ASSERT(loadImage(_trexTiff1));
  avg.start();
  CPPUNIT_ASSERT(p.findROI(&_image, threshold, expectedROIX, expectedROIY));
  avg.end();
  avg.show();

  CPPUNIT_ASSERT(p._bottomROI.x - p._topROI.x < expectedROIX * 2);
  CPPUNIT_ASSERT(p._bottomROI.y - p._topROI.y < expectedROIY * 2);
}

void ImageProcessorTest::testROI2() {
  ImageProcessor p;
  unsigned long threshold = 150; // bg pixel intensity is less than 250
  int expectedROIX = 35;
  int expectedROIY = 120;

  TimeAverage avg(1, "TREX img ROI Time");
  
  CPPUNIT_ASSERT(loadImage(_trexTiff2));
  avg.start();
  CPPUNIT_ASSERT(p.findROI(&_image, threshold, expectedROIX, expectedROIY));
  avg.end();
  avg.show();

  CPPUNIT_ASSERT(p._bottomROI.x - p._topROI.x < expectedROIX * 2);
  CPPUNIT_ASSERT(p._bottomROI.y - p._topROI.y < expectedROIY * 2);
}

void ImageProcessorTest::testROI3() {
  ImageProcessor p;
  unsigned long threshold = 150; // bg pixel intensity is less than 250
  int expectedROIX = 35;
  int expectedROIY = 120;

  TimeAverage avg(1, "TREX img ROI Time");
  
  CPPUNIT_ASSERT(loadImage(_trexTiff3));
  avg.start();
  CPPUNIT_ASSERT(p.findROI(&_image, threshold, expectedROIX, expectedROIY));
  avg.end();
  avg.show();

  CPPUNIT_ASSERT(p._bottomROI.x - p._topROI.x < expectedROIX * 2);
  CPPUNIT_ASSERT(p._bottomROI.y - p._topROI.y < expectedROIY * 2);
}

void ImageProcessorTest::testSumColumn() {
  ImageProcessor p;

  TimeAverage avg(1, "SumColumn Time");

  unsigned long expectedSum = 65;

  CPPUNIT_ASSERT(loadImage(_singleLineTiff));
  avg.start();
  unsigned long sum = p.sumColumn(&_image, 10);
  CPPUNIT_ASSERT(sum == expectedSum);
  avg.end();
  avg.show();

  sum = p.sumColumn(&_image, 5000);
  CPPUNIT_ASSERT(sum == 0);
}

void ImageProcessorTest::testSumRow() {
  ImageProcessor p;

  TimeAverage avg(1, "SumRow Time");

  CPPUNIT_ASSERT(loadImage(_singleLineTiff));
  unsigned long expectedSum = 65 * _image.dims[0].size;
  avg.start();
  unsigned long sum = p.sumRow(&_image, 0);
  CPPUNIT_ASSERT(sum == expectedSum);
  avg.end();
  avg.show();

  sum = p.sumRow(&_image, 5000);
  CPPUNIT_ASSERT(sum == 0);

  expectedSum = 0;
  sum = p.sumRow(&_image, 1);
  CPPUNIT_ASSERT(sum == expectedSum);
}

void ImageProcessorTest::testSaveROI() {
  ImageProcessor p;
  unsigned long threshold = 150; // bg pixel intensity is less than 250
  int expectedROIX = 30;
  int expectedROIY = 100;

  TimeAverage avg(1, "TREX save ROI Time");
  
  CPPUNIT_ASSERT(loadImage(_trexTiff1));
  CPPUNIT_ASSERT(p.findROI(&_image, threshold, expectedROIX, expectedROIY));
  avg.start();
  p.saveROI(&_image);
  avg.end();
  avg.show();

  /** Check if data has been copied correctly */
  char *srcPtr = (char *) (_image.pData);
  char *destPtr = (char *) (p._roiArray.pData);

  srcPtr += sizeof(epicsUInt16) * (p._topROI.y * _image.dims[0].size); // Skip first rows
  srcPtr += sizeof(epicsUInt16) * p._topROI.x; // Skip first cols
  
  int w = (p._bottomROI.x - p._topROI.x + 1);


  for (int i = p._topROI.y; i < p._bottomROI.y; ++i) {
    epicsUInt16 *srcPtr16 = (epicsUInt16 *) srcPtr;
    epicsUInt16 *destPtr16 = (epicsUInt16 *) destPtr;

    Log::getInstance() << Log::showtime << Log::flagGlobal << Log::dpInfo
		       << "SRC: " << (unsigned int) srcPtr16[0] << " "
		       << (unsigned int) srcPtr16[1] << " "
		       << (unsigned int) srcPtr16[2] << " ... "
		       << (unsigned int) srcPtr16[w/2-1] << " "
		       << (unsigned int) srcPtr16[w/2] << " "
		       << (unsigned int) srcPtr16[w/2+1] << " ... "	
		       << (unsigned int) srcPtr16[w-3] << " "
		       << (unsigned int) srcPtr16[w-2] << " "
		       << (unsigned int) srcPtr16[w-1] << " "
		       << Log::dp;

    Log::getInstance() << Log::showtime << Log::flagGlobal << Log::dpInfo
		       << "DST: " << (unsigned int) destPtr16[0] << " "
		       << (unsigned int) destPtr16[1] << " "
		       << (unsigned int) destPtr16[2] << " ... "
		       << (unsigned int) destPtr16[w/2-1] << " "
		       << (unsigned int) destPtr16[w/2] << " "
		       << (unsigned int) destPtr16[w/2+1] << " ... "	
		       << (unsigned int) destPtr16[w-3] << " "
		       << (unsigned int) destPtr16[w-2] << " "
		       << (unsigned int) destPtr16[w-1] << " "
		       << Log::dp;

    CPPUNIT_ASSERT(strncmp(srcPtr, destPtr, sizeof(epicsUInt16) * w) == 0);
    srcPtr += sizeof(epicsUInt16) * _image.dims[0].size; // Go to the next row
    destPtr += sizeof(epicsUInt16) * w;
  }

  CPPUNIT_ASSERT(p._bottomROI.x - p._topROI.x < expectedROIX * 2);
  CPPUNIT_ASSERT(p._bottomROI.y - p._topROI.y < expectedROIY * 2);
}

void ImageProcessorTest::testSliceROI1() {
  testSliceROI(30);
}

void ImageProcessorTest::testSliceROI2() {
  testSliceROI(7);
}

void ImageProcessorTest::testSliceROI3() {
  testSliceROI(13);
}

/**
 * The test image is a 300x200 horizontal striped image, with pixel
 * intensities of 0 and 19735. A row has either 0 intensity or
 * 19735 intensity pixels.
 */
void ImageProcessorTest::testSliceROI(int slices) {
  ImageProcessor p;

  TimeAverage avg(1, "Slice ROI Time");

  unsigned long intensity = 19735;

  CPPUNIT_ASSERT(loadImage(_roiStripesTiff));
  p._topROI.x = 0;
  p._topROI.y = 0;
  p._bottomROI.x = _image.dims[0].size - 1;
  p._bottomROI.y = _image.dims[1].size - 1;
  p.saveROI(&_image);

  p.setROISlices(slices);

  // Expected slice sum is 19735 multiplied be half of the rows
  // (the other half is zero) multiplied by the number of cols
  // each slice gets
  unsigned long expectedHProjection =
    (p._roiArray.dims[0].size / p._roiSlices) * intensity *
    (p._roiArray.dims[1].size / 2);

  // Each row of a slice is summed, generating a waveform
  // of either 0s or 19735 * number of cols per slice
  unsigned long expectedVProjection = 
    (p._roiArray.dims[0].size / p._roiSlices) * intensity;

  // If ROI can't be sliced in pieces of same size, last slice
  // gets more columns
  unsigned long expectedLastHProjection = expectedHProjection;
  unsigned long expectedLastVProjection = expectedVProjection;
  if (p._roiArray.dims[0].size % slices != 0) {
    expectedLastHProjection = expectedHProjection + 
      (p._roiArray.dims[0].size % p._roiSlices) * intensity *
      (p._roiArray.dims[1].size / 2);
    expectedLastVProjection = expectedVProjection +
      (p._roiArray.dims[0].size % p._roiSlices) * intensity;
  }

  avg.start();
  CPPUNIT_ASSERT(p.sliceROI());
  avg.end();
  avg.show();

  // Check all the vertical projections (_roiVerticalProjection)
  int slice = 0;
  for (Waveforms::iterator it = p._roiVerticalProjection.begin();
       slice < p._roiSlices; ++it, ++slice) {
    bool zero = false;
    Waveform v = *it;
    unsigned int row = 0;
    for (Waveform::iterator it2 = (*it).begin();
	 row < _image.dims[1].size; ++it2, ++row) {
      if (zero) {
	long expected = expectedVProjection;
	if ((p._roiArray.dims[0].size % slices != 0) &&
	    (slice == p._roiSlices - 1)) {
	  expected = expectedLastVProjection;
	}
	if (*it2 != expected) {
	  std::cout << "INFO: val=" << (*it2) << ", expected=" << expected
		    << " at row=" << row << " slice=" << slice << std::endl;
	}
	CPPUNIT_ASSERT(*it2 == expected);
      }
      else {
	if (*it2 != 0) {
	  std::cout << "INFO: val=" << (*it2) << ", expected=" << 0
		    << " at row=" << row << " slice=" << slice << std::endl;
	}
	CPPUNIT_ASSERT(*it2 == 0);
      }
      zero = !zero;
    }
  }

  // Check the horizontal projection (_roiProjection)
  CPPUNIT_ASSERT(p._roiProjection.size() == p._roiSlices);
  slice = 0;
  for (Waveform::iterator it = p._roiProjection.begin();
       slice < p._roiSlices; ++it, ++slice) {
    unsigned long val = (*it);
    unsigned long expected = expectedHProjection;
    if (slice == p._roiSlices - 1) {
      expected = expectedLastHProjection;
    }
    CPPUNIT_ASSERT(val == expected);
  }

}

void ImageProcessorTest::testLoadBackground() {
  ImageProcessor p;

  CPPUNIT_ASSERT(p.loadBackground("InvalidFile.tif") == false);
  CPPUNIT_ASSERT(p.loadBackground("TestImage-FourLines.tif"));
}

void ImageProcessorTest::testSaveBackground() {
  ImageProcessor p;

  CPPUNIT_ASSERT(loadImage(_singleLineTiff));
  CPPUNIT_ASSERT(p.saveBackground(&_image));

  char *src = (char *) _image.pData;
  char *dest = (char *) p._backgroundArray.pData;

  CPPUNIT_ASSERT(strncmp(src, dest,
			 _image.dims[0].size * _image.dims[1].size * sizeof(epicsUInt16)) == 0);
}

void ImageProcessorTest::testSaveBackgroundFile() {
  ImageProcessor p;

  CPPUNIT_ASSERT(loadImage(_trexTiff2));
  CPPUNIT_ASSERT(p.saveBackground(&_image));
  CPPUNIT_ASSERT(p.saveBackground("/background.tif") == false);
  CPPUNIT_ASSERT(p.saveBackground("./background.tif"));
}

void ImageProcessorTest::testGetSquareAverage() {
  ImageProcessor p;

  CPPUNIT_ASSERT(p.getSquareAverage(0, 0, 100) == 0);

  p._backgroundArray.ndims = 2;

  p._backgroundArray.dims[0].size = 100;
  p._backgroundArray.dims[0].offset = 0;
  p._backgroundArray.dims[0].binning = 0;
  p._backgroundArray.dims[0].reverse = 0;

  p._backgroundArray.dims[1].size = 100;
  p._backgroundArray.dims[1].offset = 0;
  p._backgroundArray.dims[1].binning = 0;
  p._backgroundArray.dims[1].reverse = 0;

  p._backgroundArray.dataType = NDUInt16;
  p._backgroundArray.dataSize = p._backgroundArray.dims[0].size *
    p._backgroundArray.dims[1].size * 4;
  p._backgroundArray.pData = realloc(p._backgroundArray.pData, p._backgroundArray.dataSize);

  memset(p._backgroundArray.pData, 0, p._backgroundArray.dataSize);

  epicsUInt16 *imagePtr = (epicsUInt16 *) p._backgroundArray.pData;
  for (int i = 0; i < 100 * 100; ++i) {
    imagePtr[i] = 100;
  }

  CPPUNIT_ASSERT(p.getSquareAverage(0, 0, 101) == 0);
  CPPUNIT_ASSERT(p.getSquareAverage(0, 0, 101) == 0);
  CPPUNIT_ASSERT(p.getSquareAverage(0, 0, 10) == 100);
}

void ImageProcessorTest::testSaveBackgroundThreshold() {
  ImageProcessor p;

  CPPUNIT_ASSERT(p.saveBackgroundThreshold() == false);

  p._backgroundArray.ndims = 2;

  p._backgroundArray.dims[0].size = 200;
  p._backgroundArray.dims[0].offset = 0;
  p._backgroundArray.dims[0].binning = 0;
  p._backgroundArray.dims[0].reverse = 0;

  p._backgroundArray.dims[1].size = 200;
  p._backgroundArray.dims[1].offset = 0;
  p._backgroundArray.dims[1].binning = 0;
  p._backgroundArray.dims[1].reverse = 0;

  p._backgroundArray.dataType = NDUInt16;
  p._backgroundArray.dataSize = p._backgroundArray.dims[0].size *
    p._backgroundArray.dims[1].size * 4;
  p._backgroundArray.pData = realloc(p._backgroundArray.pData, p._backgroundArray.dataSize);

  memset(p._backgroundArray.pData, 0, p._backgroundArray.dataSize);

  epicsUInt16 *imagePtr = (epicsUInt16 *) p._backgroundArray.pData;

  /**
   * image = 100 ... 100 200 ... 200
   *         100  .           .  200
   *         100   .         .   200
   *         300  .           .  400
   *         300 ... 300 400 ... 400
   */
  for (int i = 0; i < 200; ++i) {
    for (int j = 0; j < 200; ++j) {
      *imagePtr = 100;
      if (i >= 100) {
	*imagePtr += 100;
      }
      if (j >= 100) {
	*imagePtr += 100;
	if (i >= 100) {
	  *imagePtr += 200;
	}
      }
      imagePtr++;
    }
  }

  CPPUNIT_ASSERT(p.saveBackgroundThreshold());
  CPPUNIT_ASSERT(p._backgroundAverageIntensity == (100+200+300+400)/4);
}

void ImageProcessorTest::testSubtractBackground() {
  ImageProcessor p;

  // Create background image, filled with a counter
  p._backgroundArray.ndims = 2;

  p._backgroundArray.dims[0].size = 50;
  p._backgroundArray.dims[0].offset = 0;
  p._backgroundArray.dims[0].binning = 0;
  p._backgroundArray.dims[0].reverse = 0;

  p._backgroundArray.dims[1].size = 50;
  p._backgroundArray.dims[1].offset = 0;
  p._backgroundArray.dims[1].binning = 0;
  p._backgroundArray.dims[1].reverse = 0;

  p._backgroundArray.dataType = NDUInt16;
  p._backgroundArray.dataSize = p._backgroundArray.dims[0].size *
    p._backgroundArray.dims[1].size * 4;
  p._backgroundArray.pData = realloc(p._backgroundArray.pData, p._backgroundArray.dataSize);

  memset(p._backgroundArray.pData, 0, p._backgroundArray.dataSize);

  epicsUInt16 *imagePtr = (epicsUInt16 *) p._backgroundArray.pData;

  for (int i = 0; i < 50; ++i) {
    for (int j = 0; j < 50; ++j) {
      *imagePtr = i * 100 + j;
      imagePtr++;
    }
  }

  // Create ROI image, filled with 10s
  p._roiArray.ndims = 2;

  p._roiArray.dims[0].size = 10;
  p._roiArray.dims[0].offset = 0;
  p._roiArray.dims[0].binning = 0;
  p._roiArray.dims[0].reverse = 0;

  p._roiArray.dims[1].size = 10;
  p._roiArray.dims[1].offset = 0;
  p._roiArray.dims[1].binning = 0;
  p._roiArray.dims[1].reverse = 0;

  p._roiArray.dataType = NDUInt16;
  p._roiArray.dataSize = p._roiArray.dims[0].size *
    p._roiArray.dims[1].size * 4;
  p._roiArray.pData = realloc(p._roiArray.pData, p._roiArray.dataSize);

  epicsUInt16 pixel = 50000;
  epicsUInt16 *roiPtr = (epicsUInt16 *) p._roiArray.pData;
  for (unsigned int i = 0; i < p._roiArray.dims[0].size *
	 p._roiArray.dims[0].size; ++i) {
    roiPtr[i] = pixel + i;
  }

  p._topROI.x = 3;
  p._topROI.y = 8;

  CPPUNIT_ASSERT(p.subtractBackground());
  
  imagePtr = (epicsUInt16 *) p._roiArray.pData;
  
  int k = 0;
  for (unsigned int i = 0; i < p._roiArray.dims[0].size; ++i) {
    for (unsigned int j = 0; j < p._roiArray.dims[1].size; ++j) {
      epicsUInt16 background = ((i + p._topROI.y) * 100 + (j + p._topROI.x));
      epicsUInt16 expected = pixel + k - background;
      if (*imagePtr != expected) {
	std::cout << "ERROR: found " << *imagePtr << ", expected " << expected
		  << " at position [" << i << ", " << j << "]"
		  << std::endl;
	CPPUNIT_ASSERT(false);
      }
      k++;
      imagePtr++;
    }
  }
}

void ImageProcessorTest::testCalculateEnergy() {
  ImageProcessor p;
  /*
  CPPUNIT_ASSERT(p.loadBackground("XTCAV-background.tif"));
  CPPUNIT_ASSERT(loadImage(_xtcavOnTiff));

  p.setROISlices(10);
  p.setExpectedROISize(50, 150);
  p.setCropThreshold(152);
  p.setStrengthParam(55.3928);
  p.setCharge(9.7475e+08);

  CPPUNIT_ASSERT(p.process(&_image, false));

  int h = -1;
  int w = -1;

  Point topLeft;
  Point bottomRight;

  p.getROIPoints(topLeft, bottomRight);

  CPPUNIT_ASSERT(topLeft.x == 171);
  CPPUNIT_ASSERT(topLeft.y == 200);
  CPPUNIT_ASSERT(bottomRight.x == 201);
  CPPUNIT_ASSERT(bottomRight.y == 292);

  p.getROISize(w, h);

  CPPUNIT_ASSERT(w == 31);
  CPPUNIT_ASSERT(h == 93);
  */
  /** Load 50(w)x200(h) test image */
  /* Image has pixel intensities of 0, 655 and 1310 */
  /* Sum of one col is 6550 (6 * 655 + 2 * 1310) */
  CPPUNIT_ASSERT(loadImage(_testEnergyTiff));

  p._topROI.x = 0;
  p._topROI.y = 60;
  p._bottomROI.x = _image.dims[0].size - 1;
  p._bottomROI.y = 70;
  p.saveROI(&_image);
  unsigned short slices = 50; // two cols per slice
  p.setROISlices(slices);

  CPPUNIT_ASSERT(p.sliceROI());

  unsigned slice = 0;
  for (Waveform::iterator it = p._roiProjection.begin();
       slice < p._roiSlices; ++it, ++slice) {
    unsigned long val = (*it);
    unsigned long expected = 6550 * 1;
    if (val != expected) {
      std::cout << "ERROR: expected " << expected << ", found " << val << std::endl;
    }
    CPPUNIT_ASSERT(val == expected);
  }

  slices = 11;
  p.setROISlices(slices);
  CPPUNIT_ASSERT(p.sliceROI());

  p.setCharge(1);
  p.setStrengthParam(2);

  CPPUNIT_ASSERT(p.calculateEnergy());
}