#include <math.h>
#include <string.h>

#include <itk/raster.h>
#include <itk/rotate.h>
#include <itk/assign.h>

#ifndef DO_ITK_INLINE
#include <itk/rotate.inl>
#endif

bool priv_rotate_1 ( const Raster &src, Raster &dest, double angle);
bool priv_rotate_4 ( const Raster &src, Raster &dest, double angle);
bool priv_rotate_8 ( const Raster &src, Raster &dest, double angle);
bool priv_rotate_16( const Raster &src, Raster &dest, double angle);
bool priv_rotate_32( const Raster &src, Raster &dest, double angle);


EMPTY_CONSTRUCTOR( Rotate, Transformation)
EMPTY_DESTRUCTOR( Rotate)


Rotate::Rotate( double new_angle)
: Transformation()
, m_angle( new_angle)
{
}


bool
Rotate::do_rotate( const Raster &src, Raster &dest)
{
  return priv_rotate_1( src, dest, angle());
}


bool
priv_rotate_1( const Raster &src, Raster &dest, double angle)
{
  double radians      = -angle / (180.0 / M_PI);
  double cosval       = cos( radians);
  double sinval       = sin( radians);
  size_t old_width    = src.width ();
  size_t old_height   = src.height();
  size_t new_width    = static_cast<size_t>( fabs( old_width  * cosval)
                                           + fabs( old_height * sinval));
  size_t new_height   = static_cast<size_t>( fabs( old_height * cosval)
                                           + fabs( old_width  * sinval));
  size_t half_owidth  = old_width  / 2;
  size_t half_oheight = old_height / 2;
  size_t half_nwidth  = new_width  / 2;
  size_t half_nheight = new_height / 2;

  Raster temp_raster;
  Raster &work = (&src == &dest) ? temp_raster : dest;

  work.redimension( new_width, new_height, src.depth(), src.interp());
  memset( work.raster(), 0, work.data_size());

  Raster_Cell *sdata  = src.raster();
  Raster_Cell *ddata  = work.raster();
  Dimension    srsize = src.row_size();
  Dimension    drsize = work.row_size();
  Raster_Cell c;

  double raw_x0 = (0.0 - half_nwidth) * cosval
                + (0.0 - half_nheight) * sinval
                + half_owidth;
  double raw_y0 = (half_nwidth - 0.0) * sinval
                + (0.0 - half_nheight) * cosval
                + half_oheight;
  double raw_x1 = (1.0 - half_nwidth) * cosval
                + (0.0 - half_nheight) * sinval
                + half_owidth;
  double raw_y1 = (half_nwidth - 1.0) * sinval
                + (0.0 - half_nheight) * cosval
                + half_oheight;
  double raw_x2 = (0.0 - half_nwidth) * cosval
                + (1.0 - half_nheight) * sinval
                + half_owidth;
  double raw_y2 = (half_nwidth - 0.0) * sinval
                + (1.0 - half_nheight) * cosval
                + half_oheight;
  long   cdx = static_cast<long>( 65536 * (raw_x1 - raw_x0));
  long   cdy = static_cast<long>( 65536 * (raw_y1 - raw_y0));
  long   rdx = static_cast<long>( 65536 * (raw_x2 - raw_x0));
  long   rdy = static_cast<long>( 65536 * (raw_y2 - raw_y0));

  long   rsx = static_cast<long>( 65536 * (raw_x0));
  long   rsy = static_cast<long>( 65536 * (raw_y0));
  long   x;
  long   y;
  long   xx;
  long   yy;
  long   c3;
  Raster_Data s1, s2, s3, s4;
  unsigned char r, g, b;
  long     mant_x, mant_y, nmant_x, nmant_y;

  for (int row = 0; row < new_height; ++row) {
    x = rsx;
    y = rsy;
    for (int col = 0; col < new_width; ++col) {
      xx = (x >> 16);
      yy = (y >> 16);
      if ((xx >= 0) && (xx < old_width) &&
          (yy >= 0) && (yy < old_height)) {
        c = sdata[ yy * srsize + (xx >> 3)];
        if ((c >> (7 - (xx & 0x07))) & 0x01) {
          ddata[ row * drsize + (col >> 3)] |= (1 << (7 - (col & 0x07)));
        }
      }
      x += cdx;
      y += cdy;
    }
    rsx += rdx;
    rsy += rdy;
  }
  
  Assign().transform( work, dest);

  return true;
}