/* * 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 "mlut_sedtg.ct" */ /* * mlut_sedtg.c - 01/11/2008 * * SEDT in G(Z^n) for Medial Axis extraction. * * References: * * [1] E. Remy and E. Thiel, 2005. Exact Medial Axis with Euclidean Distance. * Image and Vision Computing, 23(2):167-175. * * [2] T. Hirata, 1996. "A unified linear-time algorithm for computing * distance maps". Information Processing Letters, 58:129-133. */ #include /*--------------------- P U B L I C - I N T E R F A C E ----------------------*/ /* * Fast Cone Distance Transform Algorithm for SED. * * Input: nCTg is a L^n image. * Output: nCTg the L^n distance image to the origin. * * Do nothing if nCTg is not ok. set not ok on error. Verbose. */ void NpicMlutCTg_SED (Npic_image *nCTg) { if (NpicImageIsOK (nCTg, __func__) != NPIC_SUCCESS) return; if (NpicImageIsSquare (nCTg, __func__) != NPIC_SUCCESS) { nCTg->gen.ok = NPIC_ERROR; return; } NpicFillWhole_i (nCTg, -1); switch (nCTg->type) { case NPIC_IMAGE_2L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_2l *CTg = NpicCastImage(nCTg); int y, x, L = CTg->xmax; CTg->pix[0][0] = 0; for (x = 1; x < L; x++) for (y = 0; y <= x; y++) { CTg->pix[y][x] = x*x + y*y; } } break; case NPIC_IMAGE_3L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_3l *CTg = NpicCastImage(nCTg); int z, y, x, L = CTg->xmax; CTg->pix[0][0][0] = 0; for (x = 1; x < L; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) { CTg->pix[z][y][x] = x*x + y*y + z*z; } } break; case NPIC_IMAGE_4L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_4l *CTg = NpicCastImage(nCTg); int t, z, y, x, L = CTg->xmax; CTg->pix[0][0][0][0] = 0; for (x = 1; x < L; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) { CTg->pix[t][z][y][x] = x*x + y*y + z*z + t*t; } } break; case NPIC_IMAGE_5L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_5l *CTg = NpicCastImage(nCTg); int s, t, z, y, x, L = CTg->xmax; CTg->pix[0][0][0][0][0] = 0; for (x = 1; x < L; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) for (s = 0; s <= t; s++) { CTg->pix[s][t][z][y][x] = x*x + y*y + z*z + t*t + s*s; } } break; case NPIC_IMAGE_6L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_6l *CTg = NpicCastImage(nCTg); int r, s, t, z, y, x, L = CTg->xmax; CTg->pix[0][0][0][0][0][0] = 0; for (x = 1; x < L; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) for (s = 0; s <= t; s++) for (r = 0; r <= s; r++) { CTg->pix[r][s][t][z][y][x] = x*x + y*y + z*z + t*t + s*s + r*r; } } break; default : nCTg->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, ""); } } /* * Fast Distance Transform in G(Z^n). * * Input: nCTg the distance cone to the origin, nDTg init to 0 * R the radius of the ball in CTg. * Output: nDTg the distance image to the background of the ball. * Assume that nCTg and nDTg are L^n images. * * Algorithm derived from Hirata's SEDT [2], with improvements [1] (mark * unpropagated dist. with -1 in cones, stop loops ASAP., store intersections). * * Do nothing if nCTg or nDTg is not ok. set not ok for nDTg on error. Verbose. * Return computation time in second.microsecond */ double NpicMlutDTg_Hirata (Npic_image *nCTg, Npic_image *nDTg, int R) { double time1, time2; if (NpicImageIsOK_DS1 (nDTg, nCTg, __func__) != NPIC_SUCCESS) return 0; if (NpicImageIsSquare (nCTg, __func__) != NPIC_SUCCESS) { nDTg->gen.ok = NPIC_ERROR; return 0; } if (NpicSameImage (nDTg, nCTg, NPIC_TYPE | NPIC_SIZE) != NPIC_SUCCESS) { nDTg->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ""); return 0; } time1 = NpicGetTime(); /* * In all the loops we must have: 0 <= t <= z <= y <= x <= xM < L */ switch (nCTg->type) { case NPIC_IMAGE_2L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_2l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, y, x, propag; Npic_l k, Cv[L]; /* compute bound xM[ and verify that xM < L */ for (xM = 0; xM < L; xM++) if (CTg->pix[0][xM] > R) break; if (xM == L) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } /* First scan */ for (x = 0; x <= xM; x++) { k = 0; propag = 0; for (y = x; y >= 0; y--) { /* outside the ball : background */ if (CTg->pix[y][x] > R) propag = 1; /* inside the ball, mark to dist. k*k from bg */ else if (propag) { k++; DTg->pix[y][x] = k*k; } /* inside the ball, no distance propagated */ else DTg->pix[y][x] = -1; } } /* Scan in x */ for (y = 0; y <= xM; y++) { for (x = y; x <= xM; x++) Cv[x] = DTg->pix[y][x]; NpicMlutHirata_scan (Cv, y, xM); for (x = y; x <= xM; x++) DTg->pix[y][x] = Cv[x]; } } break; case NPIC_IMAGE_3L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_3l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, z, y, x, propag; Npic_l k, Cv[L]; /* compute bound xM[ and verify that xM < L */ for (xM = 0; xM < L; xM++) if (CTg->pix[0][0][xM] > R) break; if (xM == L) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } /* First scan */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) { k = 0; propag = 0; for (z = y; z >= 0; z--) { /* outside the ball : background */ if (CTg->pix[z][y][x] > R) propag = 1; /* inside the ball, mark to dist. k*k from bg */ else if (propag) { k++; DTg->pix[z][y][x] = k*k; } /* inside the ball, no distance propagated */ else DTg->pix[z][y][x] = -1; } } /* Scan in y */ for (x = 0; x <= xM; x++) for (z = 0; z <= x; z++) { for (y = z; y <= x; y++) Cv[y] = DTg->pix[z][y][x]; NpicMlutHirata_scan (Cv, z, x); for (y = z; y <= x; y++) DTg->pix[z][y][x] = Cv[y]; } /* Scan in x */ for (y = 0; y <= xM; y++) for (z = 0; z <= y; z++) { for (x = y; x <= xM; x++) Cv[x] = DTg->pix[z][y][x]; NpicMlutHirata_scan (Cv, y, xM); for (x = y; x <= xM; x++) DTg->pix[z][y][x] = Cv[x]; } } break; case NPIC_IMAGE_4L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_4l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, t, z, y, x, propag; Npic_l k, Cv[L]; /* compute bound xM[ and verify that xM < L */ for (xM = 0; xM < L; xM++) if (CTg->pix[0][0][0][xM] > R) break; if (xM == L) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } /* First scan */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) { k = 0; propag = 0; for (t = z; t >= 0; t--) { /* outside the ball : background */ if (CTg->pix[t][z][y][x] > R) propag = 1; /* inside the ball, mark to dist. k*k from bg */ else if (propag) { k++; DTg->pix[t][z][y][x] = k*k; } /* inside the ball, no distance propagated */ else DTg->pix[t][z][y][x] = -1; } } /* Scan in z */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (t = 0; t <= y; t++) { for (z = t; z <= y; z++) Cv[z] = DTg->pix[t][z][y][x]; NpicMlutHirata_scan (Cv, t, y); for (z = t; z <= y; z++) DTg->pix[t][z][y][x] = Cv[z]; } /* Scan in y */ for (x = 0; x <= xM; x++) for (z = 0; z <= x; z++) for (t = 0; t <= z; t++) { for (y = z; y <= x; y++) Cv[y] = DTg->pix[t][z][y][x]; NpicMlutHirata_scan (Cv, z, x); for (y = z; y <= x; y++) DTg->pix[t][z][y][x] = Cv[y]; } /* Scan in x */ for (y = 0; y <= xM; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) { for (x = y; x <= xM; x++) Cv[x] = DTg->pix[t][z][y][x]; NpicMlutHirata_scan (Cv, y, xM); for (x = y; x <= xM; x++) DTg->pix[t][z][y][x] = Cv[x]; } } break; case NPIC_IMAGE_5L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_5l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, s, t, z, y, x, propag; Npic_l k, Cv[L]; /* compute bound xM[ and verify that xM < L */ for (xM = 0; xM < L; xM++) if (CTg->pix[0][0][0][0][xM] > R) break; if (xM == L) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } /* First scan */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) { k = 0; propag = 0; for (s = t; s >= 0; s--) { /* outside the ball : background */ if (CTg->pix[s][t][z][y][x] > R) propag = 1; /* inside the ball, mark to dist. k*k from bg */ else if (propag) { k++; DTg->pix[s][t][z][y][x] = k*k; } /* inside the ball, no distance propagated */ else DTg->pix[s][t][z][y][x] = -1; } } /* Scan in t */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (s = 0; s <= z; s++) { for (t = s; t <= z; t++) Cv[t] = DTg->pix[s][t][z][y][x]; NpicMlutHirata_scan (Cv, s, z); for (t = s; t <= z; t++) DTg->pix[s][t][z][y][x] = Cv[t]; } /* Scan in z */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (t = 0; t <= y; t++) for (s = 0; s <= t; s++) { for (z = t; z <= y; z++) Cv[z] = DTg->pix[s][t][z][y][x]; NpicMlutHirata_scan (Cv, t, y); for (z = t; z <= y; z++) DTg->pix[s][t][z][y][x] = Cv[z]; } /* Scan in y */ for (x = 0; x <= xM; x++) for (z = 0; z <= x; z++) for (t = 0; t <= z; t++) for (s = 0; s <= t; s++) { for (y = z; y <= x; y++) Cv[y] = DTg->pix[s][t][z][y][x]; NpicMlutHirata_scan (Cv, z, x); for (y = z; y <= x; y++) DTg->pix[s][t][z][y][x] = Cv[y]; } /* Scan in x */ for (y = 0; y <= xM; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) for (s = 0; s <= t; s++) { for (x = y; x <= xM; x++) Cv[x] = DTg->pix[s][t][z][y][x]; NpicMlutHirata_scan (Cv, y, xM); for (x = y; x <= xM; x++) DTg->pix[s][t][z][y][x] = Cv[x]; } } break; case NPIC_IMAGE_6L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_6l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, r, s, t, z, y, x, propag; Npic_l k, Cv[L]; /* compute bound xM[ and verify that xM < L */ for (xM = 0; xM < L; xM++) if (CTg->pix[0][0][0][0][0][xM] > R) break; if (xM == L) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } /* First scan */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) for (s = 0; s <= t; s++) { k = 0; propag = 0; for (r = s; r >= 0; r--) { /* outside the ball : background */ if (CTg->pix[r][s][t][z][y][x] > R) propag = 1; /* inside the ball, mark to dist. k*k from bg */ else if (propag) { k++; DTg->pix[r][s][t][z][y][x] = k*k; } /* inside the ball, no distance propagated */ else DTg->pix[r][s][t][z][y][x] = -1; } } /* Scan in s */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) for (r = 0; r <= t; r++) { for (s = r; s <= t; s++) Cv[s] = DTg->pix[r][s][t][z][y][x]; NpicMlutHirata_scan (Cv, r, t); for (s = r; s <= t; s++) DTg->pix[r][s][t][z][y][x] = Cv[s]; } /* Scan in t */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (z = 0; z <= y; z++) for (s = 0; s <= z; s++) for (r = 0; r <= s; r++) { for (t = s; t <= z; t++) Cv[t] = DTg->pix[r][s][t][z][y][x]; NpicMlutHirata_scan (Cv, s, z); for (t = s; t <= z; t++) DTg->pix[r][s][t][z][y][x] = Cv[t]; } /* Scan in z */ for (x = 0; x <= xM; x++) for (y = 0; y <= x; y++) for (t = 0; t <= y; t++) for (s = 0; s <= t; s++) for (r = 0; r <= s; r++) { for (z = t; z <= y; z++) Cv[z] = DTg->pix[r][s][t][z][y][x]; NpicMlutHirata_scan (Cv, t, y); for (z = t; z <= y; z++) DTg->pix[r][s][t][z][y][x] = Cv[z]; } /* Scan in y */ for (x = 0; x <= xM; x++) for (z = 0; z <= x; z++) for (t = 0; t <= z; t++) for (s = 0; s <= t; s++) for (r = 0; r <= s; r++) { for (y = z; y <= x; y++) Cv[y] = DTg->pix[r][s][t][z][y][x]; NpicMlutHirata_scan (Cv, z, x); for (y = z; y <= x; y++) DTg->pix[r][s][t][z][y][x] = Cv[y]; } /* Scan in x */ for (y = 0; y <= xM; y++) for (z = 0; z <= y; z++) for (t = 0; t <= z; t++) for (s = 0; s <= t; s++) for (r = 0; r <= s; r++) { for (x = y; x <= xM; x++) Cv[x] = DTg->pix[r][s][t][z][y][x]; NpicMlutHirata_scan (Cv, y, xM); for (x = y; x <= xM; x++) DTg->pix[r][s][t][z][y][x] = Cv[x]; } } break; default : nDTg->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, ""); } time2 = NpicGetTime(); return time2 - time1; } /* * Slow Distance Transform in G(Z^n), used to check NpicMlutDTg_Hirata. * * Input: nCTg the distance cone to the origin, * R the radius of the ball in CTg. * Output: nDTg the distance image to the background of the ball. * Assume that nCTg and nDTg are L^n images. * * Do nothing if nCTg or nDTg is not ok. set not ok for nDTg on error. Verbose. * Return computation time in second.microsecond */ double NpicMlutDTg_quadra (Npic_image *nCTg, Npic_image *nDTg, int R) { double time1, time2; if (NpicImageIsOK_DS1 (nDTg, nCTg, __func__) != NPIC_SUCCESS) return 0; if (NpicImageIsSquare (nCTg, __func__) != NPIC_SUCCESS) { nDTg->gen.ok = NPIC_ERROR; return 0; } if (NpicSameImage (nDTg, nCTg, NPIC_TYPE | NPIC_SIZE) != NPIC_SUCCESS) { nDTg->gen.ok = NpicError (__func__, NPIC_ERR_INCOMPAT, ""); return 0; } time1 = NpicGetTime(); /* * In all the loops we must have: 0 <= r <= s <= t <= z <= y <= x <= xM < L */ switch (nCTg->type) { case NPIC_IMAGE_2L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_2l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, y1, x1, y2, x2; Npic_l j, k; /* compute bound xM[ and verify that xM < L-1 */ for (xM = 0; xM < L-1; xM++) if (CTg->pix[0][xM] > R) break; if (xM == L-1) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } for (x1 = 0; x1 <= xM; x1++) for (y1 = 0; y1 <= x1; y1++) { if (CTg->pix[y1][x1] > R) continue; k = L*L*2; /* Scan up to xM+1 to find pixels > R */ for (x2 = 0; x2 <= xM+1; x2++) for (y2 = 0; y2 <= x2; y2++) if (CTg->pix[y2][x2] > R) { j = (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1); if (j < k) k = j; } DTg->pix[y1][x1] = k; } } break; case NPIC_IMAGE_3L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_3l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, z1, y1, x1, z2, y2, x2; Npic_l j, k; /* compute bound xM[ and verify that xM < L-1 */ for (xM = 0; xM < L-1; xM++) if (CTg->pix[0][0][xM] > R) break; if (xM == L-1) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } for (x1 = 0; x1 <= xM; x1++) for (y1 = 0; y1 <= x1; y1++) for (z1 = 0; z1 <= y1; z1++) { if (CTg->pix[z1][y1][x1] > R) continue; k = L*L*3; /* Scan up to xM+1 to find pixels > R */ for (x2 = 0; x2 <= xM+1; x2++) for (y2 = 0; y2 <= x2; y2++) for (z2 = 0; z2 <= y2; z2++) if (CTg->pix[z2][y2][x2] > R) { j = (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1); if (j < k) k = j; } DTg->pix[z1][y1][x1] = k; } } break; case NPIC_IMAGE_4L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_4l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, t1, z1, y1, x1, t2, z2, y2, x2; Npic_l j, k; /* compute bound xM[ and verify that xM < L-1 */ for (xM = 0; xM < L-1; xM++) if (CTg->pix[0][0][0][xM] > R) break; if (xM == L-1) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } for (x1 = 0; x1 <= xM; x1++) for (y1 = 0; y1 <= x1; y1++) for (z1 = 0; z1 <= y1; z1++) for (t1 = 0; t1 <= z1; t1++) { if (CTg->pix[t1][z1][y1][x1] > R) continue; k = L*L*4; /* Scan up to xM+1 to find pixels > R */ for (x2 = 0; x2 <= xM+1; x2++) for (y2 = 0; y2 <= x2; y2++) for (z2 = 0; z2 <= y2; z2++) for (t2 = 0; t2 <= z2; t2++) if (CTg->pix[t2][z2][y2][x2] > R) { j = (t2-t1)*(t2-t1) + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1); if (j < k) k = j; } DTg->pix[t1][z1][y1][x1] = k; } } break; case NPIC_IMAGE_5L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_5l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, s1, t1, z1, y1, x1, s2, t2, z2, y2, x2; Npic_l j, k; /* compute bound xM[ and verify that xM < L-1 */ for (xM = 0; xM < L-1; xM++) if (CTg->pix[0][0][0][0][xM] > R) break; if (xM == L-1) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } for (x1 = 0; x1 <= xM; x1++) for (y1 = 0; y1 <= x1; y1++) for (z1 = 0; z1 <= y1; z1++) for (t1 = 0; t1 <= z1; t1++) for (s1 = 0; s1 <= t1; s1++) { if (CTg->pix[s1][t1][z1][y1][x1] > R) continue; k = L*L*5; /* Scan up to xM+1 to find pixels > R */ for (x2 = 0; x2 <= xM+1; x2++) for (y2 = 0; y2 <= x2; y2++) for (z2 = 0; z2 <= y2; z2++) for (t2 = 0; t2 <= z2; t2++) for (s2 = 0; s2 <= t2; s2++) if (CTg->pix[s2][t2][z2][y2][x2] > R) { j = (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1); if (j < k) k = j; } DTg->pix[s1][t1][z1][y1][x1] = k; } } break; case NPIC_IMAGE_6L : { /* - - - - - - - - - - - - - - - - - - - - - */ Npic_image_6l *CTg = NpicCastImage(nCTg), *DTg = NpicCastImage(nDTg); int xM, L = CTg->xmax, r1, s1, t1, z1, y1, x1, r2, s2, t2, z2, y2, x2; Npic_l j, k; /* compute bound xM[ and verify that xM < L-1 */ for (xM = 0; xM < L-1; xM++) if (CTg->pix[0][0][0][0][0][xM] > R) break; if (xM == L-1) { DTg->ok = NpicError (__func__, NPIC_ERROR, ": R is too large"); return 0; } for (x1 = 0; x1 <= xM; x1++) for (y1 = 0; y1 <= x1; y1++) for (z1 = 0; z1 <= y1; z1++) for (t1 = 0; t1 <= z1; t1++) for (s1 = 0; s1 <= t1; s1++) for (r1 = 0; r1 <= s1; r1++) { if (CTg->pix[r1][s1][t1][z1][y1][x1] > R) continue; k = L*L*6; /* Scan up to xM+1 to find pixels > R */ for (x2 = 0; x2 <= xM+1; x2++) for (y2 = 0; y2 <= x2; y2++) for (z2 = 0; z2 <= y2; z2++) for (t2 = 0; t2 <= z2; t2++) for (s2 = 0; s2 <= t2; s2++) for (r2 = 0; r2 <= s2; r2++) if (CTg->pix[r2][s2][t2][z2][y2][x2] > R) { j = (r2-r1)*(r2-r1) + (s2-s1)*(s2-s1) + (t2-t1)*(t2-t1) + (z2-z1)*(z2-z1) + (y2-y1)*(y2-y1) + (x2-x1)*(x2-x1); if (j < k) k = j; } DTg->pix[r1][s1][t1][z1][y1][x1] = k; } } break; default : nDTg->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 ---------------------*/ /* * Scan column Cv[a..b] and compute dE values ; -1 means no distance. * Used in NpicMlutDTg_Hirata. * * Differences with Npic_Hirata_scan_inf: * - scan [a..b] in place of [0..cmax[ * - stop at second consecutive 0 */ void NpicMlutHirata_scan (Npic_l *Cv, int a, int b) { Npic_l k, Sv[b+1]; /* values */ int j, Si[b+1], Sn; /* index, number */ double Sr[b+1]; /* intersections */ Sn = 0; /* empty stacks */ /* Compute stacks indices Si[Sn] and values Sv[Sn] */ for (j = a; j <= b; j++) { k = Cv[j]; if (k < 0) continue; /* Non propagated value */ /* To speedup algorithm, stop at the second consecutive 0 */ if (k == 0 && j > a && Cv[j-1] == 0) break; while (Sn >= 2 && Npic_Hirata_inter (Si[Sn-1], Sv[Sn-1], j, k) < Sr[Sn-1]) Sn--; /* pop */ /* push */ Si[Sn] = j; Sv[Sn] = k; if (Sn >= 1) Sr[Sn] = Npic_Hirata_inter (Si[Sn-1], Sv[Sn-1], Si[Sn], Sv[Sn]); Sn++; } if (Sn == 0) return; /* Empty stacks, nothing to change in Cv[] */ /* Compute new DTg values using stacks */ for (j = b; j >= a; j--) { if (Cv[j] == 0) continue; while (Sn >= 2 && j < Sr[Sn-1]) Sn--; /* pop */ Cv[j] = Sv[Sn-1] + (j-Si[Sn-1])*(j-Si[Sn-1]); } } /*----------------------------------------------------------------------------*/