The Ecere SDK - sdk/blob - ecere/src/gfx/drivers/edtaa3func.ec

   1 /*
   2  * edtaa3()
   3  *
   4  * Sweep-and-update Euclidean distance transform of an
   5  * image. Positive pixels are treated as object pixels,
   6  * zero or negative pixels are treated as background.
   7  * An attempt is made to treat antialiased edges correctly.
   8  * The input image must have pixels in the range [0,1],
   9  * and the antialiased image should be a box-filter
  10  * sampling of the ideal, crisp edge.
  11  * If the antialias region is more than 1 pixel wide,
  12  * the result from this transform will be inaccurate.
  13  *
  14  * By Stefan Gustavson (stefan.gustavson@gmail.com).
  15  *
  16  * Originally written in 1994, based on a verbal
  17  * description of the SSED8 algorithm published in the
  18  * PhD dissertation of Ingemar Ragnemalm. This is his
  19  * algorithm, I only implemented it in C.
  20  *
  21  * Updated in 2004 to treat border pixels correctly,
  22  * and cleaned up the code to improve readability.
  23  *
  24  * Updated in 2009 to handle anti-aliased edges.
  25  *
  26  * Updated in 2011 to avoid a corner case infinite loop.
  27  *
  28  * Updated 2012 to change license from LGPL to MIT.
  29  */
  30
  31 /*
  32  Copyright (C) 2009-2012 Stefan Gustavson (stefan.gustavson@gmail.com)
  33  The code in this file is distributed under the MIT license:
  34
  35  Permission is hereby granted, free of charge, to any person obtaining a copy
  36  of this software and associated documentation files (the "Software"), to deal
  37  in the Software without restriction, including without limitation the rights
  38  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  39  copies of the Software, and to permit persons to whom the Software is
  40  furnished to do so, subject to the following conditions:
  41
  42  The above copyright notice and this permission notice shall be included in
  43  all copies or substantial portions of the Software.
  44
  45  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  46  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  47  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  48  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  49  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  50  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  51  THE SOFTWARE.
  52  */
  53
  54 #ifdef BUILDING_ECERE_COM
  55 import "instance"
  56 #endif
  57
  58  // Note: Changed double to float from original version
  59
  60 /*
  61  * Compute the local gradient at edge pixels using convolution filters.
  62  * The gradient is computed only at edge pixels. At other places in the
  63  * image, it is never used, and it's mostly zero anyway.
  64  */
  65 void computegradient(float *img, int w, int h, float *gx, float *gy)
  66 {
  67     int i,j,k;
  68     float glength;
  69 #define SQRT2 1.4142136f
  70     for(i = 1; i < h-1; i++) { // Avoid edges where the kernels would spill over
  71         for(j = 1; j < w-1; j++) {
  72             k = i*w + j;
  73             if((img[k]>0.0) && (img[k]<1.0f)) { // Compute gradient for edge pixels only
  74                 gx[k] = -img[k-w-1] - SQRT2*img[k-1] - img[k+w-1] + img[k-w+1] + SQRT2*img[k+1] + img[k+w+1];
  75                 gy[k] = -img[k-w-1] - SQRT2*img[k-w] - img[k+w-1] + img[k-w+1] + SQRT2*img[k+w] + img[k+w+1];
  76                 glength = gx[k]*gx[k] + gy[k]*gy[k];
  77                 if(glength > 0.0f) { // Avoid division by zero
  78                     glength = (float)sqrt(glength);
  79                     gx[k]=gx[k]/glength;
  80                     gy[k]=gy[k]/glength;
  81                 }
  82             }
  83         }
  84     }
  85     // TODO: Compute reasonable values for gx, gy also around the image edges.
  86     // (These are zero now, which reduces the accuracy for a 1-pixel wide region
  87         // around the image edge.) 2x2 kernels would be suitable for this.
  88 }
  89
  90 /*
  91  * A somewhat tricky function to approximate the distance to an edge in a
  92  * certain pixel, with consideration to either the local gradient (gx,gy)
  93  * or the direction to the pixel (dx,dy) and the pixel greyscale value a.
  94  * The latter alternative, using (dx,dy), is the metric used by edtaa2().
  95  * Using a local estimate of the edge gradient (gx,gy) yields much better
  96  * accuracy at and near edges, and reduces the error even at distant pixels
  97  * provided that the gradient direction is accurately estimated.
  98  */
  99 float edgedf(float gx, float gy, float a)
 100 {
 101     float df, glength, temp, a1;
 102
 103     if ((gx == 0) || (gy == 0)) { // Either A) gu or gv are zero, or B) both
 104         df = 0.5f-a;  // Linear approximation is A) correct or B) a fair guess
 105     } else {
 106         glength = (float)sqrt(gx*gx + gy*gy);
 107         if(glength>0) {
 108             gx = gx/glength;
 109             gy = gy/glength;
 110         }
 111         /* Everything is symmetric wrt sign and transposition,
 112          * so move to first octant (gx>=0, gy>=0, gx>=gy) to
 113          * avoid handling all possible edge directions.
 114          */
 115         gx = (float)fabs(gx);
 116         gy = (float)fabs(gy);
 117         if(gx<gy) {
 118             temp = gx;
 119             gx = gy;
 120             gy = temp;
 121         }
 122         a1 = 0.5f*gy/gx;
 123         if (a < a1) { // 0 <= a < a1
 124             df = 0.5f*(gx + gy) - (float)sqrt(2.0*gx*gy*a);
 125         } else if (a < (1.0-a1)) { // a1 <= a <= 1-a1
 126             df = (0.5f-a)*gx;
 127         } else { // 1-a1 < a <= 1
 128             df = -0.5f*(gx + gy) + (float)sqrt(2.0*gx*gy*(1.0-a));
 129         }
 130     }
 131     return df;
 132 }
 133
 134 float distaa3(float *img, float *gximg, float *gyimg, int w, int c, int xc, int yc, int xi, int yi)
 135 {
 136   float di, df, dx, dy, gx, gy, a;
 137   int closest;
 138
 139   closest = c-xc-yc*w; // Index to the edge pixel pointed to from c
 140   a = img[closest];    // Grayscale value at the edge pixel
 141   gx = gximg[closest]; // X gradient component at the edge pixel
 142   gy = gyimg[closest]; // Y gradient component at the edge pixel
 143
 144   if(a > 1.0) a = 1.0;
 145   if(a < 0.0) a = 0.0; // Clip grayscale values outside the range [0,1]
 146   if(a == 0.0) return 1000000.0; // Not an object pixel, return "very far" ("don't know yet")
 147
 148   dx = (float)xi;
 149   dy = (float)yi;
 150   di = (float)sqrt(dx*dx + dy*dy); // Length of integer vector, like a traditional EDT
 151   if(di==0) { // Use local gradient only at edges
 152       // Estimate based on local gradient only
 153       df = edgedf(gx, gy, a);
 154   } else {
 155       // Estimate gradient based on direction to edge (accurate for large di)
 156       df = edgedf(dx, dy, a);
 157   }
 158   return di + df; // Same metric as edtaa2, except at edges (where di=0)
 159 }
 160
 161 // Shorthand macro: add ubiquitous parameters dist, gx, gy, img and w and call distaa3()
 162 #define DISTAA(c,xc,yc,xi,yi) (distaa3(img, gx, gy, w, c, xc, yc, xi, yi))
 163
 164 void edtaa3(float *img, float *gx, float *gy, int w, int h, short *distx, short *disty, float *dist)
 165 {
 166   int x, y, i, c;
 167   int offset_u, offset_ur, offset_r, offset_rd,
 168   offset_d, offset_dl, offset_l, offset_lu;
 169   float olddist, newdist;
 170   int cdistx, cdisty, newdistx, newdisty;
 171   int changed;
 172   float epsilon = 1e-3;
 173
 174   /* Initialize index offsets for the current image width */
 175   offset_u = -w;
 176   offset_ur = -w+1;
 177   offset_r = 1;
 178   offset_rd = w+1;
 179   offset_d = w;
 180   offset_dl = w-1;
 181   offset_l = -1;
 182   offset_lu = -w-1;
 183
 184   /* Initialize the distance images */
 185   for(i=0; i<w*h; i++) {
 186     distx[i] = 0; // At first, all pixels point to
 187     disty[i] = 0; // themselves as the closest known.
 188     if(img[i] <= 0.0)
 189       {
 190         dist[i]= 1000000.0; // Big value, means "not set yet"
 191       }
 192     else if (img[i]<1.0) {
 193       dist[i] = edgedf(gx[i], gy[i], img[i]); // Gradient-assisted estimate
 194     }
 195     else {
 196       dist[i]= 0.0; // Inside the object
 197     }
 198   }
 199
 200   /* Perform the transformation */
 201   do
 202     {
 203       changed = 0;
 204
 205       /* Scan rows, except first row */
 206       for(y=1; y<h; y++)
 207         {
 208
 209           /* move index to leftmost pixel of current row */
 210           i = y*w;
 211
 212           /* scan right, propagate distances from above & left */
 213
 214           /* Leftmost pixel is special, has no left neighbors */
 215           olddist = dist[i];
 216           if(olddist > 0) // If non-zero distance or not set yet
 217             {
 218               c = i + offset_u; // Index of candidate for testing
 219               cdistx = distx[c];
 220               cdisty = disty[c];
 221               newdistx = cdistx;
 222               newdisty = cdisty+1;
 223               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 224               if(newdist < olddist-epsilon)
 225                 {
 226                   distx[i]=(short)newdistx;
 227                   disty[i]=(short)newdisty;
 228                             dist[i]=(short)newdist;
 229                   olddist=newdist;
 230                   changed = 1;
 231                 }
 232
 233               c = i+offset_ur;
 234               cdistx = distx[c];
 235               cdisty = disty[c];
 236               newdistx = cdistx-1;
 237               newdisty = cdisty+1;
 238               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 239               if(newdist < olddist-epsilon)
 240                 {
 241                   distx[i]=(short)newdistx;
 242                   disty[i]=(short)newdisty;
 243                             dist[i]=(short)newdist;
 244                   changed = 1;
 245                 }
 246             }
 247           i++;
 248
 249           /* Middle pixels have all neighbors */
 250           for(x=1; x<w-1; x++, i++)
 251             {
 252               olddist = dist[i];
 253               if(olddist <= 0) continue; // No need to update further
 254
 255               c = i+offset_l;
 256               cdistx = distx[c];
 257               cdisty = disty[c];
 258               newdistx = cdistx+1;
 259               newdisty = cdisty;
 260               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 261               if(newdist < olddist-epsilon)
 262                 {
 263                   distx[i]=(short)newdistx;
 264                   disty[i]=(short)newdisty;
 265                             dist[i]=(short)newdist;
 266                   olddist=newdist;
 267                   changed = 1;
 268                 }
 269
 270               c = i+offset_lu;
 271               cdistx = distx[c];
 272               cdisty = disty[c];
 273               newdistx = cdistx+1;
 274               newdisty = cdisty+1;
 275               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 276               if(newdist < olddist-epsilon)
 277                 {
 278                   distx[i]=(short)newdistx;
 279                   disty[i]=(short)newdisty;
 280                             dist[i]=(short)newdist;
 281                   olddist=newdist;
 282                   changed = 1;
 283                 }
 284
 285               c = i+offset_u;
 286               cdistx = distx[c];
 287               cdisty = disty[c];
 288               newdistx = cdistx;
 289               newdisty = cdisty+1;
 290               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 291               if(newdist < olddist-epsilon)
 292                 {
 293                   distx[i]=(short)newdistx;
 294                   disty[i]=(short)newdisty;
 295                             dist[i]=(short)newdist;
 296                   olddist=newdist;
 297                   changed = 1;
 298                 }
 299
 300               c = i+offset_ur;
 301               cdistx = distx[c];
 302               cdisty = disty[c];
 303               newdistx = cdistx-1;
 304               newdisty = cdisty+1;
 305               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 306               if(newdist < olddist-epsilon)
 307                 {
 308                   distx[i]=(short)newdistx;
 309                   disty[i]=(short)newdisty;
 310                             dist[i]=(short)newdist;
 311                   changed = 1;
 312                 }
 313             }
 314
 315           /* Rightmost pixel of row is special, has no right neighbors */
 316           olddist = dist[i];
 317           if(olddist > 0) // If not already zero distance
 318             {
 319               c = i+offset_l;
 320               cdistx = distx[c];
 321               cdisty = disty[c];
 322               newdistx = cdistx+1;
 323               newdisty = cdisty;
 324               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 325               if(newdist < olddist-epsilon)
 326                 {
 327                   distx[i]=(short)newdistx;
 328                   disty[i]=(short)newdisty;
 329                             dist[i]=(short)newdist;
 330                   olddist=newdist;
 331                   changed = 1;
 332                 }
 333
 334               c = i+offset_lu;
 335               cdistx = distx[c];
 336               cdisty = disty[c];
 337               newdistx = cdistx+1;
 338               newdisty = cdisty+1;
 339               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 340               if(newdist < olddist-epsilon)
 341                 {
 342                   distx[i]=(short)newdistx;
 343                   disty[i]=(short)newdisty;
 344                             dist[i]=(short)newdist;
 345                   olddist=newdist;
 346                   changed = 1;
 347                 }
 348
 349               c = i+offset_u;
 350               cdistx = distx[c];
 351               cdisty = disty[c];
 352               newdistx = cdistx;
 353               newdisty = cdisty+1;
 354               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 355               if(newdist < olddist-epsilon)
 356                 {
 357                   distx[i]=(short)newdistx;
 358                   disty[i]=(short)newdisty;
 359                             dist[i]=(short)newdist;
 360                   changed = 1;
 361                 }
 362             }
 363
 364           /* Move index to second rightmost pixel of current row. */
 365           /* Rightmost pixel is skipped, it has no right neighbor. */
 366           i = y*w + w-2;
 367
 368           /* scan left, propagate distance from right */
 369           for(x=w-2; x>=0; x--, i--)
 370             {
 371               olddist = dist[i];
 372               if(olddist <= 0) continue; // Already zero distance
 373
 374               c = i+offset_r;
 375               cdistx = distx[c];
 376               cdisty = disty[c];
 377               newdistx = cdistx-1;
 378               newdisty = cdisty;
 379               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 380               if(newdist < olddist-epsilon)
 381                 {
 382                   distx[i]=(short)newdistx;
 383                   disty[i]=(short)newdisty;
 384                             dist[i]=(short)newdist;
 385                   changed = 1;
 386                 }
 387             }
 388         }
 389
 390       /* Scan rows in reverse order, except last row */
 391       for(y=h-2; y>=0; y--)
 392         {
 393           /* move index to rightmost pixel of current row */
 394           i = y*w + w-1;
 395
 396           /* Scan left, propagate distances from below & right */
 397
 398           /* Rightmost pixel is special, has no right neighbors */
 399           olddist = dist[i];
 400           if(olddist > 0) // If not already zero distance
 401             {
 402               c = i+offset_d;
 403               cdistx = distx[c];
 404               cdisty = disty[c];
 405               newdistx = cdistx;
 406               newdisty = cdisty-1;
 407               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 408               if(newdist < olddist-epsilon)
 409                 {
 410                   distx[i]=(short)newdistx;
 411                   disty[i]=(short)newdisty;
 412                             dist[i]=(short)newdist;
 413                   olddist=newdist;
 414                   changed = 1;
 415                 }
 416
 417               c = i+offset_dl;
 418               cdistx = distx[c];
 419               cdisty = disty[c];
 420               newdistx = cdistx+1;
 421               newdisty = cdisty-1;
 422               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 423               if(newdist < olddist-epsilon)
 424                 {
 425                   distx[i]=(short)newdistx;
 426                   disty[i]=(short)newdisty;
 427                             dist[i]=(short)newdist;
 428                   changed = 1;
 429                 }
 430             }
 431           i--;
 432
 433           /* Middle pixels have all neighbors */
 434           for(x=w-2; x>0; x--, i--)
 435             {
 436               olddist = dist[i];
 437               if(olddist <= 0) continue; // Already zero distance
 438
 439               c = i+offset_r;
 440               cdistx = distx[c];
 441               cdisty = disty[c];
 442               newdistx = cdistx-1;
 443               newdisty = cdisty;
 444               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 445               if(newdist < olddist-epsilon)
 446                 {
 447                   distx[i]=(short)newdistx;
 448                   disty[i]=(short)newdisty;
 449                   dist[i]=(short)newdist;
 450                   olddist=newdist;
 451                   changed = 1;
 452                 }
 453
 454               c = i+offset_rd;
 455               cdistx = distx[c];
 456               cdisty = disty[c];
 457               newdistx = cdistx-1;
 458               newdisty = cdisty-1;
 459               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 460               if(newdist < olddist-epsilon)
 461                 {
 462                   distx[i]=(short)newdistx;
 463                   disty[i]=(short)newdisty;
 464                   dist[i]=(short)newdist;
 465                   olddist=newdist;
 466                   changed = 1;
 467                 }
 468
 469               c = i+offset_d;
 470               cdistx = distx[c];
 471               cdisty = disty[c];
 472               newdistx = cdistx;
 473               newdisty = cdisty-1;
 474               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 475               if(newdist < olddist-epsilon)
 476                 {
 477                   distx[i]=(short)newdistx;
 478                   disty[i]=(short)newdisty;
 479                   dist[i]=(short)newdist;
 480                   olddist=newdist;
 481                   changed = 1;
 482                 }
 483
 484               c = i+offset_dl;
 485               cdistx = distx[c];
 486               cdisty = disty[c];
 487               newdistx = cdistx+1;
 488               newdisty = cdisty-1;
 489               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 490               if(newdist < olddist-epsilon)
 491                 {
 492                   distx[i]=(short)newdistx;
 493                   disty[i]=(short)newdisty;
 494                   dist[i]=(short)newdist;
 495                   changed = 1;
 496                 }
 497             }
 498           /* Leftmost pixel is special, has no left neighbors */
 499           olddist = dist[i];
 500           if(olddist > 0) // If not already zero distance
 501             {
 502               c = i+offset_r;
 503               cdistx = distx[c];
 504               cdisty = disty[c];
 505               newdistx = cdistx-1;
 506               newdisty = cdisty;
 507               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 508               if(newdist < olddist-epsilon)
 509                 {
 510                   distx[i]=(short)newdistx;
 511                   disty[i]=(short)newdisty;
 512                   dist[i]=(short)newdist;
 513                   olddist=newdist;
 514                   changed = 1;
 515                 }
 516
 517               c = i+offset_rd;
 518               cdistx = distx[c];
 519               cdisty = disty[c];
 520               newdistx = cdistx-1;
 521               newdisty = cdisty-1;
 522               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 523               if(newdist < olddist-epsilon)
 524                 {
 525                   distx[i]=(short)newdistx;
 526                   disty[i]=(short)newdisty;
 527                   dist[i]=(short)newdist;
 528                   olddist=newdist;
 529                   changed = 1;
 530                 }
 531
 532               c = i+offset_d;
 533               cdistx = distx[c];
 534               cdisty = disty[c];
 535               newdistx = cdistx;
 536               newdisty = cdisty-1;
 537               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 538               if(newdist < olddist-epsilon)
 539                 {
 540                   distx[i]=(short)newdistx;
 541                   disty[i]=(short)newdisty;
 542                   dist[i]=(short)newdist;
 543                   changed = 1;
 544                 }
 545             }
 546
 547           /* Move index to second leftmost pixel of current row. */
 548           /* Leftmost pixel is skipped, it has no left neighbor. */
 549           i = y*w + 1;
 550           for(x=1; x<w; x++, i++)
 551             {
 552               /* scan right, propagate distance from left */
 553               olddist = dist[i];
 554               if(olddist <= 0) continue; // Already zero distance
 555
 556               c = i+offset_l;
 557               cdistx = distx[c];
 558               cdisty = disty[c];
 559               newdistx = cdistx+1;
 560               newdisty = cdisty;
 561               newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
 562               if(newdist < olddist-epsilon)
 563                 {
 564                   distx[i]=(short)newdistx;
 565                   disty[i]=(short)newdisty;
 566                   dist[i]=(short)newdist;
 567                   changed = 1;
 568                 }
 569             }
 570         }
 571     }
 572   while(changed); // Sweep until no more updates are made
 573
 574   /* The transformation is completed. */
 575
 576 }