commit 4bf69e820467582456060c25aec5e07a9eebdbb8
Author: Fabien Freling <public@ffreling.com>
Date:   Wed Jun 18 08:47:10 2014 +0200

    Initial commit

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..b87d621
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+.DS_Store
+*.dSYM
+*.swp
+
+*.ppm
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..0ecdb90
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,9 @@
+CXX     = clang++
+CXXFLAGS = -std=c++11 -W -Wall -Werror -O3
+BUILD_DIR=/tmp
+
+all: rotation.cpp
+	$(CXX) $(CXXFLAGS) $< -o $(BUILD_DIR)/rotation
+
+clean:
+	@rm -f *~ *.o .*.swp
diff --git a/TODO.md b/TODO.md
new file mode 100644
index 0000000..e69de29
diff --git a/img/luigi.jpg b/img/luigi.jpg
new file mode 100644
index 0000000..4bea30d
Binary files /dev/null and b/img/luigi.jpg differ
diff --git a/rotation.cpp b/rotation.cpp
new file mode 100644
index 0000000..0377584
--- /dev/null
+++ b/rotation.cpp
@@ -0,0 +1,486 @@
+#include <string>
+#include <fstream>
+#include <iostream>
+#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)
+  {
+    int const index = y * width + x;
+    r_chan[index] = r;
+    g_chan[index] = g;
+    b_chan[index] = b;
+  }
+
+  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;
+    }
+};
+
+
+//
+//
+// Drawing
+//
+
+void draw_line(Image& img, unsigned int x1, unsigned int y1, unsigned int x2, unsigned int y2)
+{
+
+  assert(x1 < x2);
+  double const slope = (double) (y2 - y1) / (double) (x2 - x1);
+  for (unsigned int i = x1; i <= x2; ++i)
+  {
+    unsigned int y = slope * (i - x1) + y1;
+    img.set_pixel(i, y, 255, 0, 0); // set line to red
+  }
+}
+
+void draw_line(Image& img, Point const& p1, Point const& p2)
+{
+  draw_line(img, p1.x, p1.y, p2.x, p2.y);
+}
+
+Point rotate(Image const& img, Point const& src, double radian, double const ratio)
+{
+  cout << "rotate (" << src.x << ", " << src.y << ") x " << radian << endl;
+  double x = src.x - (img.width / 2.0f);
+  x *= ratio;
+  // double y = - src.y + (src.height / 2.0f);
+  // y *= ratio;
+
+  double const angle_value = acos(x) + radian;
+  double const cos_x = cos(angle_value);
+  double const sin_x = sin(angle_value);
+  unsigned int const new_x = ceil(cos_x / ratio);
+  unsigned int const new_y = ceil(sin_x / ratio);
+
+  cout << " = (" << new_x << ", " << new_y << ")" << endl;
+  return Point(new_x, new_y);
+}
+
+
+//
+//
+// Trigonometry
+//
+
+double convert_radian(Image const& img, Point const& p, double const ratio)
+{
+  cout << "X: " << p.x << " - " << img.width / 2.0f << endl;
+  cout << "Y: " << - (int) p.y << " + " << img.height / 2.0f << endl;
+  double const centered_x = p.x - (img.width / 2.0f);
+  double const centered_y = (- (int) p.y) + (img.height / 2.0f);
+  cout << "centered point (" << centered_x << ", " << centered_y << ")" << endl;
+  double const cos_value = centered_x * ratio;
+  double const sin_value = centered_y * ratio;
+  double angle = acos(cos_value);
+  if (sin_value < 0)
+    angle = 2 * M_PI - angle;
+
+  return angle;
+}
+
+APoint convert_abs_coord(double const angle, double const ratio)
+{
+  cout << "Angle: " << angle << endl;
+  cout << "cos: " << cos(angle) << endl;
+  cout << "sin: " << sin(angle) << endl;
+  APoint tmp((int) ceil(cos(angle) / ratio), (int) ceil(sin(angle) / ratio));
+  cout << "point: " << tmp << endl;
+  return APoint((int) ceil(cos(angle) / ratio), (int) ceil(sin(angle) / ratio));
+}
+
+Point convert_img_coord(Image const& img, APoint const& p)
+{
+  cout << "image: " << img.width << " x " << img.height << endl;
+  cout << "p: " << p << endl;
+  cout << "h / 2: " << img.height / 2.0f << endl;
+  return Point(p.x + img.width / 2.0f, - p.y + img.height / 2.0f);
+}
+
+double compute_ratio(Image const& img)
+{
+  cout << "Compute ratio" << endl;
+  unsigned int const nb_points = img.width * img.height;
+  cout << " " << nb_points << " points" << endl;
+  double const square_side = sqrt(nb_points) - 1;
+  cout << " square side: " << square_side << endl;
+  double const half_side = square_side / 2;
+  cout << " half side: " << half_side << endl;
+  unsigned int const trigo_length = (unsigned int) ceil(half_side * sqrt(2));
+  cout << " trigo length: " << trigo_length << endl;
+  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);
+  int min_w = 0;
+  int max_w = 0;
+  int min_h = 0;
+  int max_h = 0;
+
+  cout << "Image dimensions: " << src.width << " x " << src.height << endl;
+  Point p(0, 0);
+  double angle = convert_radian(src, p, ratio);
+  APoint 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);
+  cout << "Rotated " << p << " (" << angle << ") = " << tl << "(" << angle + rotation << ")" << endl << endl;
+
+  p = Point(src.width - 1, 0);
+  angle = convert_radian(src, p, ratio);
+  APoint 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);
+  cout << "Rotated " << p << " (" << angle << ") = " << tr << "(" << angle + rotation << ")" << endl << endl;
+
+  p = Point(0, src.height - 1);
+  angle = convert_radian(src, p, ratio);
+  APoint 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);
+  cout << "Rotated " << p << " (" << angle << ") = " << bl << "(" << angle + rotation << ")" << endl << endl;
+
+  p = Point(src.width - 1, src.height - 1);
+  angle = convert_radian(src, p, ratio);
+  APoint 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);
+  cout << "Rotated " << p << " (" << angle << ") = " << br << "(" << angle + rotation << ")" << endl << endl;
+
+  width = max_w - min_w + 1;
+  height = max_h - min_h + 1;
+}
+
+bool check_trigo()
+{
+  Image square(500, 500);
+  double const ratio = compute_ratio(square);
+
+  // 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 (abs(angle - (3 * M_PI / 4)) > sigma)
+  {
+    cout << "Invalid angle value: " << angle << " != " << 3 * M_PI / 4 << endl;
+    return false;
+  }
+
+  // Check that we can reverse the origin point.
+  APoint const abs_reverse_point = convert_abs_coord(angle, ratio);
+  cout << "reversed abs origin: " << abs_reverse_point << endl;
+  Point 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 (!fequal(w, square.width * sqrt(2), sigma)
+//      || !fequal(h, square.height * sqrt(2), sigma))
+//  {
+//    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);
+
+  APoint const a_p45 = convert_abs_coord(angle + rotation, ratio);
+  Point const p45 = convert_img_coord(rotated, a_p45);
+  if (abs((float) (-1) - p45.x) > sigma)
+  {
+    cerr << "Rotation origin by 45 degrees:" << endl;
+    cerr << "Invalid x value: " << p45.x << " != " << -1 << endl;
+    cerr << "Absolute point: (" << a_p45.x << ", " << a_p45.y << ")" << endl;
+    return false;
+  }
+  if (abs(0.0 - p45.y) > sigma)
+  {
+    cerr << "Rotation origin by 45 degrees:" << endl;
+    cerr << "Invalid y value: " << p45.y << " != " << 0 << endl;
+    return false;
+  }
+
+  return true;
+}
+
+
+//
+//
+// Rotation
+//
+
+Image rotate(Image const& src, double angle)
+{
+  unsigned int const nb_points = src.width * src.height;
+  double const square_side = sqrt(nb_points);
+  double const half_side = square_side / 2;
+  unsigned int const trigo_length = (unsigned int) ceil(half_side * sqrt(2));
+  double const ratio = 1.0f / trigo_length;
+
+  // top left
+  double const cos_value = - (src.width / 2.0f) * ratio;
+//  double const sin_value = (src.height / 2.0f) * ratio;
+  double const angle_value = acos(cos_value);
+  cout << "top left angle: " << angle_value << endl;
+
+  double const radian = (angle / 360.0f) * (2 * M_PI);
+
+
+  // FIXME
+  Image rot(trigo_length, trigo_length);
+  Point tl = rotate(src, Point(0, 0), radian, ratio);
+  Point tr = rotate(src, Point(src.width - 1, 0), radian, ratio);
+  draw_line(rot, tl, tr);
+  // draw_line(rot, 0, 0, rot.width - 1, 0);
+  // draw_line(rot, 0, rot.height - 1, rot.width - 1, rot.height - 1);
+  // draw_line(rot, 0, 0, rot.width - 1, rot.height - 1);
+  // rot.set_pixel(rot.width - 1, rot.height - 1, 0, 255, 0);
+//  draw_line(rot, 0, 0, 0, rot.height - 1);
+//  draw_line(rot, rot.width - 1, 0, rot.width - 1, rot.height - 1);
+  return rot;
+}
+
+
+//
+//
+// Main
+//
+
+int main()
+{
+  if (!check_trigo())
+    return 1;
+
+  return 0;
+
+  Image img("img/luigi.ppm");
+
+  Image rotated = rotate(img, 30);
+  rotated.save("rotated.ppm");
+
+  return 0;
+}