//------------------------------------------------------------------------------
// s u m m a r i z e . c c  -- show the statistics and histogram of an image
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
// Compilation.
//------------------------------------------------------------------------------
// MAC COMPILE:
// g++ `pkg-config --cflags --libs opencv` summarize.cc -o summarize
//
// LINUX COMPILE via Makefile:
/*
# M a k e f i l e    f o r    O p e n C V    v 3 . 1

CC = g++
CFLAGS = -g -I/usr/include/opencv2 -L/usr/lib/
PROGRAM=summarize

(PROGRAM): $(PROGRAM).cc
	$(CC) $(CFLAGS) -o $(PROGRAM) $(PROGRAM).cc \
	-lopencv_highgui -lopencv_core -lm -lopencv_imgproc -lopencv_imgcodecs

# End of Makefile
*/
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
// Boilerplate.
//------------------------------------------------------------------------------

#include <iostream>
#include <fstream>
#include <climits>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace std;
using namespace cv;

const int MAXGREY = 256;
const int MAXDISP = 800;

//------------------------------------------------------------------------------
// Routines.
//------------------------------------------------------------------------------
// See http://stackoverflow.com/questions/10167534/how-to-find-out-what-type-of-a-mat-object-is-with-mattype-in-opencv

string type2str (int type)
{
  string r;

  uchar depth = type & CV_MAT_DEPTH_MASK;
  uchar chans = 1 + (type >> CV_CN_SHIFT);

  switch (depth) {
  case CV_8U:  r = "8U"; break;
  case CV_8S:  r = "8S"; break;
  case CV_16U: r = "16U"; break;
  case CV_16S: r = "16S"; break;
  case CV_32S: r = "32S"; break;
  case CV_32F: r = "32F"; break;
  case CV_64F: r = "64F"; break;
  default:     r = "User"; break;
  }
  r += "C";
  r += (chans+'0');  // requires < 10 channels!

  return r;
}


// immin -- return the minimum value in an image
int
immin (Mat im)
{
  int ny, nx, nc, y, x, c, v, lo = INT_MAX;
  Vec3b pixel;
  
  ny = im.rows;
  nx = im.cols;
  nc = im.channels();
 
  for (y = 0; y < ny; y++) {
    for (x = 0; x < nx; x++) {
      pixel = im.at<Vec3b>(y, x);
      for (c = 0; c < nc; c++) {
        v = int (pixel[c]);
        if (v < lo) lo = v;
      }
    }
  }
  return lo;
}


// immax -- return the maximum value in an image
int
immax (Mat im)
{
  int ny, nx, nc, y, x, c, v, hi = 0;
  Vec3b pixel;
  
  ny = im.rows;
  nx = im.cols;
  nc = im.channels();

  for (y = 0; y < ny; y++) {
    for (x = 0; x < nx; x++) {
      pixel = im.at<Vec3b>(y, x);
      for (c = 0; c < nc; c++) {
        v = int (pixel[c]);
        if (v > hi) hi = v;
      }
    }
  }
  return hi;
}

// ave -- return the mean of all pixels of an image
double
ave (Mat im)
{
  int ny, nx, nc, y, x, c, v;
  double total = 0.0;
  Vec3b pixel;
  
  ny = im.rows;
  nx = im.cols;
  nc = im.channels();

  for (y = 0; y < ny; y++) {
    for (x = 0; x < nx; x++) {
      pixel = im.at<Vec3b>(y, x);
      for (c = 0; c < nc; c++) {
        v = int (pixel[c]);
        total += v;
      }
    }
  }
  return total / ny / nx / nc;
}

// sdev -- return the standard deviation of all pixels of an image
double
sdev (Mat im)
{
  int ny, nx, nc, y, x, c;
  double v, total = 0.0;
  Vec3b pixel;
  
  double mean = ave (im);

  ny = im.rows;
  nx = im.cols;
  nc = im.channels();

  for (y = 0; y < ny; y++) {
    for (x = 0; x < nx; x++) {
      pixel = im.at<Vec3b>(y, x);
      for (c = 0; c < nc; c++) {
        v = double (pixel[c]) - mean;
        total += v * v;
      }
    }
  }
  return sqrt (total / ny / nx / nc);
}

// histogram -- work out the histogram of an image
void
histogram (Mat im, int *ab, int *h)
{
  int ny, nx, nc, y, x, c, i, v;
  Vec3b pixel;
  
  // Fill up the abscissa and zero the histogram.
  for (i = 0; i < MAXGREY; i++) {
    ab[i] = i;
    h[i] = 0;
  }

  // Work out the histogram from the image values.
  ny = im.rows;
  nx = im.cols;
  nc = im.channels();

  for (y = 0; y < ny; y++) {
    for (x = 0; x < nx; x++) {
      pixel = im.at<Vec3b>(y, x);
      for (c = 0; c < nc; c++) {
        v = int (pixel[c]);
        h[v] += 1;
      }
    }
  }

  // Plotting graphs directly from C++ code is pesky, so generate output
  // that can simply be given to Gnuplot.
  ofstream of ("histogram.dat", ios::out);
  for (i = 0; i < MAXGREY; i++)
    of << ab[i] << " " << h[i] << endl;  
  of.close();
}

// iround -- round a value > 0 to the nearest integer
int
iround (double v)
{
  return int (v + 0.5);
}

//------------------------------------------------------------------------------
// The main program.
//------------------------------------------------------------------------------

int main (int argc, char *argv[])
{
  char *fn;
  int ny, nx, nc;
  int lo, hi;
  double a;
  Mat im;
  Scalar chave;

  // Ensure we were invoked correctly.
  if (argc < 2) {
    std::cerr << "Usage:" << argv[0] << " <image>..." << std::endl;
    return 1;
  }

  // Process the files given on the command line.
  for (int ifn = 1; ifn < argc; ifn++) {
    // Read in the image and print out its dimensions.
    fn = argv[ifn];
    cout << fn << ":" << endl;
    im = imread (fn);
    ny = im.rows;
    nx = im.cols;
    nc = im.channels();
    cout << "  Dimensions: " << nx << " pixels, " << ny << " lines, "
         << nc << " channels of type " << type2str (im.type()) << endl;
   
    // Calculate and output some important statistics.
    lo = immin (im);
    hi = immax (im);
    cout << "  Range: " << lo << " to " << hi << endl;
    cout << "  Mean: " << ave (im) << " (using ave)" << endl;

    // OpenCV's mean routine returns _channel_ means, not the image mean.
    chave = mean (im);
    a = (chave[0] + chave[1] + chave[2]) / 3.0;
    cout << "  Mean: " << a << " (using OpenCV's mean)" << endl;
    cout << "  Standard deviation: " << sdev (im) << endl;

    // Work out and display the histogram.
    int x[MAXGREY], h[MAXGREY];
    histogram (im, x, h);

    // To display the image, we first create the window.
    namedWindow (fn, CV_WINDOW_AUTOSIZE);

    // Ensure the image is no more than maxdisp pixels in x or y.
    if (ny > MAXDISP or nx > MAXDISP) {
      int nmax = max (ny, nx);
      double fac = double (MAXDISP) / nmax;
      int nny = iround (ny * fac);
      int nnx = iround (nx * fac);
      Mat disp_im;
      cout <<  "  [re-sizing to " << nnx << " x " << nny <<
        " pixels for display]" << endl;
      resize (im, disp_im, Size(nnx, nny));
      imshow (fn, disp_im);
    } else {
      imshow (fn, im);
    }

    // Wait for the user to do something.
    waitKey (0);
    destroyWindow (fn);
    cout << endl;
  }
}

//------------------------------------------------------------------------------
// End of summarize.cc
//------------------------------------------------------------------------------