--- /dev/null
+
+///////////////////////////////////////////////////////////////////////////
+// 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);
+}