/*---------------------------------------------------------------------------*/ /* Program: distance.c */ /* */ /* Purpose: This program implements a four pass Euclidean distance */ /* transform (EDT) based on the algorithm in Mohamad Seaidoun's */ /* Ph.D. dissertation. All distances up to 100 are exact. */ /* The magic numbers K1..K5 and search offsets are derived */ /* by an exhaustive search for point configurations which */ /* may cause distance calculation errors (the real heart of */ /* the dissertation research focused on this process). */ /* */ /* Author: John Gauch and Mohamad Seaidoun */ /* */ /* Date: April 13, 1994 */ /* */ /* Note: Copyright (C) The University of Kansas, 1994 */ /*---------------------------------------------------------------------------*/ #include /* Neighborhood size constants */ #define K1 1 #define K2 169 /* really 13*13 */ #define K3 961 /* really 31*31 */ #define K4 2401 /* really 49*49 */ #define K5 5184 /* really 72*72 */ /* Global variables */ SHORT_TYPE **Data1; SHORT_TYPE **Xnear; SHORT_TYPE **Ynear; FLOAT_TYPE **Dist; int Xdim, Ydim, Zdim; int Debug = FALSE; int Exact = FALSE; int Pass = 4; double sqrt(); /*---------------------------------------------------------------------------*/ /* Purpose: This macro checks the specified neighbor. */ /*---------------------------------------------------------------------------*/ #define CHECK(dx,dy)\ {\ x2=x+(dx);\ y2=y+(dy);\ if ((x2>=0) && (y2>=0) && (x2<=Xdim-1) && (y2<=Ydim-1))\ if (Dist[y2][x2] < Dist[y][x])\ {\ Dx = Xnear[y2][x2] - x;\ Dy = Ynear[y2][x2] - y;\ NewDist = Dx*Dx + Dy*Dy;\ if (NewDist < Dist[y][x])\ {\ Xnear[y][x] = Xnear[y2][x2];\ Ynear[y][x] = Ynear[y2][x2];\ Dist[y][x] = (float)NewDist;\ }\ }\ } /*---------------------------------------------------------------------------*/ /* Purpose: This routine calculates the EDT of an image. */ /*---------------------------------------------------------------------------*/ void Distance() { int x, y, x2, y2, Dx, Dy, NewDist; /* Initialize output image */ if (Debug == TRUE) printf("Initialize\n"); for (y = 0; y < Ydim; y++) for (x = 0; x < Xdim; x++) { if (Data1[y][x] == 1) { Xnear[y][x] = x; Ynear[y][x] = y; Dist[y][x] = 0; } else { Xnear[y][x] = 0; Ynear[y][x] = 0; Dist[y][x] = (float) (Xdim * Xdim + Ydim * Ydim); } } /* Forward X, forward Y pass */ if (Debug == TRUE) printf("Forward X, forward Y pass\n"); if (Pass >= 1) for (y = 0; y < Ydim; y++) for (x = 0; x < Xdim; x++) { if (Dist[y][x] > K1) { CHECK(-1, 0); CHECK(-1, -1); CHECK(0, -1); } if (Exact == FALSE) continue; if (Dist[y][x] > K2) { CHECK(-2, -1); CHECK(-1, -2); } if (Dist[y][x] > K3) { CHECK(-3, -1); CHECK(-3, -2); CHECK(-2, -3); CHECK(-1, -3); } if (Dist[y][x] > K4) { CHECK(-4, -1); CHECK(-4, -3); CHECK(-3, -4); CHECK(-1, -4); } if (Dist[y][x] > K5) { CHECK(-5, -1); CHECK(-5, -2); CHECK(-5, -3); CHECK(-5, -4); CHECK(-1, -5); CHECK(-2, -5); CHECK(-3, -5); CHECK(-4, -5); } } /* Backward X, forward Y pass */ if (Debug == TRUE) printf("Backward X, forward Y pass\n"); if (Pass >= 2) for (y = 0; y < Ydim; y++) for (x = (Xdim - 1); x >= 0; x--) { if (Dist[y][x] > K1) { CHECK(1, 0); CHECK(1, -1); CHECK(0, -1); } if (Exact == FALSE) continue; if (Dist[y][x] > K2) { CHECK(2, -1); CHECK(1, -2); } if (Dist[y][x] > K3) { CHECK(3, -1); CHECK(3, -2); CHECK(2, -3); CHECK(1, -3); } if (Dist[y][x] > K4) { CHECK(4, -1); CHECK(4, -3); CHECK(3, -4); CHECK(1, -4); } if (Dist[y][x] > K5) { CHECK(5, -1); CHECK(5, -2); CHECK(5, -3); CHECK(5, -4); CHECK(1, -5); CHECK(2, -5); CHECK(3, -5); CHECK(4, -5); } } /* Forward X, backward Y pass */ if (Debug == TRUE) printf("Forward X, backward Y pass\n"); if (Pass >= 3) for (y = (Ydim - 1); y >= 0; y--) for (x = 0; x < Xdim; x++) { if (Dist[y][x] > K1) { CHECK(-1, 0); CHECK(-1, 1); CHECK(0, 1); } if (Exact == FALSE) continue; if (Dist[y][x] > K2) { CHECK(-2, 1); CHECK(-1, 2); } if (Dist[y][x] > K3) { CHECK(-3, 1); CHECK(-3, 2); CHECK(-2, 3); CHECK(-1, 3); } if (Dist[y][x] > K4) { CHECK(-4, 1); CHECK(-4, 3); CHECK(-3, 4); CHECK(-1, 4); } if (Dist[y][x] > K5) { CHECK(-5, 1); CHECK(-5, 2); CHECK(-5, 3); CHECK(-5, 4); CHECK(-1, 5); CHECK(-2, 5); CHECK(-3, 5); CHECK(-4, 5); } } /* Backward X, backward Y pass */ if (Debug == TRUE) printf("Backward X, backward Y pass\n"); if (Pass >= 4) for (y = (Ydim - 1); y >= 0; y--) for (x = (Xdim - 1); x >= 0; x--) { if (Dist[y][x] > K1) { CHECK(1, 0); CHECK(1, 1); CHECK(0, 1); } if (Exact == FALSE) continue; if (Dist[y][x] > K2) { CHECK(2, 1); CHECK(1, 2); } if (Dist[y][x] > K3) { CHECK(3, 1); CHECK(3, 2); CHECK(2, 3); CHECK(1, 3); } if (Dist[y][x] > K4) { CHECK(4, 1); CHECK(4, 3); CHECK(3, 4); CHECK(1, 4); } if (Dist[y][x] > K5) { CHECK(5, 1); CHECK(5, 2); CHECK(5, 3); CHECK(5, 4); CHECK(1, 5); CHECK(2, 5); CHECK(3, 5); CHECK(4, 5); } } /* Finalize output image */ for (y = 0; y < Ydim; y++) for (x = 0; x < Xdim; x++) Dist[y][x] = (float) sqrt((double) (Dist[y][x])); } /*---------------------------------------------------------------------------*/ /* Purpose: This is the main program. */ /*---------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { /* Image variables */ char Name1[50]; char Name2[50]; IM_TYPE *Image1; IM_TYPE *Image2; int PixType, DimCnt; /* Program variables */ int i = 0; /* Interpret program options */ printf("DISTANCE Program - KUIM Version 1.0\n\n"); while ((++i < argc) && (argv[i][0] == '-')) switch (argv[i][1]) { case 'p': if (sscanf(argv[++i], "%d", &Pass) == 0) Error("Could not get number of passes"); break; case 'd': Debug = TRUE; break; case 'e': Exact = TRUE; break; default: Error("Invalid option encountered"); break; } /* Check number of file names */ if (argc - i != 2) { fprintf(stderr, "Usage: distance [options] infile outfile\n"); fprintf(stderr, " [-d] Print debugging information\n"); fprintf(stderr, " [-e] Calculate exact distances\n"); fprintf(stderr, " [-p #] Number of passes (default=4)\n"); exit(1); } /* Get image file names from argument list */ if (sscanf(argv[i++], "%s", Name1) == 0) Error("Could not get input file name"); if (sscanf(argv[i++], "%s", Name2) == 0) Error("Could not get output file name"); /* Read input image */ Image1 = im_open(Name1, &PixType, &Xdim, &Ydim, &Zdim, &DimCnt); if (DimCnt != 2) Error("Can not process 1D or 3D images"); Data1 = (SHORT_TYPE **) im_alloc2D(Image1, SHORT); im_read(Image1, SHORT, (char *) &(Data1[0][0])); /* Create output image */ Image2 = im_create(Name2, FLOAT, Xdim, Ydim, Zdim); Xnear = (SHORT_TYPE **) im_alloc2D(Image2, SHORT); Ynear = (SHORT_TYPE **) im_alloc2D(Image2, SHORT); Dist = (FLOAT_TYPE **) im_alloc2D(Image2, FLOAT); /* Calculate distance */ Distance(); /* Write information to output image */ im_write(Image2, FLOAT, (char *) &(Dist[0][0])); im_free2D((char **) Data1); im_free2D((char **) Xnear); im_free2D((char **) Ynear); im_free2D((char **) Dist); return (0); }