/*****************************************************************************
 * ass1log.cpp
 ***************************************************************************** 
 *
 * Description: Assignment 1 (COMPSCI775)
 *
 * Topic: Edge Detection in Colour Images - Laplacian of Gaussian
 *
 * (c) 2002 Christian Graf, Uli Schroeder (Group 7)
 *
 *****************************************************************************
 *
 * $Revision: 1.0 $
 *
 */

#include <math.h>
#include <iostream.h>
#include "HalconCpp.h"

class LoG {

  public:
  
    LoG(HByteImage inImg, HByteImage *outImg, double aTheta, int aKernelSize) {
      // Store object attributes
      kernelSize = aKernelSize;
	  offset = (int)floor(kernelSize/2);
      imgWidth = inImg.Width();
	  imgHeight = inImg.Height();
      theta = aTheta;
      // Allocate memory for kernel matrix
      kernelMat = (double **) new double*[kernelSize];
      for (int i=0; i < kernelSize; i++) {
        kernelMat[i] = (double *) new double[kernelSize];
      }
      // Generate kernel matrix
      genKernel();
      // Apply kernel to source picture
      applyKernel(inImg, *outImg);
    }
    
    ~LoG() {
      // Free memory allocated for kernel matrix
      for (int i = 0; i<kernelSize; i++) {
        delete [] kernelMat[i];
      }
      delete [] kernelMat;
    }

  private:
  
    double **kernelMat, theta;
    int kernelSize, offset, imgWidth, imgHeight;

    void genKernel() {
      for(int j=0; j<kernelSize; ++j){
        for(int i=0; i<kernelSize; ++i){
	      int x = i-offset;
	      int y = j-offset;
	      kernelMat[i][j] = logElement(x,y);
        }
      }
    }

    double logElement(int x, int y) {
      double g = 0;
      double s = -((x*x)+(y*y))/(2*theta*theta);
      g = -(1/(PI*pow(theta,4)))*(1+s)*exp(s);
      return g;
    }

    void applyKernel(HByteImage inImg, HByteImage outImg) {
      // Store local parameters
      double pixsum = 0;
      // Allocate memory for picBufferIn and picBufferOut
      int **picBufferOut = NULL, **picBufferIn = NULL;
      picBufferIn = (int **) new int*[imgWidth];
      for (int i=0; i<imgWidth; i++) {
        picBufferIn[i] = (int *) new int[imgHeight];
      }
      picBufferOut = (int **) new int*[imgWidth];
      for (int i=0; i<imgWidth; i++) {
        picBufferOut[i] = (int *) new int[imgHeight];
      }
      // Copy input picture to picBufferIn
      for (int y=0; y<imgHeight; y++) {
        for (int x=0; x<imgWidth; x++) {
	        picBufferIn[x][y] = (int)inImg.GetPixVal(x, y);
        }
      }
	  // Process convolution
	  double max = 0, min = 0;
	  for (int y=offset; y<(imgHeight-offset); y++) {
        for (int x=offset; x<(imgWidth-offset); x++) {
          pixsum = 0;
          // Apply kernel to current position
          for (int k_y = -offset; k_y<=offset; k_y++) {
            for (int k_x = -offset; k_x<=offset; k_x++) {
		      pixsum += picBufferIn[x+k_x][y+k_y]*kernelMat[offset+k_x][offset+k_y];
            }          
          }
          picBufferOut[x][y] = (int)floor(pixsum+0.5);
		  if (pixsum > max) max = pixsum; 
		  if (pixsum < min) min = pixsum;
        }
      }
/*
      // Scale values to the range 0-255
      double factor = 255/(max-min);
      for (int y=0; y<picBufferOut_h; y++) {
        for (int x=0; x<picBufferOut_w; x++) {
          picBufferOut[x][y]=(int)floor(((picBufferOut[x][y]-min)*factor)+0.5);
        }
      }
*/
      // Check for zero crossings
      for (int y=offset; y<(imgHeight-offset); y++) {
        for (int x=offset; x<(imgWidth-offset); x++) {
          if (((picBufferOut[x][y]*picBufferOut[x+1][y])<0) ||
              ((picBufferOut[x][y]*picBufferOut[x][y+1])<0)) {
	        picBufferOut[x][y] = 0;
          } else {
	        picBufferOut[x][y] = 255;
          }
        }
      }
      // Copy picBufferOut to output picture
      for (int y=offset; y<(imgHeight-offset); y++) {
        for (int x=offset; x<(imgWidth-offset); x++) {
	        outImg.SetPixVal(x, y, picBufferOut[x][y]);
        }
      }
      // Free memory allocated for picBufferIn and picBufferOut
      for (int i = 0; i<imgWidth; i++) {
        delete [] picBufferIn[i];
      }
      delete [] picBufferIn;
      for (int i = 0; i<imgWidth; i++) {
        delete [] picBufferOut[i];
      }
      delete [] picBufferOut;
    }

};

main(int argc, char* argv[]) {
  cout << "---------------------------------" << endl;
  cout << " Assignment 1 (COMPSCI775):      " << endl;
  cout << " Edge Detection in Colour Images " << endl;
  cout << " (Laplacian of Gaussian)         " << endl;
  cout << "---------------------------------" << endl;
  cout << " Group 7:                        " << endl;
  cout << " Christian Graf                  " << endl;
  cout << " Uli Schroeder                   " << endl;
  cout << "---------------------------------" << endl << endl;
  if (argc != 5) {
    cout << " Wrong number of arguments!" << endl << endl;
    cout << " Program usage:" << endl;
    cout << " ./ass1log inputimage outputimage theta kernelsize" << endl << endl;
    return(-1);
  }
  cout << " Processing image " << argv[1] << "..." << endl;
  HImage     myImage(argv[1]);
  HByteImage r, g, b, r_edges, g_edges, b_edges, result;
  r = myImage.Decompose3(&g, &b);
  g = myImage.Decompose3(&r, &b);
  b = myImage.Decompose3(&r, &g);
  r_edges = HByteImage(myImage.Width(),myImage.Height());
  g_edges = HByteImage(myImage.Width(),myImage.Height());
  b_edges = HByteImage(myImage.Width(),myImage.Height());
  result = HByteImage(myImage.Width(),myImage.Height());
  // Apply Laplacian of Gaussian
  new LoG(r, &r_edges, atof(argv[3]), atoi(argv[4]));
  new LoG(g, &g_edges, atof(argv[3]), atoi(argv[4]));
  new LoG(b, &b_edges, atof(argv[3]), atoi(argv[4]));
  // Compose  result (An edge has to be in at least two pictures.)
  for (int x=0; x<result.Width(); x++) {
    for (int y=0; y<result.Height(); y++) {
      if (r_edges(x,y)+g_edges(x,y)+b_edges(x,y)<256) {
        result(x,y) = 0;
      } else {
        result(x,y) = 255;
      }
    }
  }
  cout << " Image processing finished..." << endl;
  // Write result to file
  result.WriteImage("jpg 50", 0xffffff, argv[2]);
/*
  // Display result
  HWindow img1;
  result.Display (img1);
  img1.Click();
*/
  cout << " The End!" << endl << endl;
  return(0);
}