automatic commit
[folded-ctf.git] / rich_image.cc
diff --git a/rich_image.cc b/rich_image.cc
new file mode 100644 (file)
index 0000000..f3b4b54
--- /dev/null
@@ -0,0 +1,453 @@
+
+///////////////////////////////////////////////////////////////////////////
+// This program is free software: you can redistribute it and/or modify  //
+// it under the terms of the version 3 of the GNU General Public License //
+// as published by the Free Software Foundation.                         //
+//                                                                       //
+// This program is distributed in the hope that it will be useful, but   //
+// WITHOUT ANY WARRANTY; without even the implied warranty of            //
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU      //
+// General Public License for more details.                              //
+//                                                                       //
+// You should have received a copy of the GNU General Public License     //
+// along with this program. If not, see <http://www.gnu.org/licenses/>.  //
+//                                                                       //
+// Written by Francois Fleuret, (C) IDIAP                                //
+// Contact <francois.fleuret@idiap.ch> for comments & bug reports        //
+///////////////////////////////////////////////////////////////////////////
+
+#include <string.h>
+
+#include "rich_image.h"
+
+#define PIXEL_DELTA(a, b) ((a) >= (b) ? (a) - (b) : (b) - (a))
+
+static inline bool edge(                   unsigned char v0, unsigned char v1,
+                         unsigned char v2, unsigned char v3, unsigned char v4, unsigned char v5,
+                                           unsigned char v6, unsigned char v7) {
+  unsigned char g = PIXEL_DELTA(v3, v4);
+
+  return
+    g > 8 &&
+    g > PIXEL_DELTA(v0, v3) &&
+    g > PIXEL_DELTA(v1, v4) &&
+    g > PIXEL_DELTA(v2, v3) &&
+    g > PIXEL_DELTA(v4, v5) &&
+    g > PIXEL_DELTA(v3, v6) &&
+    g > PIXEL_DELTA(v4, v7);
+}
+
+void RichImage::free() {
+  delete[] _edge_map;
+  _edge_map = 0;
+  delete[] _line_edge_map;
+  _line_edge_map = 0;
+  delete[] _tag_edge_map;
+  _tag_edge_map = 0;
+  delete[] _scale_edge_map;
+  _scale_edge_map = 0;
+  delete[] _width_at_scale;
+  _width_at_scale = 0;
+  delete[] _height_at_scale;
+  _height_at_scale = 0;
+}
+
+void RichImage::compute_one_scale_edge_maps(int width, int height,
+                                            unsigned int ***scale_edge_map,
+                                            unsigned char *pixel_map,
+                                            unsigned int *sum_pixel_map,
+                                            unsigned int *sum_sq_pixel_map) {
+
+  unsigned char d00, d01, d02, d03, d04, d05, d06, d07;
+  unsigned char d08, d09, d10, d11, d12, d13, d14, d15;
+  unsigned char *local_pixel_map;
+  unsigned int *local_sum_pixel_map, *local_sum_sq_pixel_map;
+
+  // Compute the integral images
+
+  {
+    unsigned int *sp = sum_pixel_map, *ssp = sum_sq_pixel_map;
+    unsigned char *p = pixel_map;
+    for(int x = 0; x < width; x++) {
+      *sp++ = 0;
+      *ssp++ = 0;
+      p++;
+    }
+
+    for(int y = 1; y < height; y++) {
+      *sp++ = 0;
+      *ssp++ = 0;
+      p++;
+      for(int x = 1; x < width; x++) {
+        *sp = *(sp - width) + *(sp - 1) - *(sp - width - 1) + ((unsigned int) (*p));
+        *ssp = *(ssp - width) + *(ssp - 1) - *(ssp - width - 1) + sq((unsigned int) (*p));
+        sp++;
+        ssp++;
+        p++;
+      }
+    }
+  }
+
+  const unsigned int var_square_size = 16;
+
+  int k00 = - 2 + width * (- 2);
+  int k01 = - 1 + width * (- 2);
+  int k02 = + 0 + width * (- 2);
+  int k03 = + 1 + width * (- 2);
+  int k04 = - 2 + width * (- 1);
+  int k05 = - 1 + width * (- 1);
+  int k06 = + 0 + width * (- 1);
+  int k07 = + 1 + width * (- 1);
+  int k08 = - 2 + width * (+ 0);
+  int k09 = - 1 + width * (+ 0);
+  int k10 = + 0 + width * (+ 0);
+  int k11 = + 1 + width * (+ 0);
+  int k12 = - 2 + width * (+ 1);
+  int k13 = - 1 + width * (+ 1);
+  int k14 = + 0 + width * (+ 1);
+  int k15 = + 1 + width * (+ 1);
+
+  unsigned int *edge_map0 = scale_edge_map[0][0];
+  unsigned int *edge_map1 = scale_edge_map[1][0];
+  unsigned int *edge_map2 = scale_edge_map[2][0];
+  unsigned int *edge_map3 = scale_edge_map[3][0];
+  unsigned int *edge_map4 = scale_edge_map[4][0];
+  unsigned int *edge_map5 = scale_edge_map[5][0];
+  unsigned int *edge_map6 = scale_edge_map[6][0];
+  unsigned int *edge_map7 = scale_edge_map[7][0];
+  unsigned int *variance_map = scale_edge_map[8][0];
+
+  int a = -1 - width, b = - width, c = -1, d = 0;
+
+  local_pixel_map = pixel_map;
+  local_sum_pixel_map = sum_pixel_map;
+  local_sum_sq_pixel_map = sum_sq_pixel_map;
+
+  const int gray_bin_width = 256 / nb_gray_tags;
+
+  for(int y = 0; y < height; y++) {
+
+    for(int x = 0; x < width; x++) {
+
+      if(x == 0 || y == 0) {
+        for(int e = 0; e < _nb_tags; e++)
+          scale_edge_map[e][0][d] = 0;
+      } else {
+        for(int e = 0; e < _nb_tags; e++)
+          scale_edge_map[e][0][d] =
+            scale_edge_map[e][0][b] +
+            scale_edge_map[e][0][c] -
+            scale_edge_map[e][0][a];
+      }
+
+      scale_edge_map[first_gray_tag +
+                     (local_pixel_map[0] / gray_bin_width)][0][d]++;
+
+      if(x - int(var_square_size/2) >= 0 &&
+         x + int(var_square_size/2) < width &&
+         y - int(var_square_size/2) >= 0 &&
+         y + int(var_square_size/2) < height) {
+
+        unsigned int s =
+          + local_sum_pixel_map[ - var_square_size/2 + width * ( - var_square_size / 2)]
+          + local_sum_pixel_map[ + var_square_size/2 + width * ( + var_square_size / 2)]
+          - local_sum_pixel_map[ - var_square_size/2 + width * ( + var_square_size / 2)]
+          - local_sum_pixel_map[ + var_square_size/2 + width * ( - var_square_size / 2)];
+
+        unsigned int s_sq =
+          + local_sum_sq_pixel_map[ - var_square_size/2 + width * ( - var_square_size / 2)]
+          + local_sum_sq_pixel_map[ + var_square_size/2 + width * ( + var_square_size / 2)]
+          - local_sum_sq_pixel_map[ - var_square_size/2 + width * ( + var_square_size / 2)]
+          - local_sum_sq_pixel_map[ + var_square_size/2 + width * ( - var_square_size / 2)];
+
+        if(sq(var_square_size) * s_sq - sq(s) >=
+           100 * sq(var_square_size) * (sq(var_square_size) - 1)) {
+
+          d00 = local_pixel_map[k00];
+          d01 = local_pixel_map[k01];
+          d02 = local_pixel_map[k02];
+          d03 = local_pixel_map[k03];
+
+          d04 = local_pixel_map[k04];
+          d05 = local_pixel_map[k05];
+          d06 = local_pixel_map[k06];
+          d07 = local_pixel_map[k07];
+
+          d08 = local_pixel_map[k08];
+          d09 = local_pixel_map[k09];
+          d10 = local_pixel_map[k10];
+          d11 = local_pixel_map[k11];
+
+          d12 = local_pixel_map[k12];
+          d13 = local_pixel_map[k13];
+          d14 = local_pixel_map[k14];
+          d15 = local_pixel_map[k15];
+
+          /*
+
+              #0         #1         #2         #3
+
+            XXXXXX     .XXXXX     ...XXX     .....X
+            XXXXXX     ..XXXX     ...XXX     ....XX
+            ......     ...XXX     ...XXX     ...XXX
+            ......     ....XX     ...XXX     ..XXXX
+
+              #4         #5         #6         #7
+
+            ......     X.....     XXX...     XXXXX.
+            ......     XX....     XXX...     XXXX..
+            XXXXXX     XXX...     XXX...     XXX...
+            XXXXXX     XXXX..     XXX...     XX....
+
+          */
+
+          if(edge(d04, d08, d01, d05, d09, d13, d06, d10)) {
+            if(d05 < d09) edge_map0[d]++;
+            else          edge_map4[d]++;
+          }
+
+          if(edge(d02, d07, d00, d05, d10, d15, d08, d13)) {
+            if(d05 < d10) edge_map7[d]++;
+            else          edge_map3[d]++;
+          }
+
+          if(edge(d01, d02, d04, d05, d06, d07, d09, d10)) {
+            if(d05 < d06) edge_map6[d]++;
+            else          edge_map2[d]++;
+          }
+
+          if(edge(d01, d04, d03, d06, d09, d12, d11, d14)) {
+            if(d06 < d09) edge_map1[d]++;
+            else          edge_map5[d]++;
+          }
+
+          variance_map[d]++;
+        }
+      }
+
+      a++; b++; c++; d++;
+      local_pixel_map++;
+      local_sum_pixel_map++;
+      local_sum_sq_pixel_map++;
+    }
+  }
+}
+
+void RichImage::compute_rich_structure() {
+  free();
+
+  unsigned char *pixel_maps;
+  unsigned char **line_pixel_maps;
+  unsigned char ***scale_pixel_maps;
+
+  _nb_scales =
+    int(global.nb_scales_per_power_of_two
+        * log(scalar_t(min(_width, _height))/scalar_t(_image_size_min)) / log(2)) + 1;
+
+  _width_at_scale = new int[_nb_scales];
+  _height_at_scale = new int[_nb_scales];
+  int total_surface = 0, total_lines = 0;
+
+  // Compute the sizes of the image at various scales
+
+  scalar_t rho = exp(log(0.5) / scalar_t(global.nb_scales_per_power_of_two));
+  scalar_t factor = 1.0;
+
+  for(int s = 0; s < _nb_scales; s++) {
+    _width_at_scale[s] = int(scalar_t(_width) * factor);
+    _height_at_scale[s] = int(scalar_t(_height) * factor);
+    total_surface += _width_at_scale[s] * _height_at_scale[s];
+    total_lines += _height_at_scale[s];
+    factor *= rho;
+  }
+
+  // Allocate the memory for the various scales and set up the pointer
+  // arrays to speed-up accesses
+
+  pixel_maps = new unsigned char[total_surface];
+  memset(pixel_maps, 0, sizeof(unsigned char) * total_surface);
+  line_pixel_maps = new unsigned char *[total_lines];
+  scale_pixel_maps = new unsigned char **[_nb_scales];
+
+  int sum_surfaces, sum_heights;
+
+  sum_surfaces = 0;
+  sum_heights = 0;
+
+  for(int s = 0; s < _nb_scales; s++) {
+    scale_pixel_maps[s] = line_pixel_maps + sum_heights;
+    for(int y = 0; y < _height_at_scale[s]; y++) {
+      scale_pixel_maps[s][y] = pixel_maps + sum_surfaces;
+      sum_surfaces += _width_at_scale[s];
+    }
+    sum_heights += _height_at_scale[s];
+  }
+
+  _edge_map = new unsigned int[total_surface * _nb_tags];
+  memset(_edge_map, 0, sizeof(unsigned int) * total_surface * _nb_tags);
+  _line_edge_map = new unsigned int *[total_lines * _nb_tags];
+  _tag_edge_map = new unsigned int **[_nb_tags * _nb_scales];
+  _scale_edge_map = new unsigned int ***[_nb_scales];
+
+  sum_surfaces = 0;
+  sum_heights = 0;
+
+  for(int s = 0; s < _nb_scales; s++) {
+    _scale_edge_map[s] = _tag_edge_map + s * _nb_tags;
+    for(int t = 0; t < _nb_tags; t++) {
+      _scale_edge_map[s][t] = _line_edge_map + sum_heights;
+      for(int y = 0; y < _height_at_scale[s]; y++) {
+        _scale_edge_map[s][t][y] = _edge_map + sum_surfaces;
+        sum_surfaces += _width_at_scale[s];
+      }
+      sum_heights += _height_at_scale[s];
+    }
+  }
+
+  // Compute the scaled images per se
+
+  for(int s = 0; s < _nb_scales; s++) {
+
+    if(s < global.nb_scales_per_power_of_two) {
+
+      if(s == 0) {
+
+        // The first image is not rescaled
+
+        for(int y = 0; y < _height_at_scale[s]; y++)
+          for(int x = 0; x < _width_at_scale[s]; x++)
+            scale_pixel_maps[s][y][x] = _content[x + _width * y];
+
+      } else {
+
+        // The nb_scales_per_power_of_two - 1 first images are
+        // generated with a 'complex' scaling, and then we just
+        // downscale by a factor of two
+
+        const int ld = 3;
+
+        int delta_column[_width_at_scale[s] << ld];
+        for(int xx = 0; xx < _width_at_scale[s] << ld; xx++)
+          delta_column[xx] = (xx * _width_at_scale[0]) / (_width_at_scale[s] << ld);
+
+        unsigned char *line[_height_at_scale[s] << ld];
+        for(int yy = 0; yy < _height_at_scale[s] << ld; yy++)
+          line[yy] = scale_pixel_maps[0][(yy * _height_at_scale[0]) / (_height_at_scale[s] << ld)];
+
+        for(int y = 0; y < _height_at_scale[s]; y++) {
+          unsigned char *dest_line = scale_pixel_maps[s][y];
+          for(int x = 0; x < _width_at_scale[s]; x++) {
+            int u = 0;
+            for(int yy = (y << ld); yy < ((y+1) << ld); yy++)
+              for(int xx = (x << ld); xx < ((x+1) << ld); xx++)
+                u += line[yy][delta_column[xx]];
+            dest_line[x] = (u >> (2 * ld));
+          }
+        }
+      }
+    } else {
+
+      // Other scales are computed by halfing one of the already
+      // computed scales
+
+      int r = s - global.nb_scales_per_power_of_two;
+      for(int y = 0; y < _height_at_scale[s]; y++)
+        for(int x = 0; x < _width_at_scale[s]; x++)
+          scale_pixel_maps[s][y][x] =
+            (int(scale_pixel_maps[r][y*2 + 0][x*2 + 0]) +
+             int(scale_pixel_maps[r][y*2 + 0][x*2 + 1]) +
+             int(scale_pixel_maps[r][y*2 + 1][x*2 + 0]) +
+             int(scale_pixel_maps[r][y*2 + 1][x*2 + 1])) >> 2;
+    }
+
+  }
+
+  // Allocate the memory for the edge maps at various scales and set
+  // up the pointer arrays to speed-up accesses
+
+  unsigned int *sum_pixel_map = new unsigned int[_width * _height];
+  unsigned int *sum_sq_pixel_map = new unsigned int[_width * _height];
+
+  for(int s = 0; s < _nb_scales; s++)
+    compute_one_scale_edge_maps(_width_at_scale[s],
+                                _height_at_scale[s],
+                                _scale_edge_map[s],
+                                scale_pixel_maps[s][0],
+                                sum_pixel_map,
+                                sum_sq_pixel_map);
+
+  delete[] sum_pixel_map;
+  delete[] sum_sq_pixel_map;
+
+  delete[] pixel_maps;
+  delete[] line_pixel_maps;
+  delete[] scale_pixel_maps;
+}
+
+void RichImage::crop(int xmin, int ymin, int width, int height) {
+  free();
+  Image::crop(xmin, ymin, width, height);
+}
+
+RichImage::RichImage() : Image() {
+  _width_at_scale = 0;
+  _height_at_scale = 0;
+  _edge_map = 0;
+  _line_edge_map = 0;
+  _tag_edge_map = 0;
+  _scale_edge_map = 0;
+}
+
+RichImage::RichImage(int width, int height) : Image(width, height) {
+  _width_at_scale = 0;
+  _height_at_scale = 0;
+  _edge_map = 0;
+  _line_edge_map = 0;
+  _tag_edge_map = 0;
+  _scale_edge_map = 0;
+}
+
+RichImage::~RichImage() {
+  delete[] _edge_map;
+  delete[] _line_edge_map;
+  delete[] _tag_edge_map;
+  delete[] _scale_edge_map;
+  delete[] _width_at_scale;
+  delete[] _height_at_scale;
+}
+
+void RichImage::write(ostream *out) {
+  Image::write(out);
+}
+
+void RichImage::read(istream *in) {
+  Image::read(in);
+}
+
+void RichImage::write_tag_png(const char *filename) {
+  const int nb_cols= 4;
+  const int nb_rows = ((nb_tags() + nb_cols - 1) / nb_cols);
+
+  int height_total = 0;
+  for(int s = 0; s < _nb_scales; s++) height_total += _height_at_scale[s];
+
+  RGBImage image(_width * nb_cols, height_total * nb_rows);
+
+  for(int y = 0; y < image.height(); y++) for(int x = 0; x < image.width(); x++)
+                                            image.set_pixel(x, y, 0, 255, 0);
+
+  int sum_height = 0;
+
+  for(int s = 0; s < _nb_scales; s++) {
+    for(int t = 0; t < nb_tags(); t++) {
+      int x0 = (t%nb_cols) * _width, y0 = sum_height + (t/nb_cols) * _height_at_scale[s];
+      for(int y = 0; y < _height_at_scale[s]; y++) for(int x = 0; x < _width_at_scale[s]; x++) {
+          int c = 255 - min(255, max(0, int(255 * nb_tags_in_window(s, t, x, y, x+1, y+1))));
+          image.set_pixel(x0 + x, y0 + y, c, c, c);
+        }
+    }
+    sum_height += nb_rows * _height_at_scale[s];
+  }
+
+  image.write_png(filename);
+}