/* cleanskel.c Written by by Zinovy Pugach November, 1997. This function takes in a skeletonized image (of a digit or a character) and removes imperfections (short lines that are stemming from the actual skeleton) that may be caused by imperfection in the sceletonizing function. The input arguments are the uint8 black and white image where the skeleton is white, and the tailsize that is considered as imperfection (everithing shorter than tailsize stemming from anything will be removed). This procedure is repeated until no changes to the picture are made -- in this fashion a Y-stem coming from the skeleton will be removed. For some strange reason, the returned image is not recognized as bwimage. To correct this problem, negate(~) the matrix twice. */ #include #include #include "mex.h" #include "matrix.h" int deltail(int x, int y, int oldx, int oldy); int cleanup(int x, int y, int adjpix); int adjpixels(int x, int y); void findpath(int *x, int *y, int oldx, int oldy); #define min(a,b) ((a < b) ? a : b) #define max(a,b) ((a > b) ? a : b) int ma, na, tailsize; int depth; char *ap; char *temp; void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int changed = 1; int x,y,adjpix; int dims[2]; char *lptr, *rptr; if (nrhs < 2) { mexErrMsgTxt("Usage: newimage = cleanskel(bwimage,tailsize)"); } ma = mxGetM(prhs[0]); na = mxGetN(prhs[0]); rptr = (char *) mxGetPr(prhs[0]); ap = (char *) calloc(ma * na, sizeof(char)); memcpy(ap, rptr, ma * na * sizeof(char)); /* pointer to the bw picture */ tailsize = *mxGetPr(prhs[1]); temp = (char *) calloc(ma * na, sizeof(char)); memcpy(temp, ap, ma * na * sizeof(char)); while (changed) { changed = 0; for(y = 0; y < ma; y++) { for(x = 0; x < na; x++) { if(ap[(x * ma) + y]) { adjpix=adjpixels(x,y); if (adjpix < 2) { depth = 0; /* printf("End point at (%d,%d)\n",x,y); */ if (deltail(x,y, -1, -1)) { /* printf("Tail erased at %d, %d\n",x,y); */ changed = 1; } } else if (adjpix < 4) /* adjpix = 2 or 3 */ if (cleanup(x,y,adjpix)) changed = 1; } } } memcpy(ap, temp, ma * na * sizeof(char)); /* refresh ap picture */ } dims[1] = na; dims[0] = ma; plhs[0] = mxCreateNumericArray(2, dims, mxUINT8_CLASS, mxREAL); lptr = (char *) mxGetPr(plhs[0]); memcpy(lptr, temp, ma * na * sizeof(char)); free(temp); free(ap); } int deltail(int x, int y, int oldx, int oldy) { int xn, yn, adjpix; depth++; /* printf(" %d: (%d,%d) ", depth,x,y); */ if (depth < tailsize) { adjpix = adjpixels(x,y); if (adjpix > 2) { return 1; } else { if (adjpix == 0) { temp[x*ma + y] = 0; return 1; } else { if ((adjpix == 1) && (oldx != -1)) { temp[x*ma + y] = 0; return 1; } else { xn = x; yn = y; findpath(&xn,&yn, oldx, oldy); if (deltail(xn,yn,x,y)) { temp[x*ma + y] = 0; return 1; } else return 0; } } } } else return 0; } int cleanup(int x, int y, int adjpix) /* 000 000 010 ==> 000 011 011 000 000 011 ==> 001 010 010 000 000 111 ==> 101 010 010 000 000 010 ==> 000 111 111 000 000 011 ==> 001 011 011 */ { if (2 == adjpix) { if (temp[(x-1)*ma+y]) { if (temp[x*ma+(y-1)] || temp[x*ma+(y+1)]) /* corner */ temp[x*ma+y] = 0; if (temp[(x-1)*ma+(y-1)] || temp[(x-1)*ma+(y+1)]) /* stump */ temp[x*ma+y] = 0; } if (temp[(x+1)*ma+y]) { if (temp[x*ma+(y-1)] || temp[x*ma+(y+1)]) temp[x*ma+y] = 0; if (temp[(x+1)*ma+(y-1)] || temp[(x+1)*ma+(y+1)]) temp[x*ma+y] = 0; } if (temp[x*ma+(y-1)]) { if (temp[(x-1)*ma+y] || temp[(x+1)*ma+y]) temp[x*ma+y] = 0; if (temp[(x-1)*ma+(y-1)] || temp[(x+1)*ma+(y-1)]) temp[x*ma+y] = 0; } if (temp[x*ma+(y+1)]) { if (temp[(x-1)*ma+y] || temp[(x+1)*ma+y]) temp[x*ma+y] = 0; if (temp[(x-1)*ma+(y+1)] || temp[(x+1)*ma+(y+1)]) temp[x*ma+y] = 0; } } if (3 == adjpix) { if (temp[x*ma+(y-1)] && temp[x*ma+(y+1)] && (temp[(x-1)*ma+y] || temp[(x+1)*ma+y])) temp[x*ma+y] = 0; if (temp[(x-1)*ma+y] && temp[(x+1)*ma+y] && (temp[x*ma+(y-1)] || temp[x*ma+(y+1)])) temp[x*ma+y] = 0; if (temp[(x-1)*ma+(y-1)] && temp[x*ma+(y-1)] && temp[(x+1)*ma+(y-1)]) temp[x*ma+y] = 0; if (temp[(x-1)*ma+(y+1)] && temp[x*ma+(y+1)] && temp[(x+1)*ma+(y+1)]) temp[x*ma+y] = 0; if (temp[(x-1)*ma+(y-1)] && temp[(x-1)*ma+y] && temp[(x-1)*ma+(y+1)]) temp[x*ma+y] = 0; if (temp[(x+1)*ma+(y-1)] && temp[(x+1)*ma+y] && temp[(x+1)*ma+(y+1)]) temp[x*ma+y] = 0; if (temp[(x+1)*ma+y] && temp[(x+1)*ma+(y+1)] && temp[x*ma+(y+1)]) temp[x*ma+y] = 0; if (temp[(x-1)*ma+y] && temp[(x-1)*ma+(y+1)] && temp[x*ma+(y+1)]) temp[x*ma+y] = 0; if (temp[(x-1)*ma+y] && temp[(x-1)*ma+(y-1)] && temp[x*ma+(y-1)]) temp[x*ma+y] = 0; if (temp[(x+1)*ma+y] && temp[(x+1)*ma+(y-1)] && temp[x*ma+(y-1)]) temp[x*ma+y] = 0; } return !(temp[x*ma+y]); } int adjpixels(int x, int y) { int a; a = ap[((x-1)*ma)+(y-1)] + ap[((x+0)*ma)+(y-1)] + ap[((x+1)*ma)+(y-1)] + ap[((x-1)*ma)+(y+1)] + ap[((x+0)*ma)+(y+1)] + ap[((x+1)*ma)+(y+1)] + ap[((x-1)*ma)+(y+0)] + ap[((x+1)*ma)+(y-0)]; return a; } void findpath(int *x, int *y, int oldx, int oldy) { if (ap[((*x-1)*ma)+(*y-1)]) if (((*x-1) != oldx) || ((*y-1) != oldy)) { *x = *x - 1; *y = *y - 1; return; } if (ap[((*x+0)*ma)+(*y-1)]) if (((*x+0) != oldx) || ((*y-1) != oldy)) { *x = *x + 0; *y = *y - 1; return; } if (ap[((*x+1)*ma)+(*y-1)]) if (((*x+1) != oldx) || ((*y-1) != oldy)) { *x = *x + 1; *y = *y - 1; return; } if (ap[((*x-1)*ma)+(*y+1)]) if (((*x-1) != oldx) || ((*y+1) != oldy)) { *x = *x - 1; *y = *y + 1; return; } if (ap[((*x+0)*ma)+(*y+1)]) if (((*x+0) != oldx) || ((*y+1) != oldy)) { *x = *x + 0; *y = *y + 1; return; } if (ap[((*x+1)*ma)+(*y+1)]) if (((*x+1) != oldx) || ((*y+1) != oldy)) { *x = *x + 1; *y = *y + 1; return; } if (ap[((*x-1)*ma)+(*y+0)]) if (((*x-1) != oldx) || ((*y+0) != oldy)) { *x = *x - 1; *y = *y + 0; return; } if (ap[((*x+1)*ma)+(*y+0)]) if (((*x+1) != oldx) || ((*y+0) != oldy)) { *x = *x + 1; *y = *y + 0; return; } printf("Function findpath is screwing up!\n"); }