diff --git a/image.cpp b/image.cpp
index 55e9bc6..5a74749 100644
--- a/image.cpp
+++ b/image.cpp
@@ -35,7 +35,6 @@ Image::Image(unsigned int w, unsigned int h, pnm::Format type)
   }
 
   buffer = new pvalue_t[width * height * pixel_size];
-  memset(buffer, 0, width * height * pixel_size * sizeof (pvalue_t));
 }
 
 Image::Image(string const& path)
diff --git a/rotation.cpp b/rotation.cpp
index f43b468..0ab9efd 100644
--- a/rotation.cpp
+++ b/rotation.cpp
@@ -123,6 +123,15 @@ void compute_output_size(Image const& src, double const rotation, unsigned int&
   height = (int) (max_h - min_h) + 1;
 }
 
+DPoint get_mapped_point(Image const& src, Point const& p, double const rotation)
+{
+  DPoint const d = convert_grid_coord(src, p);
+  double p_angle = 0;
+  double dist = 0;
+  convert_abs_to_polar_coord(d, p_angle, dist);
+  return convert_polar_to_grid_coord(p_angle + rotation, dist);
+}
+
 
 
 //
@@ -149,16 +158,79 @@ bool fequal(float a, float b, float sigma)
 
 //
 //
-// Image rotation
+// Padding
 //
 
-DPoint get_mapped_point(Image const& src, Point const& p, double const rotation)
+uint16_t* generate_padding_table(Image const& src,
+                                Image const& rotated,
+                                double rotation,
+                                int q_pos)
 {
-  DPoint const d = convert_grid_coord(src, p);
-  double p_angle = 0;
-  double dist = 0;
-  convert_abs_to_polar_coord(d, p_angle, dist);
-  return convert_polar_to_grid_coord(p_angle + rotation, dist);
+  uint16_t* padding_table = new uint16_t[rotated.height];
+
+  // We suppose the image is square.
+  double const sigma = 1.0e-5;
+  if (fequal(rotation, 0, sigma)
+      || fequal(rotation, M_PI / 2, sigma)
+      || fequal(rotation, M_PI, sigma)
+      || fequal(rotation, 3 * M_PI / 2, sigma))
+  {
+    memset(padding_table, 0, sizeof (uint16_t) * rotated.height);
+    return padding_table;
+  }
+
+  double padding_rotation = 0.0;
+  if (rotation < M_PI / 2)
+  {
+    padding_rotation = rotation;
+  } else if (rotation < M_PI)
+  {
+    padding_rotation = rotation - M_PI / 2;
+  } else if (rotation < 3 * M_PI / 2)
+  {
+    padding_rotation = rotation - M_PI;
+  }
+  else
+  {
+    padding_rotation = rotation - 3 * M_PI / 2;
+  }
+
+  // Quantized delta for padding
+  // TODO: we should be able to infer these values from qdx and qdy
+  DPoint const top_left_grid = get_mapped_point(src, Point(0, 0), padding_rotation);
+  Point const top_left = convert_img_coord(rotated, top_left_grid);
+  cout << "Padding top left: " << top_left << endl;
+
+  DPoint const top_right_grid = get_mapped_point(src, Point(src.width - 1, 0), padding_rotation);
+  DPoint const top_right = convert_img_coord_precision(rotated, top_right_grid);
+  Point const q_top_right(top_right.x * q_pos, top_right.y * q_pos);
+
+  DPoint const bottom_left_grid = get_mapped_point(src, Point(0, src.height), padding_rotation);
+  DPoint const bottom_left = convert_img_coord_precision(rotated, bottom_left_grid);
+  Point const q_bottom_left(bottom_left.x * q_pos, bottom_left.y * q_pos);
+
+  int i = 0;
+  int const q_top_width = q_top_right.x;
+  int const q_top_step = q_top_width / top_left.y;
+  int previous_padding = 0;
+  for (; i <= top_left.y; ++i)
+  {
+    padding_table[i] = (q_top_width - i * q_top_step) / q_pos;
+    previous_padding = padding_table[i];
+  }
+
+  int const remaining_height = rotated.height - 1 - i;
+  int const q_bottom_width = q_bottom_left.x;
+  int const q_bottom_step = q_bottom_width / remaining_height;
+  int remaining_index = 1;
+  for (; i < (int) rotated.height; ++i, ++remaining_index)
+  {
+    padding_table[i] = min((unsigned int) (remaining_index * q_bottom_step) / q_pos,
+                           (unsigned int) rotated.height - 1 - padding_table[rotated.height - 1 - i]);
+  }
+
+
+  return padding_table;
 }
 
 inline