#include <math.h>
#include <common/fixed.h>
#include <ocr/get_skew.h>
#include <itk/raster.h>
#include <itk/bit_count.h>
#include <itk/rotate.h>

#ifndef DO_OCR_INLINE
#include <ocr/get_skew.inl>
#endif

inline double sqr( double x) { return x * x; }
inline double abs( double x) { return (x > 0) ? x : -x; }

double min( double d1, double d2) { return (d1 < d2) ? d1 : d2; }
int    min( int    i1, int    i2) { return (i1 < i2) ? i1 : i2; }
double max( double d1, double d2) { return (d1 > d2) ? d1 : d2; }
int    max( int    i1, int    i2) { return (i1 > i2) ? i1 : i2; }


Get_Skew::Get_Skew( unsigned new_max_angle, int new_sample_skip)
: m_max_angle  ( new_max_angle)
, m_sample_skip( new_sample_skip)
, m_skew       ( 0.0)
, m_max_rows   ( 0)
, m_rows       ( 0)
{
}


Get_Skew::~Get_Skew()
{
  delete []m_rows;
}


double
Get_Skew::get_skew( const Raster & image)
{
  double var_test;
  double var_opt     = 0.0;
  Rotate rotate;
  double step        = 10.0;
  double angle       = 0.0;
  double opt         = 0.0;
  int    i;
  int    j;
  int    iter        = 0;
  int    low         = -25;
  int    high        = 25;

  for (i = 0; i < 2; ++i) {
    angle = opt + low * step;
    for (j = low; j <= high; ++j) {
      var_test = histogram( image, angle/10);
      if (var_test > var_opt) {
        var_opt = var_test;
        opt    = angle;
      } else {
        /* do nothing - I'd like to short-circuit, but it doesn't work :( */
      }
      angle += step;
    }
    step /= 10.0;
    low  = -10;
    high = 10;
  }
    
  m_skew = opt / 10.0;
  return m_skew;
}


unsigned
Get_Skew::cast_ray( const Raster &raster, int row, const Fixed &dx, int dy)
{
  unsigned  bits  = 0;
  Fixed     start = 0;
  Fixed     end   = 0;
  Fixed     width = static_cast<int>( raster.width());
  int       r     = row;

  while (start < width) {
    end = start + dx;
    if (end > width) {
      end = width;
    }
    bits += bit_count( raster[ r]
                     , static_cast<int>( start)
                     , static_cast<int>( end));
    start = end;
    r += dy;
  }
  
  return bits;
}


double
Get_Skew::histogram( const Raster &raster, double angle)
{
  int    i;
  double sum;
  double mean;
  double angle_diff  = tan(angle / (180.0 / M_PI));
  int    diff_y      = -static_cast<int>( raster.width() * angle_diff);
  int    min_y       = max( 0, diff_y);
  int    max_y       = min( static_cast<int>( raster.height())
                          , raster.height() + diff_y);
  int    num_rows;
  Fixed  dx;
  int    dy;


  if (raster.height() > m_max_rows) {
    delete []m_rows;
    m_rows = new unsigned[ raster.height()];
    m_max_rows = raster.height();
  }
  num_rows = (max_y - min_y) / m_sample_skip + 1;

  if (angle < 0) {
    dy = -1;
  } else {
    dy = +1;
  }
  if ((-0.05 < angle) && (angle < 0.05)) {
    dx = static_cast<int>( raster.width());
  } else {
    dx = dy / (tan( angle / (180.0 / M_PI)));
  }

  for (i = 0; i < num_rows; ++i) {
    m_rows[ i] = cast_ray( raster, min_y + i * m_sample_skip, dx, dy);
  }

  sum = 0.0;
  for (i = 0; i < num_rows; ++i) {
    sum += m_rows[ i];
  }
  mean = sum / num_rows;

  sum = 0.0;
  for (i = 0; i < num_rows; ++i) {
    sum += sqr( m_rows[ i] - mean);
  }

  return sum / num_rows;
}