rotate-me-fast/rotation.cpp

#include <string>
#include <fstream>
#include <iostream>
#include <sstream>
#include <cmath>
#include <cassert>
#include <cstring>

using namespace std;


template <typename T>
struct TPoint {
  T x;
  T y;

  TPoint(T a, T b)
  : x(a)
  , y(b)
  {}
};
typedef TPoint<unsigned int> Point;
typedef TPoint<int> APoint; // absolute point, can be negative
typedef TPoint<double> DPoint; // absolute point, can be negative
template<typename Elem, typename Traits, typename T>
std::basic_ostream<Elem, Traits>& operator << (std::basic_ostream<Elem, Traits>& o, TPoint<T> const& p)
{
  o << "(" << p.x << ", " << p.y << ")";
  return o;
}

struct Image {
  unsigned int width;
  unsigned int height;
  uint8_t* r_chan;
  uint8_t* g_chan;
  uint8_t* b_chan;

  Image()
  : width(0)
  , height(0)
  , r_chan(NULL)
  , g_chan(NULL)
  , b_chan(NULL)
  {}

  Image(unsigned int w, unsigned int h)
  {
    this->width = w;
    this->height = h;
    r_chan = new uint8_t[width * height];
    memset(r_chan, 0, width * height * sizeof (uint8_t));
    g_chan = new uint8_t[width * height];
    memset(g_chan, 0, width * height * sizeof (uint8_t));
    b_chan = new uint8_t[width * height];
    memset(b_chan, 0, width * height * sizeof (uint8_t));
  }

  Image(string const& path)
  : Image()
  {
    ifstream is(path);
    if (!is.is_open())
    {
      cerr << "Cannot open file '" << path << "'" << endl;
      abort();
    }

    if (!this->read_header(is))
    {
      cerr << "Invalid header." << endl;
      abort();
    }

    if (!this->read_body(is))
    {
      delete r_chan;
      r_chan = nullptr;
      delete g_chan;
      r_chan = nullptr;
      delete b_chan;
      r_chan = nullptr;

      cerr << "Invalid header." << endl;
      abort();
    }
  }

  bool save(string const& path)
  {
    ofstream os(path);
    if (!os.is_open())
    {
      cerr << "Cannot open file '" << path << "'" << endl;
      return false;
    }
    this->write_header(os);
    this->write_body(os);
    return true;
  }

  void set_pixel(unsigned int x, unsigned int y, uint8_t r, uint8_t g, uint8_t b)
  {
    if (x >= width || y >= height)
    {
        cerr << "Point (" << x << ", " << y << ") out of bounds" << endl;
        cerr << " Image dimensions: " << width << " x " << height << endl;
        assert(false);
    }
    int const index = y * width + x;
    r_chan[index] = r;
    g_chan[index] = g;
    b_chan[index] = b;
  }

  void set_pixel(Point const& p, uint8_t r, uint8_t g, uint8_t b)
  {
    this->set_pixel(p.x, p.y, r, g, b);
  }

  void get_pixel(APoint const& p, uint8_t& r, uint8_t& g, uint8_t& b) const
  {
    if (p.x < 0 || p.x >= (int) width || p.y < 0 || p.y >= (int) height)
    {
      // set out of domain pixels to black
      r = 0;
      g = 0;
      b = 0;
      return;
    }

    unsigned int const index = p.y * width + p.x;
    r = r_chan[index];
    g = g_chan[index];
    b = b_chan[index];
  }

  private:
    bool read_header(std::ifstream& istr)
    {
      // check magic
      if (istr.get() != 'P' )
      {
        return false;
      }

      char type = static_cast<char>(istr.get());
      if (type != '6')
      {
        return false;
      }

      if (istr.get() != '\n')
      {
        return false;
      }

      // skip comments
      while (istr.peek() == '#')
      {
        std::string line;
        std::getline(istr, line);
      }

      // get size
      istr >> width >> height;
      if (width == 0 || height == 0)
      {
        return false;
      }

      // get maxvalue
      if (istr.get() != '\n')
      {
        return false;
      }

      int max_value = -1;
      istr >> max_value;
      if (max_value > 255)
      {
        return false;
      }

      if (istr.get() != '\n')
      {
        return false;
      }

      // cout << "width: " << width << endl;
      // cout << "height: " << height << endl;

      return true;
    }

    bool write_header(std::ofstream& ostr)
    {
      ostr << "P6" << endl;
      ostr << width << " " << height << endl;
      ostr << "255" << endl;
      return true;
    }

    bool read_body(std::ifstream& istr)
    {
      r_chan = new uint8_t[width * height];
      g_chan = new uint8_t[width * height];
      b_chan = new uint8_t[width * height];

      for (unsigned int row = 0; row < height; ++row)
      {
        for (unsigned int col = 0; col < width; ++col)
        {
          int index = row * width + col;
          r_chan[index] = istr.get();
          g_chan[index] = istr.get();
          b_chan[index] = istr.get();
        }
      }

      return true;
    }

    bool write_body(std::ofstream& ostr)
    {
      for (unsigned int row = 0; row < height; ++row)
      {
        for (unsigned int col = 0; col < width; ++col)
        {
          int index = row * width + col;
          ostr << (char) r_chan[index];
          ostr << (char) g_chan[index];
          ostr << (char) b_chan[index];
        }
      }
      return true;
    }
};



//
//
// Trigonometry
//

DPoint convert_grid_coord(Image const& img, Point const& p)
{
  return DPoint(p.x - img.width / 2.0f + 0.5, p.y - img.height / 2.0f + 0.5);
}

double convert_radian(Image const& img, Point const& p, double const ratio)
{
  DPoint centered = convert_grid_coord(img, p);
  cout << "-> grid " << centered << endl;
  double const cos_value = centered.x * ratio;
  cout << "cos = " << cos_value << endl;
  double const sin_value = - (centered.y * ratio);
  cout << "sin = " << sin_value << endl;
  double angle = acos(cos_value);
  if (sin_value < 0)
  {
    cout << "mirror angle" << endl;
    angle = (2 * M_PI) - angle;
  }
  cout << "radian = " << angle << endl;
  cout << "acos = " << acos(cos_value) << endl;
  cout << "asin = " << asin(sin_value) << endl;

  cout << "revert cos = " << cos(angle) << endl;
  cout << "revert sin = " << sin(angle) << endl;
  cout << "full cycle sin: " << sin(asin(sin_value)) << endl;

  return angle;
}

DPoint convert_abs_coord(double const angle, double const ratio)
{
  return DPoint(cos(angle) / ratio, sin(angle) / ratio);
}

APoint convert_img_coord(Image const& img, DPoint const& p)
{
  int x = round(p.x + (img.width / 2.0f) - 0.5);
  int y = round(p.y + (img.height / 2.0f) - 0.5);
  APoint p_mapped(x, y);
//  if (p_mapped.x >= img.width || p_mapped.y >= img.height)
//  {
//    cerr << "Point coord mapping" << endl;
//    cerr << " Input point:" << p << endl;
//    cerr << " Point " << p_mapped << " out of bounds" << endl;
//    cerr << " Image dimensions: " << img.width << " x " << img.height << endl;
//    assert(false);
//  }

  return p_mapped;
}

double compute_ratio(Image const& img)
{
  double const trigo_length = (sqrt(img.width * img.width + img.height * img.height) - 1) / 2;
  return 1.0f / trigo_length;
}

inline
bool fequal(float a, float b, float sigma)
{
  return abs(a - b) < sigma;
}

void compute_output_size(Image const& src, double const rotation, unsigned int& width, unsigned int& height)
{
  double const ratio = compute_ratio(src);
  double min_w = 0;
  double max_w = 0;
  double min_h = 0;
  double max_h = 0;

  //cout << "Image dimensions: " << src.width << " x " << src.height << endl;
  Point p(0, 0);
  double angle = convert_radian(src, p, ratio);
  DPoint tl = convert_abs_coord(angle + rotation, ratio);
  min_w = min(min_w, tl.x);
  max_w = max(max_w, tl.x);
  min_h = min(min_h, tl.y);
  max_h = max(max_h, tl.y);
  // debug print
  if (rotation == 0.0)
  {
    cout << "Rotated " << p << "  = " << tl << endl << endl;
  }

  p = Point(src.width - 1, 0);
  angle = convert_radian(src, p, ratio);
  DPoint tr = convert_abs_coord(angle + rotation, ratio);
  min_w = min(min_w, tr.x);
  max_w = max(max_w, tr.x);
  min_h = min(min_h, tr.y);
  max_h = max(max_h, tr.y);
  if (rotation == 0.0)
  {
    cout << "Rotated " << p << "  = " << tr << endl << endl;
  }

  p = Point(0, src.height - 1);
  angle = convert_radian(src, p, ratio);
  DPoint bl = convert_abs_coord(angle + rotation, ratio);
  min_w = min(min_w, bl.x);
  max_w = max(max_w, bl.x);
  min_h = min(min_h, bl.y);
  max_h = max(max_h, bl.y);
  if (rotation == 0.0)
  {
    cout << "Rotated " << p << "  = " << bl << endl << endl;
  }

  p = Point(src.width - 1, src.height - 1);
  angle = convert_radian(src, p, ratio);
  DPoint br = convert_abs_coord(angle + rotation, ratio);
  min_w = min(min_w, br.x);
  max_w = max(max_w, br.x);
  min_h = min(min_h, br.y);
  max_h = max(max_h, br.y);
  if (rotation == 0.0)
  {
    cout << "Rotated " << p << "  = " << br << endl << endl;
  }

  width = (int) (max_w - min_w) + 1;
  height = (int) (max_h - min_h) + 1;
}

void convert_abs_to_polar_coord(DPoint const& p, double const ratio, double& angle, double& dist)
{
  double const cos_value = p.x * ratio;
  double const sin_value = - (p.y * ratio);
  angle = acos(cos_value);
  if (sin_value < 0)
    angle = 2 * M_PI - angle;

  dist = sqrt(p.x * p.x + p.y * p.y);
}

DPoint convert_polar_to_grid_coord(double const angle, double const distance)
{
  return DPoint(cos(angle) * distance, - (sin(angle) * distance));
}


//
//
// Point rotation
//

APoint rotate(Image const& src, Point const& p, double const ratio, double const rotation, Image const& rotated)
{
  double angle = convert_radian(src, p, ratio);
  DPoint a_point = convert_abs_coord(angle + rotation, ratio);
  return convert_img_coord(rotated, a_point);
}

//
//
// Drawing
//

void draw_line(Image& img, unsigned int x1, unsigned int y1, unsigned int x2, unsigned int y2)
{
  int x_inc = x1 <= x2 ? 1 : -1;

  unsigned int const y_min = min(y1, y2);
  unsigned int const y_max = max(y1, y2);
  double slope = (double) y2 - y1;
  if (x1 != x2)
    slope = ((double) y2 - y1) / abs((double) x2 - x1);
  int y_inc = slope > 0 ? 1 : -1;

  if (x1 == x2)
  {
    for (unsigned int runner = y1; runner != y2; runner+= y_inc)
      img.set_pixel(x1, runner, 255, 0, 0); // set line to red
    return;
  }

  if (y1 == y2)
  {
    for (unsigned int runner = x1; runner != x2; runner+= x_inc)
      img.set_pixel(runner, y1, 255, 0, 0); // set line to red
    return;
  }

  unsigned int previous_y = y1;
  for (unsigned int i = x1, steps = 0; i != x2; i += x_inc, ++steps)
  {
    unsigned int y = slope * steps + y1;
    y = min(y, y_max);
    y = max(y, y_min);
    for (unsigned int runner = previous_y; runner != y; runner+= y_inc)
      img.set_pixel(i, runner, 255, 0, 0); // set line to red
    previous_y = y;
  }
}

void draw_line(Image& img, APoint const& p1, APoint const& p2)
{
  draw_line(img, p1.x, p1.y, p2.x, p2.y);
}

void draw_outline(Image const& input, unsigned int degrees, string const& name)
{
  double const rotation = (degrees / 180.0f) * M_PI;
  unsigned int w = 0;
  unsigned int h = 0;
  compute_output_size(input, rotation, w, h);
  cout << "rotation(" << degrees << ") -> " << w << " x " << h << endl;
  Image rotated(w, h);

  double const ratio = compute_ratio(input);
  APoint tl = rotate(input, Point(0, 0), ratio, rotation, rotated);
  APoint tr = rotate(input, Point(input.width - 1, 0), ratio, rotation, rotated);
  APoint bl = rotate(input, Point(0, input.height - 1), ratio, rotation, rotated);
  APoint br = rotate(input, Point(input.width - 1, input.height - 1), ratio, rotation, rotated);
  cout << tl << " " << tr << " " << bl << " " << br << endl;

  draw_line(rotated, tl, tr);
  draw_line(rotated, tr, br);
  draw_line(rotated, br, bl);
  draw_line(rotated, bl, tl);

  stringstream ss;
  ss << "check_lines_" << name << "_" << degrees << ".ppm";
  rotated.save(ss.str());
}



//
//
// Image rotation
//

Image rotate(Image const& src, double angle)
{
  double const rotation = (angle / 180.0f) * M_PI;
  unsigned int w = 0;
  unsigned int h = 0;
  compute_output_size(src, rotation, w, h);
  Image rotated(w, h);

  // debug print
  if (rotation == 0.0)
  {
    cout << "src dimensions: " << src.width << " x " << src.height << endl;
    cout << "rotated dimensions: " << w << " x " << h << endl;
  }

  double const ratio = compute_ratio(src);

  for (int y = 0; y < (int) rotated.height; ++y)
  {
    for (int x = 0; x < (int) rotated.width; ++x)
    {
      Point const p(x, y);
      DPoint const d = convert_grid_coord(rotated, p);
      double p_angle = 0;
      double dist = 0;
      convert_abs_to_polar_coord(d, ratio, p_angle, dist);
      DPoint const src_rotated_point = convert_polar_to_grid_coord(p_angle - rotation, dist);

      // FIXME: get source points
      APoint src_p = convert_img_coord(src, src_rotated_point);
      uint8_t r = 0;
      uint8_t g = 0;
      uint8_t b = 0;
      src.get_pixel(src_p, r, g, b);
      rotated.set_pixel(p, r, g, b);
    }
  }

  return rotated;
}



//
//
// Check
//

bool check_points()
{
  Image five(5, 5);
  Point origin(0, 0);
  DPoint d1 = convert_grid_coord(five, origin);
  assert(d1.x == -2);
  assert(d1.y == -2);

  return true;
}

bool check_trigo()
{
  Image square(500, 500);
  double const ratio = compute_ratio(square);
  cout << "ratio: " << ratio << endl;

  // Check that the origin of a square image is at sqrt(2) / 2
  double const angle = convert_radian(square, Point(0, 0), ratio);

  double const sigma = 1.0e-2;
  if (!fequal(angle, 3 * M_PI / 4, sigma))
  {
    cout << __LINE__ << " | Invalid angle value: " << angle << " != " << 3 * M_PI / 4 << endl;
    return false;
  }

  // Check that we can reverse the origin point.
  DPoint const abs_reverse_point = convert_abs_coord(angle, ratio);
  cout << "reversed abs origin: " << abs_reverse_point << endl;
  APoint const reverse_point = convert_img_coord(square, abs_reverse_point);
  cout << "reversed origin in square: " << reverse_point << endl;
  if (abs(0.0 - reverse_point.x) > sigma)
  {
    cerr << "Reverse origin:" << endl;
    cout << "Invalid x value: " << reverse_point.x << " != " << 0 << endl;
    return false;
  }
  if (abs(0.0 - reverse_point.y) > sigma)
  {
    cerr << "Reverse origin:" << endl;
    cout << "Invalid y value: " << reverse_point.y << " != " << 0 << endl;
    return false;
  }

  // Check that when rotating the origin by 45 degrees
  double const rotation = M_PI / 4; // 45 degrees
  unsigned int w = 0;
  unsigned int h = 0;
  compute_output_size(square, rotation, w, h);

  // failed check: is precision an issue?
  if (true)
  {
    if (!fequal(w, square.width * sqrt(2), sigma * square.width)
        || !fequal(h, square.height * sqrt(2), sigma * square.height))
    {
      cerr << "Invalid rotated image dimensions " << w << " x " << h << endl;
      cerr << " expected " << (int) ceil(square.width * sqrt(2)) << " x " << (int) ceil(square.height * sqrt(2)) << endl;
      return false;
    }


    Image rotated(w, h);

    DPoint const a_p45 = convert_abs_coord(angle + rotation, ratio);
    APoint const p45 = convert_img_coord(rotated, a_p45);
    if (!fequal(0, p45.x, sigma))
    {
      cerr << __LINE__ << " > Rotation origin by 45 degrees:" << endl;
      cerr << " invalid x value: " << p45.x << " != " << 0 << endl;
      cerr << " absolute point: " << a_p45 << endl;
      cerr << " relative point: " << p45 << endl;
      return false;
    }
    if (!fequal(p45.y, h / 2.0f - 1, sigma))
    {
      cerr << __LINE__ << " > Rotation origin by 45 degrees:" << endl;
      cerr << "Invalid y value: " << p45.y << " != " << h / 2.0f - 1 << endl;
      cerr << " absolute point: " << a_p45 << endl;
      cerr << " relative point: " << p45 << endl;
      return false;
    }
  }

  return true;
}

void check_lines()
{
  Image const square(500, 500);
  draw_outline(square, 5, "square");
  draw_outline(square, 12, "square");
  draw_outline(square, 22, "square");
  draw_outline(square, 33, "square");
  draw_outline(square, 45, "square");
  draw_outline(square, 60, "square");
  draw_outline(square, 75, "square");
  draw_outline(square, 90, "square");

  Image const rect1(640, 480);
  draw_outline(rect1, 22, "rect1");
  draw_outline(rect1, 33, "rect1");
  draw_outline(rect1, 45, "rect1");
  draw_outline(rect1, 90, "rect1");
}



//
//
// Main
//

int main()
{
  bool perform_check = false;

  if (perform_check)
  {
    if (!check_points())
      return 1;

    if (!check_trigo())
      return 1;

    check_lines();
  }

  Image img("img/luigi.ppm");

  for (double rotation : {0, 15, 30})
  {
    auto const before = chrono::high_resolution_clock::now();

    Image rotated = rotate(img, rotation);

    auto const after = chrono::high_resolution_clock::now();
    auto const duration_ms = std::chrono::duration_cast<std::chrono::milliseconds>(after - before);
    cout << "rotate(): " << duration_ms.count() << " ms" << endl;

    stringstream filename;
    filename << "rotated_" << rotation << ".ppm";
    rotated.save(filename.str());
  }

  return 0;
}