123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582 |
- #ifdef __cplusplus
- extern "C" {
- #endif
- #include <math.h>
- void computegradient(double *img, int w, int h, double *gx, double *gy)
- {
- int i,j,k;
- double glength;
- #define SQRT2 1.4142136
- for(i = 1; i < h-1; i++) {
- for(j = 1; j < w-1; j++) {
- k = i*w + j;
- if((img[k]>0.0) && (img[k]<1.0)) {
- 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];
- 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];
- glength = gx[k]*gx[k] + gy[k]*gy[k];
- if(glength > 0.0) {
- glength = sqrt(glength);
- gx[k]=gx[k]/glength;
- gy[k]=gy[k]/glength;
- }
- }
- }
- }
-
-
-
- }
- double edgedf(double gx, double gy, double a)
- {
- double df, glength, temp, a1;
- if ((gx == 0) || (gy == 0)) {
- df = 0.5-a;
- } else {
- glength = sqrt(gx*gx + gy*gy);
- if(glength>0) {
- gx = gx/glength;
- gy = gy/glength;
- }
-
- gx = fabs(gx);
- gy = fabs(gy);
- if(gx<gy) {
- temp = gx;
- gx = gy;
- gy = temp;
- }
- a1 = 0.5*gy/gx;
- if (a < a1) {
- df = 0.5*(gx + gy) - sqrt(2.0*gx*gy*a);
- } else if (a < (1.0-a1)) {
- df = (0.5-a)*gx;
- } else {
- df = -0.5*(gx + gy) + sqrt(2.0*gx*gy*(1.0-a));
- }
- }
- return df;
- }
- double distaa3(double *img, double *gximg, double *gyimg, int w, int c, int xc, int yc, int xi, int yi)
- {
- double di, df, dx, dy, gx, gy, a;
- int closest;
-
- closest = c-xc-yc*w;
- a = img[closest];
- gx = gximg[closest];
- gy = gyimg[closest];
-
- if(a > 1.0) a = 1.0;
- if(a < 0.0) a = 0.0;
- if(a == 0.0) return 1000000.0;
- dx = (double)xi;
- dy = (double)yi;
- di = sqrt(dx*dx + dy*dy);
- if(di==0) {
-
- df = edgedf(gx, gy, a);
- } else {
-
- df = edgedf(dx, dy, a);
- }
- return di + df;
- }
- #define DISTAA(c,xc,yc,xi,yi) (distaa3(img, gx, gy, w, c, xc, yc, xi, yi))
- void edtaa3(double *img, double *gx, double *gy, int w, int h, short *distx, short *disty, double *dist)
- {
- int x, y, i, c;
- int offset_u, offset_ur, offset_r, offset_rd,
- offset_d, offset_dl, offset_l, offset_lu;
- double olddist, newdist;
- int cdistx, cdisty, newdistx, newdisty;
- int changed;
- double epsilon = 1e-3;
-
- offset_u = -w;
- offset_ur = -w+1;
- offset_r = 1;
- offset_rd = w+1;
- offset_d = w;
- offset_dl = w-1;
- offset_l = -1;
- offset_lu = -w-1;
-
- for(i=0; i<w*h; i++) {
- distx[i] = 0;
- disty[i] = 0;
- if(img[i] <= 0.0)
- {
- dist[i]= 1000000.0;
- }
- else if (img[i]<1.0) {
- dist[i] = edgedf(gx[i], gy[i], img[i]);
- }
- else {
- dist[i]= 0.0;
- }
- }
-
- do
- {
- changed = 0;
-
- for(y=1; y<h; y++)
- {
-
- i = y*w;
-
-
- olddist = dist[i];
- if(olddist > 0)
- {
- c = i + offset_u;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx;
- newdisty = cdisty+1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_ur;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx-1;
- newdisty = cdisty+1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
- i++;
-
- for(x=1; x<w-1; x++, i++)
- {
- olddist = dist[i];
- if(olddist <= 0) continue;
- c = i+offset_l;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx+1;
- newdisty = cdisty;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_lu;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx+1;
- newdisty = cdisty+1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_u;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx;
- newdisty = cdisty+1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_ur;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx-1;
- newdisty = cdisty+1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
-
- olddist = dist[i];
- if(olddist > 0)
- {
- c = i+offset_l;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx+1;
- newdisty = cdisty;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_lu;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx+1;
- newdisty = cdisty+1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_u;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx;
- newdisty = cdisty+1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
-
-
- i = y*w + w-2;
-
- for(x=w-2; x>=0; x--, i--)
- {
- olddist = dist[i];
- if(olddist <= 0) continue;
- c = i+offset_r;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx-1;
- newdisty = cdisty;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
- }
-
-
- for(y=h-2; y>=0; y--)
- {
-
- i = y*w + w-1;
-
-
- olddist = dist[i];
- if(olddist > 0)
- {
- c = i+offset_d;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx;
- newdisty = cdisty-1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_dl;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx+1;
- newdisty = cdisty-1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
- i--;
-
- for(x=w-2; x>0; x--, i--)
- {
- olddist = dist[i];
- if(olddist <= 0) continue;
- c = i+offset_r;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx-1;
- newdisty = cdisty;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_rd;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx-1;
- newdisty = cdisty-1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_d;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx;
- newdisty = cdisty-1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_dl;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx+1;
- newdisty = cdisty-1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
-
- olddist = dist[i];
- if(olddist > 0)
- {
- c = i+offset_r;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx-1;
- newdisty = cdisty;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_rd;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx-1;
- newdisty = cdisty-1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- olddist=newdist;
- changed = 1;
- }
- c = i+offset_d;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx;
- newdisty = cdisty-1;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
-
-
- i = y*w + 1;
- for(x=1; x<w; x++, i++)
- {
-
- olddist = dist[i];
- if(olddist <= 0) continue;
- c = i+offset_l;
- cdistx = distx[c];
- cdisty = disty[c];
- newdistx = cdistx+1;
- newdisty = cdisty;
- newdist = DISTAA(c, cdistx, cdisty, newdistx, newdisty);
- if(newdist < olddist-epsilon)
- {
- distx[i]=newdistx;
- disty[i]=newdisty;
- dist[i]=newdist;
- changed = 1;
- }
- }
- }
- }
- while(changed);
-
- }
- #ifdef __cplusplus
- }
- #endif
|