edges.c

   1 /* Copyright (C) <2010> Douglas Bagnall <douglas@halo.gen.nz>
   2  *
   3  * This library is free software; you can redistribute it and/or
   4  * modify it under the terms of the GNU Library General Public
   5  * License as published by the Free Software Foundation; either
   6  * version 2 of the License, or (at your option) any later version.
   7  *
   8  * This library is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11  * Library General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Library General Public
  14  * License along with this library; if not, write to the
  15  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  16  * Boston, MA 02111-1307, USA.
  17  */
  18
  19
  20 #include "sparrow.h"
  21 #include "gstsparrow.h"
  22 #include "edges.h"
  23
  24 #include <string.h>
  25 #include <math.h>
  26
  27 #include "cv.h"
  28
  29
  30 /*
  31 */
  32
  33 #define SIG_WEIGHT 5
  34
  35 /*3 pixels manhatten distance makes you an outlier */
  36 #define OUTLIER_THRESHOLD 3 << (SPARROW_FIXED_POINT)
  37 #define OUTLIER_PENALTY 8
  38
  39 /*
  40 */
  41
  42
  43 #define SPARROW_MAP_LUT_SHIFT 1
  44 #define SPARROW_FP_2_LUT (SPARROW_FIXED_POINT - SPARROW_MAP_LUT_SHIFT)
  45
  46 #define OFFSET(x, y, w)((((y) * (w)) >> SPARROW_FIXED_POINT) + ((x) >> SPARROW_FIXED_POINT))
  47
  48 static void corners_to_lut(GstSparrow *sparrow, sparrow_find_lines_t *fl){
  49   sparrow_map_t *map = &sparrow->map; /*rows in sparrow->out */
  50   guint8 *mask = sparrow->screenmask; /*mask in sparrow->in */
  51   sparrow_corner_t *mesh = fl->mesh;   /*maps regular points in ->out to points in ->in */
  52
  53   size_t point_memsize = (sizeof(sparrow_map_point_t) * sparrow->out.pixcount / LINE_PERIOD) + 1;
  54   size_t row_memsize = sizeof(sparrow_map_row_t) * sparrow->out.height + 1;
  55   map->point_mem = malloc_aligned_or_die(point_memsize);
  56   map->rows = zalloc_aligned_or_die(row_memsize);
  57
  58   int mesh_w = fl->n_vlines;
  59   int mesh_h = fl->n_hlines;
  60   int in_w = sparrow->in.width;
  61   int mcy, mmy, mcx; /*Mesh Corner|Modulus X|Y*/
  62
  63   int x;
  64
  65   sparrow_map_row_t *row = map->rows;
  66   sparrow_map_point_t *p = map->point_mem;
  67   sparrow_corner_t *mesh_row = mesh;
  68   for(mcy = 0; mcy < mesh_h; mcy++){
  69     for (mmy = 0; mmy < LINE_PERIOD; mmy++){
  70       sparrow_corner_t *mesh_square = mesh_row;
  71       row->points = p;
  72       row->start = 0;
  73       row->end = 0;
  74       for(mcx = 0; mcx < mesh_w; mcx++){
  75         if (mesh_square->used){
  76           int iy = mesh_square->in_y + mmy * mesh_square->dyv;
  77           int ix = mesh_square->in_x + mmy * mesh_square->dxv;
  78           int ii = OFFSET(ix, iy, in_w);
  79           int ii_end = OFFSET(ix + (LINE_PERIOD - 1) * mesh_square->dxh,
  80               iy + (LINE_PERIOD - 1) * mesh_square->dyh, in_w);
  81           int start_on = mask[ii];
  82           int end_on = mask[ii_end];
  83           if(start_on && end_on){
  84             /*add the point, maybe switch on */
  85             if (row->start == row->end){/* if both are 0 */
  86               row->start = mcx * LINE_PERIOD;
  87             }
  88             p->x = ix;
  89             p->y = iy;
  90             p->dx = mesh_square->dxh;
  91             p->dy = mesh_square->dyh;
  92             p++;
  93           }
  94           else if (start_on){
  95             /*add the point, switch off somewhere in the middle*/
  96             for (x = 1; x < LINE_PERIOD; x++){
  97               iy += mesh_square->dyh;
  98               ix += mesh_square->dxh;
  99               ii = OFFSET(ix, iy, in_w);
 100               if (mask[ii]){
 101                 /*point is not in the same column with the others,
 102                   but sparrow knows this because the row->start says so */
 103                 row->start = mcx + x;
 104                 p->x = ix;
 105                 p->y = iy;
 106                 p->dx = mesh_square->dxh;
 107                 p->dy = mesh_square->dyh;
 108                 p++;
 109                 break;
 110               }
 111             }
 112           }
 113           else if (end_on){
 114             /* add some, switch off */
 115             for (x = 1; x < LINE_PERIOD; x++){
 116               iy += mesh_square->dyh;
 117               ix += mesh_square->dxh;
 118               ii = OFFSET(ix, iy, in_w);
 119               if (! mask[ii]){
 120                 row->end = mcx + x;
 121                 break;
 122               }
 123             }
 124           }
 125           else {
 126             /*3 cases:
 127               start > end: this is first off pixel.
 128               start == end: row hasn't started (both 0)
 129               start < end: both are set -- row is done
 130             */
 131             if (row->start > row->end){
 132               row->end = mcx * LINE_PERIOD;
 133             }
 134             else if (row->start < row->end){
 135               break;
 136             }
 137           }
 138         }
 139         mesh_square++;
 140       }
 141       row++;
 142     }
 143     mesh_row += mesh_w;
 144   }
 145 }
 146
 147 UNUSED static void
 148 corners_to_full_lut(GstSparrow *sparrow, sparrow_find_lines_t *fl){
 149   sparrow_corner_t *mesh = fl->mesh;   /*maps regular points in ->out to points in ->in */
 150   size_t lutsize = sizeof(sparrow_map_lut_t) * sparrow->out.pixcount;
 151   sparrow_map_lut_t *map_lut = malloc_aligned_or_die(lutsize);
 152
 153   int mesh_w = fl->n_vlines;
 154   int mesh_h = fl->n_hlines;
 155   int mcy, mmy, mcx, mmx; /*Mesh Corner|Modulus X|Y*/
 156   int i = 0;
 157   sparrow_corner_t *mesh_row = mesh;
 158   for(mcy = 0; mcy < mesh_h; mcy++){
 159     for (mmy = 0; mmy < LINE_PERIOD; mmy++){
 160       sparrow_corner_t *mesh_square = mesh_row;
 161       for(mcx = 0; mcx < mesh_w; mcx++){
 162         int iy = mesh_square->in_y + mmy * mesh_square->dyv;
 163         int ix = mesh_square->in_x + mmy * mesh_square->dxv;
 164         for (mmx = 0; mmx < LINE_PERIOD; mmx++, i++){
 165           map_lut[i].x = ix >> SPARROW_FP_2_LUT;
 166           map_lut[i].y = iy >> SPARROW_FP_2_LUT;
 167           ix += mesh_square->dxh;
 168           iy += mesh_square->dyh;
 169         }
 170         mesh_square++;
 171       }
 172     }
 173     mesh_row += mesh_w;
 174   }
 175   sparrow->map_lut = map_lut;
 176 }
 177
 178
 179
 180 /*create the mesh */
 181
 182 static void
 183 find_corners(GstSparrow *sparrow, guint8 *in, sparrow_find_lines_t *fl){
 184   int i;
 185   int width = fl->n_vlines;
 186   int height = fl->n_hlines;
 187
 188   sparrow_cluster_t *clusters = malloc_or_die(height * width * sizeof(sparrow_cluster_t));
 189   sparrow_corner_t *corners = zalloc_aligned_or_die(height * width * sizeof(sparrow_corner_t));
 190   gint x, y;
 191
 192   for (y = 0; y < sparrow->in.height; y++){
 193     for (x = 0; x < sparrow->in.width; x++){
 194       sparrow_intersect_t *p = &fl->map[y * sparrow->in.width + x];
 195       /*remembering that 0 is valid as a line no, but not as a signal */
 196       if (! p->signal[SPARROW_HORIZONTAL] ||
 197           ! p->signal[SPARROW_VERTICAL]){
 198         continue;
 199       }
 200       /*This one is lobbying for the position of a corner.*/
 201
 202       /*XXX what to do in the case that there is no intersection?  often cases
 203         this will happen in the dark bits and be OK. But if it happens in the
 204         light?*/
 205       /*linearise the xy coordinates*/
 206       int vline = p->lines[SPARROW_VERTICAL];
 207       int hline = p->lines[SPARROW_HORIZONTAL];
 208       sparrow_cluster_t *cluster = &clusters[vline * height + hline];
 209
 210       int n = cluster->n;
 211       if (n < 8){
 212         cluster->voters[n].x = x << SPARROW_FIXED_POINT;
 213         cluster->voters[n].y = y << SPARROW_FIXED_POINT;
 214         cluster->voters[n].signal = (SIG_WEIGHT + p->signal[SPARROW_HORIZONTAL]) *
 215           (SIG_WEIGHT + p->signal[SPARROW_VERTICAL]);
 216         cluster->n++;
 217         /*these next two could of course be computed from the offset */
 218       }
 219       else {
 220         GST_DEBUG("more than 8 pixels at cluster for corner %d, %d\n",
 221             vline, hline);
 222         /*if this message ever turns up, replace the weakest signals or add
 223           more slots.*/
 224       }
 225     }
 226   }
 227
 228   i = 0;
 229   for (y = 0; y < height; y++){
 230     for (x = 0; x < width; x++, i++){
 231       /* how to do this?
 232          1. centre of gravity (x,y, weighted average)
 233          2. discard outliers? look for connectedness? but if 2 are outliers?
 234       */
 235       sparrow_cluster_t *cluster = clusters + i;
 236       if (cluster->n == 0){
 237         continue;
 238       }
 239       int xsum, ysum;
 240       int xmean, ymean;
 241       int votes;
 242       while(1) {
 243         int j;
 244         xsum = 0;
 245         ysum = 0;
 246         votes = 0;
 247         for (j = 0; j < cluster->n; j++){
 248           votes += cluster->voters[j].signal;
 249           ysum += cluster->voters[j].y * cluster->voters[j].signal;
 250           xsum += cluster->voters[j].x * cluster->voters[j].signal;
 251         }
 252         if (votes == 0){
 253           /* don't diminish signal altogether. The previous iteration's means
 254              will be used. */
 255           break;
 256         }
 257         xmean = xsum / votes;
 258         ymean = ysum / votes;
 259         int worst = -1;
 260         int worstn;
 261         int devsum = 0;
 262         for (j = 0; j < cluster->n; j++){
 263           int xdiff = abs(cluster->voters[j].x - xmean);
 264           int ydiff = abs(cluster->voters[j].y - ymean);
 265           devsum += xdiff + ydiff;
 266           if (xdiff + ydiff > worst){
 267             worst = xdiff + ydiff;
 268             worstn = j;
 269           }
 270         }
 271         /*a bad outlier has significantly greater than average deviation
 272           (but how much is bad? median deviation would be more useful)*/
 273         if (worst > 3 * devsum / cluster->n){
 274           /* reduce the worst ones weight. it is a silly aberration. */
 275           cluster->voters[worstn].signal /= OUTLIER_PENALTY;
 276           GST_DEBUG("dropping outlier at %s,%s (mean %s,%s)\n",
 277               cluster->voters[worstn].x, cluster->voters[worstn].y, xmean, ymean);
 278           continue;
 279         }
 280         break;
 281       }
 282       corners[i].out_x = x * LINE_PERIOD;
 283       corners[i].out_y = y * LINE_PERIOD;
 284       corners[i].in_x = xmean;
 285       corners[i].in_y = ymean;
 286       corners[i].used = TRUE;
 287       double div = (double)(1 << SPARROW_FIXED_POINT); /*for printf only*/
 288       GST_DEBUG("found corner %d (%d,%d) at (%3f, %3f)\n",
 289           i, corners[i].out_x, corners[i].out_y,
 290           xmean / div, ymean / div);
 291     }
 292   }
 293   free(clusters);
 294
 295   /* calculate deltas toward adjacent corners */
 296   /* try to extrapolate left and up, if possible, so need to go backwards. */
 297   for (y = height - 2; y >= 0; y--){
 298     for (x = width - 2; x >= 0; x--){
 299       i = y * width + x;
 300       if (corners[i].used){
 301         corners[i].dxh = (corners[i + 1].in_x - corners[i].in_x) / LINE_PERIOD;
 302         corners[i].dyh = (corners[i + 1].in_y - corners[i].in_y) / LINE_PERIOD;
 303         corners[i].dxv = (corners[i + width].in_x - corners[i].in_x) / LINE_PERIOD;
 304         corners[i].dyv = (corners[i + width].in_y - corners[i].in_y) / LINE_PERIOD;
 305       }
 306       else {
 307           /*prefer copy from left unless it is itself reconstructed,
 308             (for no great reason)
 309             A mixed copy would be possible and better */
 310         if(corners[i + 1].used &&
 311             (corners[i + 1].used < corners[i + width].used)){
 312           corners[i].dxh = corners[i + 1].dxh;
 313           corners[i].dyh = corners[i + 1].dyh;
 314           corners[i].dxv = corners[i + 1].dxv;
 315           corners[i].dyv = corners[i + 1].dyv;
 316           corners[i].in_x = corners[i + 1].in_x - corners[i + 1].dxh * LINE_PERIOD;
 317           corners[i].in_y = corners[i + 1].in_y - corners[i + 1].dyh * LINE_PERIOD;
 318           corners[i].out_x = corners[i + 1].out_x - 1;
 319           corners[i].out_y = corners[i + 1].out_y;
 320           corners[i].used = corners[i + 1].used + 1;
 321         }
 322         else if(corners[i + width].used){
 323           corners[i].dxh = corners[i + width].dxh;
 324           corners[i].dyh = corners[i + width].dyh;
 325           corners[i].dxv = corners[i + width].dxv;
 326           corners[i].dyv = corners[i + width].dyv;
 327           corners[i].in_x = corners[i + width].in_x - corners[i + width].dxv * LINE_PERIOD;
 328           corners[i].in_y = corners[i + width].in_y - corners[i + width].dyv * LINE_PERIOD;
 329           corners[i].out_x = corners[i + width].out_x;
 330           corners[i].out_y = corners[i + width].out_y - 1;
 331           corners[i].used = corners[i + width].used + 1;
 332         }
 333       }
 334     }
 335   }
 336
 337   fl->mesh = corners;
 338 }
 339
 340
 341 /* With no line drawn (in our colour) look at the background noise.  Any real
 342    signal has to be stringer than this.
 343
 344    XXX looking for simple maximum -- maybe heap or histogram might be better,
 345    so as to be less susceptible to wierd outliers (e.g., bad pixels).  */
 346 static void
 347 look_for_threshold(GstSparrow *sparrow, guint8 *in, sparrow_find_lines_t *fl){
 348   int i;
 349   guint32 colour;
 350   guint32 signal;
 351   guint32 *in32 = (guint32 *)in;
 352   guint32 highest = 0;
 353   for (i = 0;  i < (int)sparrow->in.size; i++){
 354     colour = in32[i] & sparrow->colour;
 355     signal = ((colour >> fl->shift1) +
 356         (colour >> fl->shift2)) & 0x1ff;
 357     if (signal > highest){
 358       highest = signal;
 359     }
 360   }
 361   fl->threshold = highest + 10;
 362   GST_DEBUG("found maximum noise of %s, using threshold %s\n", highest, fl->threshold);
 363 }
 364
 365
 366 static void
 367 look_for_line(GstSparrow *sparrow, guint8 *in, sparrow_find_lines_t *fl,
 368     sparrow_line_t *line){
 369   guint i;
 370   guint32 colour;
 371   int signal;
 372   guint32 *in32 = (guint32 *)in;
 373   for (i = 0; i < sparrow->in.size; i++){
 374     colour = in32[i] & sparrow->colour;
 375     signal = ((colour >> fl->shift1) +
 376         (colour >> fl->shift2)) & 0x1ff;
 377     if (signal > fl->threshold){
 378       if (fl->map[i].lines[line->dir]){
 379         GST_DEBUG("HEY, expected point %d to be in line %d (dir %d)"
 380             "and thus empty, but it is also in line %d\n",
 381             i, line->index, line->dir, fl->map[i].lines[line->dir]);
 382       }
 383       fl->map[i].lines[line->dir] = line->index;
 384       fl->map[i].signal[line->dir] = signal;
 385     }
 386   }
 387 }
 388
 389
 390 /* draw the line (in sparrow->colour) */
 391 static inline void
 392 draw_line(GstSparrow * sparrow, sparrow_line_t *line, guint8 *out){
 393   guint32 *p = (guint32 *)out;
 394   int i;
 395   if (line->dir == SPARROW_HORIZONTAL){
 396     p += line->offset * sparrow->out.width;
 397     for (i = 0; i < sparrow->out.width; i++){
 398       p[i] = sparrow->colour;
 399     }
 400   }
 401   else {
 402     guint32 *p = (guint32 *)out;
 403     p += line->offset;
 404     for(i = 0; i < sparrow->out.height; i++){
 405       *p = sparrow->colour;
 406       p += sparrow->out.width;
 407     }
 408   }
 409 }
 410
 411 /* show each line for 2 frames, then wait sparrow->lag frames, leaving line on
 412    until last one.
 413  */
 414
 415 INVISIBLE sparrow_state
 416 mode_find_edges(GstSparrow *sparrow, guint8 *in, guint8 *out){
 417   sparrow_find_lines_t *fl = (sparrow_find_lines_t *)&sparrow->helper_struct;
 418   sparrow_line_t *line = fl->shuffled_lines[fl->current];
 419   sparrow->countdown--;
 420   memset(out, 0, sparrow->out.size);
 421   if (sparrow->countdown){
 422     /* show the line except on the first round, when we find a threshold*/
 423     if (fl->threshold){
 424       draw_line(sparrow, line, out);
 425     }
 426   }
 427   else{
 428       /*show nothing, look for result */
 429     if (fl->threshold){
 430       if (fl->current == fl->n_lines){
 431         goto done;
 432       }
 433       look_for_line(sparrow, in, fl, line);
 434       fl->current++;
 435     }
 436     else {
 437       look_for_threshold(sparrow, in, fl);
 438     }
 439     sparrow->countdown = sparrow->lag + 2;
 440   }
 441   return SPARROW_STATUS_QUO;
 442  done:
 443   /*match up lines and find corners */
 444   find_corners(sparrow, in, fl);
 445   corners_to_lut(sparrow, fl);
 446
 447   /*free stuff!, including fl, and reset pointer to NULL*/
 448
 449   return SPARROW_NEXT_STATE;
 450 }
 451
 452
 453 static void
 454 setup_colour_shifts(GstSparrow *sparrow, sparrow_find_lines_t *fl){
 455   switch (sparrow->colour){
 456   case SPARROW_WHITE:
 457   case SPARROW_GREEN:
 458     fl->shift1 = sparrow->in.gshift;
 459     fl->shift2 = sparrow->in.gshift;
 460     break;
 461   case SPARROW_MAGENTA:
 462     fl->shift1 = sparrow->in.rshift;
 463     fl->shift2 = sparrow->in.bshift;
 464     break;
 465   }
 466 }
 467
 468 INVISIBLE void
 469 init_find_edges(GstSparrow *sparrow){
 470   gint32 w = sparrow->out.width;
 471   gint32 h = sparrow->out.height;
 472   gint i;
 473   sparrow_find_lines_t *fl = zalloc_aligned_or_die(sizeof(sparrow_find_lines_t));
 474   sparrow->helper_struct = (void *)fl;
 475
 476   gint h_lines = (h + LINE_PERIOD - 1) / LINE_PERIOD;
 477   gint v_lines = (w + LINE_PERIOD - 1) / LINE_PERIOD;
 478   gint n_lines = (h_lines + v_lines);
 479   fl->n_hlines = h_lines;
 480   fl->n_vlines = v_lines;
 481   fl->n_lines = n_lines;
 482
 483   fl->h_lines = malloc_aligned_or_die(sizeof(sparrow_line_t) * n_lines);
 484   fl->shuffled_lines = malloc_aligned_or_die(sizeof(sparrow_line_t*) * n_lines);
 485   fl->map = zalloc_aligned_or_die(sizeof(sparrow_intersect_t) * sparrow->in.pixcount);
 486
 487   sparrow_line_t *line = fl->h_lines;
 488   sparrow_line_t **sline = fl->shuffled_lines;
 489   int offset = LINE_PERIOD / 2;
 490
 491   for (i = 0, offset = LINE_PERIOD / 2; offset < h;
 492        i++, offset += LINE_PERIOD){
 493     line->offset = offset;
 494     line->dir = SPARROW_HORIZONTAL;
 495     line->index = i;
 496     *sline = line;
 497     line++;
 498     sline++;
 499   }
 500
 501   /*now add the vertical lines */
 502   fl->v_lines = line;
 503   for (i = 0, offset = LINE_PERIOD / 2; offset < w;
 504        i++, offset += LINE_PERIOD){
 505     line->offset = offset;
 506     line->dir = SPARROW_VERTICAL;
 507     line->index = i;
 508     *sline = line;
 509     line++;
 510     sline++;
 511   }
 512
 513   GST_DEBUG("allocated %s lines, used %s\n", n_lines, line - fl->h_lines);
 514
 515   /*now shuffle (triangluar, to no particular advantage) */
 516   for (i = 0; i < n_lines - 1; i++){
 517     int j = RANDINT(sparrow, i + 1, n_lines);
 518     sparrow_line_t *tmp = fl->shuffled_lines[j];
 519     fl->shuffled_lines[j] = fl->shuffled_lines[i];
 520     fl->shuffled_lines[i] = tmp;
 521   }
 522
 523   setup_colour_shifts(sparrow, fl);
 524   sparrow->countdown = sparrow->lag + 2;
 525
 526   fl->mesh = malloc_aligned_or_die(sizeof(sparrow_corner_t) * h_lines * v_lines);
 527 }