/* * The Npic library * * Copyright (C) 2003 Edouard Thiel * * This library is free software under the terms of the * GNU Lesser General Public License (LGPL) version 2.1. */ /* DO NOT EDIT !!! Generated by npic-templa from "dist_wdt.ct" */ /* * dist_wdt.c - 13/01/2004 * * Weighted distance Transforms (a.k.a chamfer DTs). * * References: * * [1] A. Rosenfeld and J.L. Pfaltz, 1966. Sequential operations in digital * picture processing. Journal of ACM, 13(4):471-494. * * [2] U. Montanari, 1968. A method for obtaining skeletons using a * quasi-euclidean distance. Journal of ACM, 15:600-624. * * [3] G. Borgefors, 1984. Distance transformations in arbitrary dimensions. * Computer Vision, Graphics and Image Processing, 27:321-345. */ #include /*--------------------- P U B L I C - I N T E R F A C E ----------------------*/ /* * Compute the weighted distance transform (WDT) on image np * for the half chamfer mask mh. * * Do nothing if np is not ok. set not ok on error. Verbose. * Return computation time in second.microsecond */ double NpicWDT (Npic_image *np, Npic_mask *mh) { int maxcoord; double time1, time2; if (NpicWDTCheckData (np, mh, __func__) != NPIC_SUCCESS) return 0; /* Compute largest coordinate in mask */ maxcoord = NpicMaskLargestCoord (mh); /* Change the border width of np to maxcoord if it was smaller */ NpicSetBorderWidthMin (np, maxcoord); /* Set border pixels to 0 */ switch (np->type) { case NPIC_IMAGE_2L : case NPIC_IMAGE_3L : case NPIC_IMAGE_4L : case NPIC_IMAGE_5L : case NPIC_IMAGE_6L : NpicFillBorder_i (np, 0); break; case NPIC_IMAGE_2D : case NPIC_IMAGE_3D : case NPIC_IMAGE_4D : case NPIC_IMAGE_5D : case NPIC_IMAGE_6D : NpicFillBorder_d (np, 0); break; } if (np->gen.ok != NPIC_SUCCESS) return 0; time1 = NpicGetTime(); /* Forward and backward scans */ switch (np->type) { case NPIC_IMAGE_2L : { Npic_image_2l *p = NpicCastImage(np); Npic_mask_2l *m = NpicCastMask(mh); Npic_l a, b; int i, y, x; /* Forward scan */ for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[y][x] != 0) { a = p->pix[y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[y][x] = a; } /* Backward scan */ for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[y][x] != 0) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[y][x] = a; } break; } case NPIC_IMAGE_2D : { Npic_image_2d *p = NpicCastImage(np); Npic_mask_2d *m = NpicCastMask(mh); Npic_d a, b; int i, y, x; /* Forward scan */ for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[y][x] != 0) { a = p->pix[y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[y][x] = a; } /* Backward scan */ for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[y][x] != 0) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[y][x] = a; } break; } case NPIC_IMAGE_3L : { Npic_image_3l *p = NpicCastImage(np); Npic_mask_3l *m = NpicCastMask(mh); Npic_l a, b; int i, z, y, x; /* Forward scan */ for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[z][y][x] != 0) { a = p->pix[z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[z][y][x] = a; } /* Backward scan */ for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[z][y][x] != 0) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[z][y][x] = a; } break; } case NPIC_IMAGE_3D : { Npic_image_3d *p = NpicCastImage(np); Npic_mask_3d *m = NpicCastMask(mh); Npic_d a, b; int i, z, y, x; /* Forward scan */ for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[z][y][x] != 0) { a = p->pix[z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[z][y][x] = a; } /* Backward scan */ for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[z][y][x] != 0) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[z][y][x] = a; } break; } case NPIC_IMAGE_4L : { Npic_image_4l *p = NpicCastImage(np); Npic_mask_4l *m = NpicCastMask(mh); Npic_l a, b; int i, t, z, y, x; /* Forward scan */ for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[t][z][y][x] != 0) { a = p->pix[t - m->v[0].t][z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[t][z][y][x] = a; } /* Backward scan */ for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[t][z][y][x] != 0) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[t][z][y][x] = a; } break; } case NPIC_IMAGE_4D : { Npic_image_4d *p = NpicCastImage(np); Npic_mask_4d *m = NpicCastMask(mh); Npic_d a, b; int i, t, z, y, x; /* Forward scan */ for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[t][z][y][x] != 0) { a = p->pix[t - m->v[0].t][z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[t][z][y][x] = a; } /* Backward scan */ for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[t][z][y][x] != 0) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[t][z][y][x] = a; } break; } case NPIC_IMAGE_5L : { Npic_image_5l *p = NpicCastImage(np); Npic_mask_5l *m = NpicCastMask(mh); Npic_l a, b; int i, s, t, z, y, x; /* Forward scan */ for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[s][t][z][y][x] != 0) { a = p->pix[s - m->v[0].s] [t - m->v[0].t][z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[s][t][z][y][x] = a; } /* Backward scan */ for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[s][t][z][y][x] != 0) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[s][t][z][y][x] = a; } break; } case NPIC_IMAGE_5D : { Npic_image_5d *p = NpicCastImage(np); Npic_mask_5d *m = NpicCastMask(mh); Npic_d a, b; int i, s, t, z, y, x; /* Forward scan */ for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[s][t][z][y][x] != 0) { a = p->pix[s - m->v[0].s] [t - m->v[0].t][z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[s][t][z][y][x] = a; } /* Backward scan */ for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[s][t][z][y][x] != 0) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[s][t][z][y][x] = a; } break; } case NPIC_IMAGE_6L : { Npic_image_6l *p = NpicCastImage(np); Npic_mask_6l *m = NpicCastMask(mh); Npic_l a, b; int i, r, s, t, z, y, x; /* Forward scan */ for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[r][s][t][z][y][x] != 0) { a = p->pix[r - m->v[0].r][s - m->v[0].s] [t - m->v[0].t][z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[r - m->v[i].r][s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[r][s][t][z][y][x] = a; } /* Backward scan */ for (r = p->rmax-1; r >= 0 ; r--) for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[r][s][t][z][y][x] != 0) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r + m->v[i].r][s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[r][s][t][z][y][x] = a; } break; } case NPIC_IMAGE_6D : { Npic_image_6d *p = NpicCastImage(np); Npic_mask_6d *m = NpicCastMask(mh); Npic_d a, b; int i, r, s, t, z, y, x; /* Forward scan */ for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[r][s][t][z][y][x] != 0) { a = p->pix[r - m->v[0].r][s - m->v[0].s] [t - m->v[0].t][z - m->v[0].z] [y - m->v[0].y][x - m->v[0].x] + m->v[0].h; for (i = 1; i < m->nb; i++) { b = p->pix[r - m->v[i].r][s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[r][s][t][z][y][x] = a; } /* Backward scan */ for (r = p->rmax-1; r >= 0 ; r--) for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[r][s][t][z][y][x] != 0) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r + m->v[i].r][s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] + m->v[i].h; if (b < a) a = b; } p->pix[r][s][t][z][y][x] = a; } break; } default : np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, ""); } time2 = NpicGetTime(); return time2-time1; } /* * Compute the weighted distance transform (WDT_inf) on image np * for the half chamfer mask mh, with "infinite border", i.e border pixels * propagating no distance value. * * Do nothing if np is not ok. set not ok on error. Verbose. * Return computation time in second.microsecond */ double NpicWDT_inf (Npic_image *np, Npic_mask *mh) { int maxcoord; double time1, time2; if (NpicWDTCheckData (np, mh, __func__) != NPIC_SUCCESS) return 0; /* Compute largest coordinate in mask */ maxcoord = NpicMaskLargestCoord (mh); /* Change the border width of np to maxcoord if it was smaller */ NpicSetBorderWidthMin (np, maxcoord); /* Set border pixels to -1 */ switch (np->type) { case NPIC_IMAGE_2L : case NPIC_IMAGE_3L : case NPIC_IMAGE_4L : case NPIC_IMAGE_5L : case NPIC_IMAGE_6L : NpicFillBorder_i (np, -1); break; case NPIC_IMAGE_2D : case NPIC_IMAGE_3D : case NPIC_IMAGE_4D : case NPIC_IMAGE_5D : case NPIC_IMAGE_6D : NpicFillBorder_d (np, -1); break; } if (np->gen.ok != NPIC_SUCCESS) return 0; time1 = NpicGetTime(); /* Forward and backward scans */ switch (np->type) { case NPIC_IMAGE_2L : { Npic_image_2l *p = NpicCastImage(np); Npic_mask_2l *m = NpicCastMask(mh); Npic_l a, b, c; int i, y, x; /* Forward scan */ for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[y][x] = a; } /* Backward scan */ for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[y][x] != 0) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[y][x] = a; } break; } case NPIC_IMAGE_2D : { Npic_image_2d *p = NpicCastImage(np); Npic_mask_2d *m = NpicCastMask(mh); Npic_d a, b, c; int i, y, x; /* Forward scan */ for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[y][x] = a; } /* Backward scan */ for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[y][x] != 0) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[y][x] = a; } break; } case NPIC_IMAGE_3L : { Npic_image_3l *p = NpicCastImage(np); Npic_mask_3l *m = NpicCastMask(mh); Npic_l a, b, c; int i, z, y, x; /* Forward scan */ for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[z][y][x] = a; } /* Backward scan */ for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[z][y][x] != 0) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[z][y][x] = a; } break; } case NPIC_IMAGE_3D : { Npic_image_3d *p = NpicCastImage(np); Npic_mask_3d *m = NpicCastMask(mh); Npic_d a, b, c; int i, z, y, x; /* Forward scan */ for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[z][y][x] = a; } /* Backward scan */ for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[z][y][x] != 0) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[z][y][x] = a; } break; } case NPIC_IMAGE_4L : { Npic_image_4l *p = NpicCastImage(np); Npic_mask_4l *m = NpicCastMask(mh); Npic_l a, b, c; int i, t, z, y, x; /* Forward scan */ for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[t][z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[t][z][y][x] = a; } /* Backward scan */ for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[t][z][y][x] != 0) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[t][z][y][x] = a; } break; } case NPIC_IMAGE_4D : { Npic_image_4d *p = NpicCastImage(np); Npic_mask_4d *m = NpicCastMask(mh); Npic_d a, b, c; int i, t, z, y, x; /* Forward scan */ for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[t][z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[t][z][y][x] = a; } /* Backward scan */ for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[t][z][y][x] != 0) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[t][z][y][x] = a; } break; } case NPIC_IMAGE_5L : { Npic_image_5l *p = NpicCastImage(np); Npic_mask_5l *m = NpicCastMask(mh); Npic_l a, b, c; int i, s, t, z, y, x; /* Forward scan */ for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[s][t][z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[s][t][z][y][x] = a; } /* Backward scan */ for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[s][t][z][y][x] != 0) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[s][t][z][y][x] = a; } break; } case NPIC_IMAGE_5D : { Npic_image_5d *p = NpicCastImage(np); Npic_mask_5d *m = NpicCastMask(mh); Npic_d a, b, c; int i, s, t, z, y, x; /* Forward scan */ for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[s][t][z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[s][t][z][y][x] = a; } /* Backward scan */ for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[s][t][z][y][x] != 0) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[s][t][z][y][x] = a; } break; } case NPIC_IMAGE_6L : { Npic_image_6l *p = NpicCastImage(np); Npic_mask_6l *m = NpicCastMask(mh); Npic_l a, b, c; int i, r, s, t, z, y, x; /* Forward scan */ for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[r][s][t][z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[r - m->v[i].r][s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[r][s][t][z][y][x] = a; } /* Backward scan */ for (r = p->rmax-1; r >= 0 ; r--) for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[r][s][t][z][y][x] != 0) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r + m->v[i].r][s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[r][s][t][z][y][x] = a; } break; } case NPIC_IMAGE_6D : { Npic_image_6d *p = NpicCastImage(np); Npic_mask_6d *m = NpicCastMask(mh); Npic_d a, b, c; int i, r, s, t, z, y, x; /* Forward scan */ for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) if (p->pix[r][s][t][z][y][x] != 0) { a = -1; for (i = 0; i < m->nb; i++) { b = p->pix[r - m->v[i].r][s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[r][s][t][z][y][x] = a; } /* Backward scan */ for (r = p->rmax-1; r >= 0 ; r--) for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) if (p->pix[r][s][t][z][y][x] != 0) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r + m->v[i].r][s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x]; if (b >= 0) { c = b + m->v[i].h; if (a < 0 || c < a) a = c; } } p->pix[r][s][t][z][y][x] = a; } break; } default : np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, ""); } time2 = NpicGetTime(); return time2-time1; } /* * Compute the reverse weighted distance transform (WDT_rev) on image np * for the half chamfer mask mh. * * Do nothing if np is not ok. set not ok on error. Verbose. * Return computation time in second.microsecond */ double NpicWDT_rev (Npic_image *np, Npic_mask *mh) { int maxcoord; double time1, time2; if (NpicWDTCheckData (np, mh, __func__) != NPIC_SUCCESS) return 0; /* Compute largest coordinate in mask */ maxcoord = NpicMaskLargestCoord (mh); /* Change the border width of np to maxcoord if it was smaller */ NpicSetBorderWidthMin (np, maxcoord); /* Set border pixels to 0 */ switch (np->type) { case NPIC_IMAGE_2L : case NPIC_IMAGE_3L : case NPIC_IMAGE_4L : case NPIC_IMAGE_5L : case NPIC_IMAGE_6L : NpicFillBorder_i (np, 0); break; case NPIC_IMAGE_2D : case NPIC_IMAGE_3D : case NPIC_IMAGE_4D : case NPIC_IMAGE_5D : case NPIC_IMAGE_6D : NpicFillBorder_d (np, 0); break; } if (np->gen.ok != NPIC_SUCCESS) return 0; time1 = NpicGetTime(); /* Forward and backward scans */ switch (np->type) { case NPIC_IMAGE_2L : { Npic_image_2l *p = NpicCastImage(np); Npic_mask_2l *m = NpicCastMask(mh); Npic_l a, b; int i, y, x; /* Forward scan */ for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[y][x] = a; } /* Backward scan */ for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[y][x] = a; } break; } case NPIC_IMAGE_2D : { Npic_image_2d *p = NpicCastImage(np); Npic_mask_2d *m = NpicCastMask(mh); Npic_d a, b; int i, y, x; /* Forward scan */ for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[y][x] = a; } /* Backward scan */ for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[y][x] = a; } break; } case NPIC_IMAGE_3L : { Npic_image_3l *p = NpicCastImage(np); Npic_mask_3l *m = NpicCastMask(mh); Npic_l a, b; int i, z, y, x; /* Forward scan */ for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[z][y][x] = a; } /* Backward scan */ for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[z][y][x] = a; } break; } case NPIC_IMAGE_3D : { Npic_image_3d *p = NpicCastImage(np); Npic_mask_3d *m = NpicCastMask(mh); Npic_d a, b; int i, z, y, x; /* Forward scan */ for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[z][y][x] = a; } /* Backward scan */ for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[z][y][x] = a; } break; } case NPIC_IMAGE_4L : { Npic_image_4l *p = NpicCastImage(np); Npic_mask_4l *m = NpicCastMask(mh); Npic_l a, b; int i, t, z, y, x; /* Forward scan */ for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[t][z][y][x] = a; } /* Backward scan */ for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[t][z][y][x] = a; } break; } case NPIC_IMAGE_4D : { Npic_image_4d *p = NpicCastImage(np); Npic_mask_4d *m = NpicCastMask(mh); Npic_d a, b; int i, t, z, y, x; /* Forward scan */ for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[t][z][y][x] = a; } /* Backward scan */ for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[t][z][y][x] = a; } break; } case NPIC_IMAGE_5L : { Npic_image_5l *p = NpicCastImage(np); Npic_mask_5l *m = NpicCastMask(mh); Npic_l a, b; int i, s, t, z, y, x; /* Forward scan */ for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[s][t][z][y][x] = a; } /* Backward scan */ for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[s][t][z][y][x] = a; } break; } case NPIC_IMAGE_5D : { Npic_image_5d *p = NpicCastImage(np); Npic_mask_5d *m = NpicCastMask(mh); Npic_d a, b; int i, s, t, z, y, x; /* Forward scan */ for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[s][t][z][y][x] = a; } /* Backward scan */ for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[s][t][z][y][x] = a; } break; } case NPIC_IMAGE_6L : { Npic_image_6l *p = NpicCastImage(np); Npic_mask_6l *m = NpicCastMask(mh); Npic_l a, b; int i, r, s, t, z, y, x; /* Forward scan */ for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r - m->v[i].r][s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[r][s][t][z][y][x] = a; } /* Backward scan */ for (r = p->rmax-1; r >= 0 ; r--) for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r + m->v[i].r][s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[r][s][t][z][y][x] = a; } break; } case NPIC_IMAGE_6D : { Npic_image_6d *p = NpicCastImage(np); Npic_mask_6d *m = NpicCastMask(mh); Npic_d a, b; int i, r, s, t, z, y, x; /* Forward scan */ for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r - m->v[i].r][s - m->v[i].s] [t - m->v[i].t][z - m->v[i].z] [y - m->v[i].y][x - m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[r][s][t][z][y][x] = a; } /* Backward scan */ for (r = p->rmax-1; r >= 0 ; r--) for (s = p->smax-1; s >= 0 ; s--) for (t = p->tmax-1; t >= 0 ; t--) for (z = p->zmax-1; z >= 0 ; z--) for (y = p->ymax-1; y >= 0 ; y--) for (x = p->xmax-1; x >= 0 ; x--) { a = p->pix[r][s][t][z][y][x]; for (i = 0; i < m->nb; i++) { b = p->pix[r + m->v[i].r][s + m->v[i].s] [t + m->v[i].t][z + m->v[i].z] [y + m->v[i].y][x + m->v[i].x] - m->v[i].h; if (b > a) a = b; } p->pix[r][s][t][z][y][x] = a; } break; } default : np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, ""); } time2 = NpicGetTime(); return time2-time1; } /*-------------------- P R I V A T E - F U N C T I O N S ---------------------*/ /* * Check if images np and mask mh are ok and are compatible for DT and RDT * computations. * Return NPIC_SUCCESS, else error code. * set not ok for np on error. Verbose. */ int NpicWDTCheckData (Npic_image *np, Npic_mask *mh, const char *funcname) { if (NpicImageIsOK (np, funcname) != NPIC_SUCCESS) return NPIC_ERR_NOT_OK; if (NpicMaskIsOK (mh, funcname) != NPIC_SUCCESS) return np->gen.ok = NPIC_ERR_NOT_OK; if (NpicMaskCompat (mh, np) != NPIC_SUCCESS) return np->gen.ok = NpicError (funcname, NPIC_ERR_INCOMPAT, ": dest image and dist mask"); if (mh->gen.nb == 0) return np->gen.ok = NpicError (funcname, NPIC_ERROR, ": empty mask"); return NPIC_SUCCESS; } /*----------------------------------------------------------------------------*/