/* * 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_sedt_hirata.ct" */ /* * dist_sedt_hirata.c - 09/12/2008 * * Squared Euclidean Distance Transforms. * * The Hirata's [2] and Saito's et al [1] algorithms belongs to the family of * separable transforms: they operate independently in each dimension. * [1] was the first n-dimensional, exact and efficient SEDT, but still a * non-linear algorithm. [2] is the linearized version of [1], and [3] is an * extension of [2]. * For SEDT and SEDT_inf, [1] is a little faster than [2] for small images, * while [2] becomes much faster than [1] for large images. * For SEDT_rev, [2] is generally much faster than [1] for all sizes. * The quadratic algorithms are very slow, but given for checks. * * References: * * [1] T. Saito et J.I. Toriwaki, 1994. "New algorithms for Euclidean distance * transformation of an n-dimensional digitized picture with applications". * Pattern Recognition, 27(11):1551-1565. * * [2] T. Hirata, 1996. "A unified linear-time algorithm for computing * distance maps". Information Processing Letters, 58:129-133. * * [3] D. Coeurjolly, 2003. "d-Dimensional Reverse Euclidean Distance * Transformation and Euclidean Medial Axis Extraction in Optimal Time". * In 11th DGCI, LNCS vol. 2886, Springer-Verlag, p. 327–337, Naples, Italy. */ #include /*--------------------- P U B L I C - I N T E R F A C E ----------------------*/ /* * Compute the Squared Euclidean Distance Transform (SEDT) on image np, * with the optimal in time linear algorithm from Hirata [2], and some fixes. * * Do nothing if np is not ok. Set ok on error. Verbose. * Return computation time in second.microsecond */ double NpicSEDT_Hirata (Npic_image *np) { double time1, time2; int L; if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return 0; /* Change the border width of np to 1 if it was smaller */ NpicSetBorderWidthMin (np, 1); /* 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; } L = np->gen.ymax; if (np->gen.dim >= 3 && np->gen.zmax > L) L = np->gen.zmax; if (np->gen.dim >= 4 && np->gen.tmax > L) L = np->gen.tmax; if (np->gen.dim >= 5 && np->gen.smax > L) L = np->gen.smax; if (np->gen.dim >= 6 && np->gen.rmax > L) L = np->gen.rmax; L += 2; /* for border pixels */ time1 = NpicGetTime(); switch (np->type) { case NPIC_IMAGE_2L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_2l *p = NpicCastImage(np); Npic_l k, j, Cw[L], *Cv=Cw+1; /* values */ int y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(y, x, k, j, Cw, Cv) { Cv=Cw+1; /* again because Cv is private */ # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (y = 0; y < p->ymax; y++) { k = 0; for (x = 0; x <= p->xmax; x++) /* up to border, 0 filled */ { if (p->pix[y][x] != 0) { k++; p->pix[y][x] = k*k; } /* forward */ else if (k != 0) { for (j = k/2; j >= 1; j--) p->pix[y][x-j] = j*j; /* backward */ k = 0; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (x = 0; x < p->xmax; x++) { for (y = -1; y <= p->ymax; y++) /* up to border, 0 filled */ Cv[y] = p->pix[y][x]; Npic_Hirata_scan (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[y][x] = Cv[y]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_3L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_3l *p = NpicCastImage(np); Npic_l k, j, Cw[L], *Cv=Cw+1; /* values */ int z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(z, y, x, k, j, Cw, Cv) { Cv=Cw+1; /* again because Cv is private */ # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) { k = 0; for (x = 0; x <= p->xmax; x++) /* up to border, 0 filled */ { if (p->pix[z][y][x] != 0) { k++; p->pix[z][y][x] = k*k; } /* forward */ else if (k != 0) { for (j = k/2; j >= 1; j--) p->pix[z][y][x-j] = j*j; /* backward */ k = 0; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = -1; y <= p->ymax; y++) /* up to border, 0 filled */ Cv[y] = p->pix[z][y][x]; Npic_Hirata_scan (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = -1; z <= p->zmax; z++) /* up to border, 0 filled */ Cv[z] = p->pix[z][y][x]; Npic_Hirata_scan (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[z][y][x] = Cv[z]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_4L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_4l *p = NpicCastImage(np); Npic_l k, j, Cw[L], *Cv=Cw+1; /* values */ int t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(t, z, y, x, k, j, Cw, Cv) { Cv=Cw+1; /* again because Cv is private */ # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) { k = 0; for (x = 0; x <= p->xmax; x++) /* up to border, 0 filled */ { if (p->pix[t][z][y][x] != 0) { k++; p->pix[t][z][y][x] = k*k; } /* forward */ else if (k != 0) { for (j = k/2; j >= 1; j--) p->pix[t][z][y][x-j] = j*j; /* backward */ k = 0; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = -1; y <= p->ymax; y++) /* up to border, 0 filled */ Cv[y] = p->pix[t][z][y][x]; Npic_Hirata_scan (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = -1; z <= p->zmax; z++) /* up to border, 0 filled */ Cv[z] = p->pix[t][z][y][x]; Npic_Hirata_scan (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = -1; t <= p->tmax; t++) /* up to border, 0 filled */ Cv[t] = p->pix[t][z][y][x]; Npic_Hirata_scan (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[t][z][y][x] = Cv[t]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_5L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_5l *p = NpicCastImage(np); Npic_l k, j, Cw[L], *Cv=Cw+1; /* values */ int s, t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(s, t, z, y, x, k, j, Cw, Cv) { Cv=Cw+1; /* again because Cv is private */ # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { k = 0; for (x = 0; x <= p->xmax; x++) /* up to border, 0 filled */ { if (p->pix[s][t][z][y][x] != 0) { k++; p->pix[s][t][z][y][x] = k*k; } /* forward */ else if (k != 0) { for (j = k/2; j >= 1; j--) p->pix[s][t][z][y][x-j] = j*j; /* backward */ k = 0; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = -1; y <= p->ymax; y++) /* up to border, 0 filled */ Cv[y] = p->pix[s][t][z][y][x]; Npic_Hirata_scan (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[s][t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = -1; z <= p->zmax; z++) /* up to border, 0 filled */ Cv[z] = p->pix[s][t][z][y][x]; Npic_Hirata_scan (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[s][t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = -1; t <= p->tmax; t++) /* up to border, 0 filled */ Cv[t] = p->pix[s][t][z][y][x]; Npic_Hirata_scan (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[s][t][z][y][x] = Cv[t]; } /* * Scan in s */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { for (s = -1; s <= p->smax; s++) /* up to border, 0 filled */ Cv[s] = p->pix[s][t][z][y][x]; Npic_Hirata_scan (Cv, p->smax); for (s = 0; s < p->smax; s++) p->pix[s][t][z][y][x] = Cv[s]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_6L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_6l *p = NpicCastImage(np); Npic_l k, j, Cw[L], *Cv=Cw+1; /* values */ int r, s, t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(r, s, t, z, y, x, k, j, Cw, Cv) { Cv=Cw+1; /* again because Cv is private */ # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { k = 0; for (x = 0; x <= p->xmax; x++) /* up to border, 0 filled */ { if (p->pix[r][s][t][z][y][x] != 0) { k++; p->pix[r][s][t][z][y][x] = k*k; } /* forward */ else if (k != 0) { for (j = k/2; j >= 1; j--) p->pix[r][s][t][z][y][x-j] = j*j; /* backward */ k = 0; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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 (x = 0; x < p->xmax; x++) { for (y = -1; y <= p->ymax; y++) /* up to border, 0 filled */ Cv[y] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[r][s][t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = -1; z <= p->zmax; z++) /* up to border, 0 filled */ Cv[z] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[r][s][t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = -1; t <= p->tmax; t++) /* up to border, 0 filled */ Cv[t] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[r][s][t][z][y][x] = Cv[t]; } /* * Scan in s */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) 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++) { for (s = -1; s <= p->smax; s++) /* up to border, 0 filled */ Cv[s] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan (Cv, p->smax); for (s = 0; s < p->smax; s++) p->pix[r][s][t][z][y][x] = Cv[s]; } /* * Scan in r */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { for (r = -1; r <= p->rmax; r++) /* up to border, 0 filled */ Cv[r] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan (Cv, p->rmax); for (r = 0; r < p->rmax; r++) p->pix[r][s][t][z][y][x] = Cv[r]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } default : np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, ""); } time2 = NpicGetTime(); return time2-time1; } /*****************************************************************************/ /* * Compute the Squared Euclidean Distance Transform (SEDT_inf) on image np * with "infinite border", i.e border pixels propagating no distance value. * Derived from the optimal in time linear algorithm from Hirata [2]. * * Do nothing if np is not ok. Set ok on error. Verbose. * Return computation time in second.microsecond */ double NpicSEDT_Hirata_inf (Npic_image *np) { double time1, time2; int L; if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return 0; /* Change the border width of np to 1 if it was smaller */ NpicSetBorderWidthMin (np, 1); /* Set border pixels to -1, which means "don't propagate distance" */ 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; } L = np->gen.ymax; if (np->gen.dim >= 3 && np->gen.zmax > L) L = np->gen.zmax; if (np->gen.dim >= 4 && np->gen.tmax > L) L = np->gen.tmax; if (np->gen.dim >= 5 && np->gen.smax > L) L = np->gen.smax; if (np->gen.dim >= 6 && np->gen.rmax > L) L = np->gen.rmax; time1 = NpicGetTime(); switch (np->type) { case NPIC_IMAGE_2L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_2l *p = NpicCastImage(np); Npic_l k, k2, j, Cv[L]; /* values */ int y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(y, x, k, k2, j, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (y = 0; y < p->ymax; y++) { /* Forward */ k = -1; for (x = 0; x < p->xmax; x++) { if (p->pix[y][x] == 0) k = 0; else if (k == -1) p->pix[y][x] = -1; /* no distance */ else { k++; p->pix[y][x] = k*k; } } /* Backward */ k = -1; for (x = p->xmax-1; x >= 0; x--) { if (p->pix[y][x] == 0) k = 0; else if (k != -1) { k++; k2 = k*k; j = p->pix[y][x]; if (k2 < j || j == -1) p->pix[y][x] = k2; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[y][x]; Npic_Hirata_scan_inf (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[y][x] = Cv[y]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_3L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_3l *p = NpicCastImage(np); Npic_l k, k2, j, Cv[L]; /* values */ int z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(z, y, x, k, k2, j, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) { /* Forward */ k = -1; for (x = 0; x < p->xmax; x++) { if (p->pix[z][y][x] == 0) k = 0; else if (k == -1) p->pix[z][y][x] = -1; /* no distance */ else { k++; p->pix[z][y][x] = k*k; } } /* Backward */ k = -1; for (x = p->xmax-1; x >= 0; x--) { if (p->pix[z][y][x] == 0) k = 0; else if (k != -1) { k++; k2 = k*k; j = p->pix[z][y][x]; if (k2 < j || j == -1) p->pix[z][y][x] = k2; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[z][y][x]; Npic_Hirata_scan_inf (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[z][y][x]; Npic_Hirata_scan_inf (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[z][y][x] = Cv[z]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_4L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_4l *p = NpicCastImage(np); Npic_l k, k2, j, Cv[L]; /* values */ int t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(t, z, y, x, k, k2, j, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) { /* Forward */ k = -1; for (x = 0; x < p->xmax; x++) { if (p->pix[t][z][y][x] == 0) k = 0; else if (k == -1) p->pix[t][z][y][x] = -1; /* no distance */ else { k++; p->pix[t][z][y][x] = k*k; } } /* Backward */ k = -1; for (x = p->xmax-1; x >= 0; x--) { if (p->pix[t][z][y][x] == 0) k = 0; else if (k != -1) { k++; k2 = k*k; j = p->pix[t][z][y][x]; if (k2 < j || j == -1) p->pix[t][z][y][x] = k2; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = 0; t < p->tmax; t++) Cv[t] = p->pix[t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[t][z][y][x] = Cv[t]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_5L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_5l *p = NpicCastImage(np); Npic_l k, k2, j, Cv[L]; /* values */ int s, t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(s, t, z, y, x, k, k2, j, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { /* Forward */ k = -1; for (x = 0; x < p->xmax; x++) { if (p->pix[s][t][z][y][x] == 0) k = 0; else if (k == -1) p->pix[s][t][z][y][x] = -1; /* no distance */ else { k++; p->pix[s][t][z][y][x] = k*k; } } /* Backward */ k = -1; for (x = p->xmax-1; x >= 0; x--) { if (p->pix[s][t][z][y][x] == 0) k = 0; else if (k != -1) { k++; k2 = k*k; j = p->pix[s][t][z][y][x]; if (k2 < j || j == -1) p->pix[s][t][z][y][x] = k2; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[s][t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[s][t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = 0; t < p->tmax; t++) Cv[t] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[s][t][z][y][x] = Cv[t]; } /* * Scan in s */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { for (s = 0; s < p->smax; s++) Cv[s] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->smax); for (s = 0; s < p->smax; s++) p->pix[s][t][z][y][x] = Cv[s]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_6L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_6l *p = NpicCastImage(np); Npic_l k, k2, j, Cv[L]; /* values */ int r, s, t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(r, s, t, z, y, x, k, k2, j, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { /* Forward */ k = -1; for (x = 0; x < p->xmax; x++) { if (p->pix[r][s][t][z][y][x] == 0) k = 0; else if (k == -1) p->pix[r][s][t][z][y][x] = -1; /* no distance */ else { k++; p->pix[r][s][t][z][y][x] = k*k; } } /* Backward */ k = -1; for (x = p->xmax-1; x >= 0; x--) { if (p->pix[r][s][t][z][y][x] == 0) k = 0; else if (k != -1) { k++; k2 = k*k; j = p->pix[r][s][t][z][y][x]; if (k2 < j || j == -1) p->pix[r][s][t][z][y][x] = k2; } } } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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 (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[r][s][t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[r][s][t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = 0; t < p->tmax; t++) Cv[t] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[r][s][t][z][y][x] = Cv[t]; } /* * Scan in s */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) 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++) { for (s = 0; s < p->smax; s++) Cv[s] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->smax); for (s = 0; s < p->smax; s++) p->pix[r][s][t][z][y][x] = Cv[s]; } /* * Scan in r */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { for (r = 0; r < p->rmax; r++) Cv[r] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_inf (Cv, p->rmax); for (r = 0; r < p->rmax; r++) p->pix[r][s][t][z][y][x] = Cv[r]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } default : np->gen.ok = NpicError (__func__, NPIC_ERR_UNEX_NPIC, ""); } time2 = NpicGetTime(); return time2-time1; } /*****************************************************************************/ /* * Compute the Squared Euclidean Reverse Distance Transform (SEDT_rev) on image np. * Derived from the optimal in time linear algorithm from Hirata [2] (see also [3]). * * Do nothing if np is not ok. Set ok on error. Verbose. * Return computation time in second.microsecond */ double NpicSEDT_Hirata_rev (Npic_image *np) { double time1, time2; int L; if (NpicImageIsOK (np, __func__) != NPIC_SUCCESS) return 0; L = np->gen.xmax; if (np->gen.ymax > L) L = np->gen.ymax; if (np->gen.dim >= 3 && np->gen.zmax > L) L = np->gen.zmax; if (np->gen.dim >= 4 && np->gen.tmax > L) L = np->gen.tmax; if (np->gen.dim >= 5 && np->gen.smax > L) L = np->gen.smax; if (np->gen.dim >= 6 && np->gen.rmax > L) L = np->gen.rmax; time1 = NpicGetTime(); switch (np->type) { case NPIC_IMAGE_2L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_2l *p = NpicCastImage(np); Npic_l Cv[L]; /* values */ int y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(y, x, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (y = 0; y < p->ymax; y++) { for (x = 0; x < p->xmax; x++) Cv[x] = p->pix[y][x]; Npic_Hirata_scan_rev (Cv, p->xmax); for (x = 0; x < p->xmax; x++) p->pix[y][x] = Cv[x]; } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[y][x]; Npic_Hirata_scan_rev (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[y][x] = Cv[y]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_3L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_3l *p = NpicCastImage(np); Npic_l Cv[L]; /* values */ int z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(z, y, x, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) { for (x = 0; x < p->xmax; x++) Cv[x] = p->pix[z][y][x]; Npic_Hirata_scan_rev (Cv, p->xmax); for (x = 0; x < p->xmax; x++) p->pix[z][y][x] = Cv[x]; } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[z][y][x]; Npic_Hirata_scan_rev (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[z][y][x]; Npic_Hirata_scan_rev (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[z][y][x] = Cv[z]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_4L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_4l *p = NpicCastImage(np); Npic_l Cv[L]; /* values */ int t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(t, z, y, x, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) Cv[x] = p->pix[t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->xmax); for (x = 0; x < p->xmax; x++) p->pix[t][z][y][x] = Cv[x]; } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = 0; t < p->tmax; t++) Cv[t] = p->pix[t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[t][z][y][x] = Cv[t]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_5L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_5l *p = NpicCastImage(np); Npic_l Cv[L]; /* values */ int s, t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(s, t, z, y, x, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) Cv[x] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->xmax); for (x = 0; x < p->xmax; x++) p->pix[s][t][z][y][x] = Cv[x]; } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (z = 0; z < p->zmax; z++) for (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[s][t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[s][t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (s = 0; s < p->smax; s++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = 0; t < p->tmax; t++) Cv[t] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[s][t][z][y][x] = Cv[t]; } /* * Scan in s */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { for (s = 0; s < p->smax; s++) Cv[s] = p->pix[s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->smax); for (s = 0; s < p->smax; s++) p->pix[s][t][z][y][x] = Cv[s]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif break; } case NPIC_IMAGE_6L : { /* - - - - - - - - - - - - - - - - - */ Npic_image_6l *p = NpicCastImage(np); Npic_l Cv[L]; /* values */ int r, s, t, z, y, x; # ifdef NPIC_USE_OMP # pragma omp parallel private(r, s, t, z, y, x, Cv) { # endif /* * First scan in x */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) Cv[x] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->xmax); for (x = 0; x < p->xmax; x++) p->pix[r][s][t][z][y][x] = Cv[x]; } /* * Scan in y */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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 (x = 0; x < p->xmax; x++) { for (y = 0; y < p->ymax; y++) Cv[y] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->ymax); for (y = 0; y < p->ymax; y++) p->pix[r][s][t][z][y][x] = Cv[y]; } /* * Scan in z */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (t = 0; t < p->tmax; t++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (z = 0; z < p->zmax; z++) Cv[z] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->zmax); for (z = 0; z < p->zmax; z++) p->pix[r][s][t][z][y][x] = Cv[z]; } /* * Scan in t */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) for (s = 0; s < p->smax; s++) for (z = 0; z < p->zmax; z++) for (y = 0; y < p->ymax; y++) for (x = 0; x < p->xmax; x++) { for (t = 0; t < p->tmax; t++) Cv[t] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->tmax); for (t = 0; t < p->tmax; t++) p->pix[r][s][t][z][y][x] = Cv[t]; } /* * Scan in s */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif for (r = 0; r < p->rmax; r++) 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++) { for (s = 0; s < p->smax; s++) Cv[s] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->smax); for (s = 0; s < p->smax; s++) p->pix[r][s][t][z][y][x] = Cv[s]; } /* * Scan in r */ # ifdef NPIC_USE_OMP # pragma omp for schedule(dynamic) # endif 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++) { for (r = 0; r < p->rmax; r++) Cv[r] = p->pix[r][s][t][z][y][x]; Npic_Hirata_scan_rev (Cv, p->rmax); for (r = 0; r < p->rmax; r++) p->pix[r][s][t][z][y][x] = Cv[r]; } # ifdef NPIC_USE_OMP } /* omp parallel */ # endif 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 ---------------------*/ /* * Return intersection abscissa of parabolas */ double Npic_Hirata_inter (int u, int gu, int v, int gv) { return (u+v + ((double)gu-gv)/(u-v)) / 2; } /* * Scan column Cv[-1..cmax] and compute dE values. * Used in NpicSEDT_Hirata. */ void Npic_Hirata_scan (Npic_l *Cv, int cmax) { Npic_l k, Sv[cmax+2]; /* values */ int j, Si[cmax+2], Sn; /* index, number */ double Sr[cmax+2]; /* intersections */ Sn = 0; /* Empty stacks */ /* Compute stacks indices Si[Sn] and SEDT values Sv[Sn] */ for (j = -1; j <= cmax; j++) /* up to border, 0 filled */ { k = Cv[j]; 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 SEDT values using stacks */ for (j = cmax-1; j >= 0; 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]); } } /* * Scan column Cv[0..cmax-1] and compute dE values ; -1 means no distance. * Used in NpicSEDT_Hirata_inf. */ void Npic_Hirata_scan_inf (Npic_l *Cv, int cmax) { Npic_l k, Sv[cmax]; /* values */ int j, Si[cmax], Sn; /* index, number */ double Sr[cmax]; /* intersections */ Sn = 0; /* Empty stacks */ /* Compute stacks indices Si[Sn] and SEDT values Sv[Sn] */ for (j = 0; j < cmax; j++) { k = Cv[j]; if (k == -1) continue; /* no distance */ 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 SEDT values using stacks */ for (j = cmax-1; j >= 0; 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]); } } /* * Scan column Cv[0..cmax-1] and compute dE values. * Used in NpicSEDT_Hirata_rev. */ void Npic_Hirata_scan_rev (Npic_l *Cv, int cmax) { Npic_l k, Sv[cmax]; /* values */ int j, Si[cmax], Sn; /* index, number */ double Sr[cmax]; /* intersections */ Sn = 0; /* Empty stacks */ /* Compute stacks indices Si[Sn] and RevSEDT values Sv[Sn] */ for (j = 0; j < cmax; j++) { k = Cv[j]; 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 RevSEDT values using stacks */ for (j = cmax-1; j >= 0; j--) { while (Sn >= 2 && j < Sr[Sn-1]) Sn--; /* pop */ k = Sv[Sn-1] - (j-Si[Sn-1])*(j-Si[Sn-1]); Cv[j] = k > 0 ? k : 0; /* see: k seems to be never negative */ } } /*----------------------------------------------------------------------------*/