diff --git a/bindings/python/py-distance-transform.cpp b/bindings/python/py-distance-transform.cpp index a0d9cbb..dda4a96 100644 --- a/bindings/python/py-distance-transform.cpp +++ b/bindings/python/py-distance-transform.cpp @@ -16,7 +16,7 @@ namespace py = pybind11; template py::array_t DistanceTransformImpl( py::array_t maskarray, - bool computeSquareDistance) { + bool computeSquareDistance, const std::vector& alphas) { // Get input shape typename dope::DopeVector::IndexD masksize, pymasksize; std::copy_n(maskarray.shape(), N, pymasksize.begin()); @@ -31,7 +31,7 @@ py::array_t DistanceTransformImpl( dope::Grid dopefield(masksize); dt::DistanceTransform::distanceTransformL2(dopemask, dopefield, - computeSquareDistance); + computeSquareDistance, alphas); return py::array_t(pymasksize, dopefield.data()); } @@ -39,27 +39,47 @@ py::array_t DistanceTransformImpl( template py::array_t DistanceTransform( py::array_t maskarray, - bool computeSquareDistance = false) { + bool computeSquareDistance = false, + std::optional > + optalphas = py::none()) { // std::cout<<"Got input with dtype="<(maskarray.ndim(), 1.0); + } + switch (maskarray.ndim()) { case 1: { - return DistanceTransformImpl(maskarray, computeSquareDistance); + return DistanceTransformImpl(maskarray, computeSquareDistance, + alphas); } case 2: { - return DistanceTransformImpl(maskarray, computeSquareDistance); + return DistanceTransformImpl(maskarray, computeSquareDistance, + alphas); } case 3: { - return DistanceTransformImpl(maskarray, computeSquareDistance); + return DistanceTransformImpl(maskarray, computeSquareDistance, + alphas); } case 4: { - return DistanceTransformImpl(maskarray, computeSquareDistance); + return DistanceTransformImpl(maskarray, computeSquareDistance, + alphas); } case 5: { - return DistanceTransformImpl(maskarray, computeSquareDistance); + return DistanceTransformImpl(maskarray, computeSquareDistance, + alphas); } case 6: { - return DistanceTransformImpl(maskarray, computeSquareDistance); + return DistanceTransformImpl(maskarray, computeSquareDistance, + alphas); } default: { throw std::out_of_range("Dimension " + std::to_string(maskarray.ndim()) + @@ -82,5 +102,6 @@ PYBIND11_MODULE(py_distance_transform, m) { Compute the distance transform https://github.com/tvercaut/distance_transform )pbdoc", - py::arg("maskarray"), py::arg("computeSquareDistance") = false); + py::arg("maskarray"), py::arg("computeSquareDistance") = false, + py::arg("alphas") = py::none()); } diff --git a/example/distance_transform_example.cpp b/example/distance_transform_example.cpp index e5fdb93..bd99232 100644 --- a/example/distance_transform_example.cpp +++ b/example/distance_transform_example.cpp @@ -92,7 +92,8 @@ int runexample() { std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); - dt::DistanceTransform::distanceTransformL2(f, f, indices, false, 1); + dt::DistanceTransform::distanceTransformL2(f, f, indices, false, + std::vector(2, 1.0), 1); std::cout << std::endl << "2D distance function computed in: " << std::chrono::duration_cast( @@ -137,7 +138,8 @@ int runexample() { } start = std::chrono::steady_clock::now(); - dt::DistanceTransform::distanceTransformL2(fWin, fWin, indicesWin, true, 1); + dt::DistanceTransform::distanceTransformL2(fWin, fWin, indicesWin, true, + std::vector(2, 1.0), 1); std::cout << std::endl << "2D distance function computed on the window in: " << std::chrono::duration_cast( @@ -168,7 +170,8 @@ int runexample() { f2D[i][j] = std::numeric_limits::max(); f2D[0][0] = 0.0f; start = std::chrono::steady_clock::now(); - dt::DistanceTransform::distanceTransformL2(f2D, f2D, false, 1); + dt::DistanceTransform::distanceTransformL2(f2D, f2D, false, + std::vector(2, 1.0), 1); std::cout << std::endl << size[0] << 'x' << size[1] << " distance function computed in: " << std::chrono::duration_cast( @@ -185,7 +188,8 @@ int runexample() { f3D[i][j][k] = std::numeric_limits::max(); f3D[0][0][0] = 0.0f; start = std::chrono::steady_clock::now(); - dt::DistanceTransform::distanceTransformL2(f3D, f3D, false, 1); + dt::DistanceTransform::distanceTransformL2(f3D, f3D, false, + std::vector(3, 1.0), 1); std::cout << std::endl << size3D[0] << 'x' << size3D[1] << 'x' << size3D[2] << " distance function computed in: " @@ -202,7 +206,8 @@ int runexample() { f3D[0][0][0] = 0.0f; start = std::chrono::steady_clock::now(); dt::DistanceTransform::distanceTransformL2( - f3D, f3D, false, std::thread::hardware_concurrency()); + f3D, f3D, false, std::vector(3, 1.0), + std::thread::hardware_concurrency()); std::cout << std::endl << size3D[0] << 'x' << size3D[1] << 'x' << size3D[2] << " distance function (concurrently) computed in: " @@ -220,7 +225,8 @@ int runexample() { f2DBig[i][j] = std::numeric_limits::max(); f2DBig[0][0] = 0.0f; start = std::chrono::steady_clock::now(); - dt::DistanceTransform::distanceTransformL2(f2DBig, f2DBig, false, 1); + dt::DistanceTransform::distanceTransformL2(f2DBig, f2DBig, false, + std::vector(2, 1.0), 1); std::cout << std::endl << size[0] << 'x' << size[1] << " distance function computed in: " << std::chrono::duration_cast( @@ -235,7 +241,8 @@ int runexample() { f2DBig[0][0] = 0.0f; start = std::chrono::steady_clock::now(); dt::DistanceTransform::distanceTransformL2( - f2DBig, f2DBig, false, std::thread::hardware_concurrency()); + f2DBig, f2DBig, false, std::vector(2, 1.0), + std::thread::hardware_concurrency()); std::cout << std::endl << size[0] << 'x' << size[1] << " distance function (concurrently) computed in: " @@ -257,7 +264,8 @@ int runexample() { f6D[i][j][k][l][m][n] = std::numeric_limits::max(); f6D[0][0][0][0][0][0] = 0.0f; start = std::chrono::steady_clock::now(); - dt::DistanceTransform::distanceTransformL2(f6D, f6D, false, 1); + dt::DistanceTransform::distanceTransformL2(f6D, f6D, false, + std::vector(6, 1.0), 1); std::cout << std::endl << size6D[0] << 'x' << size6D[1] << 'x' << size6D[2] << 'x' << size6D[3] << 'x' << size6D[4] << 'x' << size6D[5] diff --git a/example/py-distance-transform-minimal-example.py b/example/py-distance-transform-minimal-example.py index bb22a9d..9c1e70e 100644 --- a/example/py-distance-transform-minimal-example.py +++ b/example/py-distance-transform-minimal-example.py @@ -4,11 +4,13 @@ def distance_transform_example(): - mask = np.ones((51,100), dtype=float)*1e2 - mask[10,23] = 0 - mask[35,84] = 0 + mask = np.ones((51, 100), dtype=float) * 1e3 + mask[10, 23] = 0 + mask[35, 84] = 0 - distance_map = dt.distance_transform(mask,True) + spacings = [2, 1.0] + + distance_map = dt.distance_transform(mask, False, spacings) print("Max dist:", np.max(distance_map[:])) print(distance_map.dtype) diff --git a/include/distance_transform/distance_transform.h b/include/distance_transform/distance_transform.h index 5d0e44a..fc8935e 100644 --- a/include/distance_transform/distance_transform.h +++ b/include/distance_transform/distance_transform.h @@ -14,6 +14,7 @@ #define INCLUDE_DISTANCE_TRANSFORM_DISTANCE_TRANSFORM_H_ #include +#include #include "dope_vector/Grid.h" @@ -32,6 +33,8 @@ class DistanceTransform { * @param D The resulting distance field of f. * @param squared Compute squared distances (L2)^2 - avoiding to * compute square roots - (true) or keep them normal (false - default). + * @param alphas Weighting factor for anisotropic distances (square + * of pixel/voxel spacing) * @param nThreads The number of threads for parallel computation. If * <= 1, the computation will be sequential. * @note Arrays f and D can also be the same. @@ -40,6 +43,7 @@ class DistanceTransform { inline static void distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, const bool squared = false, + std::vector alphas = std::vector(DIM, 1.0), const std::size_t nThreads = std::thread::hardware_concurrency()); /** @@ -49,6 +53,8 @@ class DistanceTransform { * @param D The resulting distance field of f. * @param squared Compute squared distances (L2)^2 - avoiding to * compute square roots - (true) or keep them normal (false - default). + * @param alphas Weighting factor for anisotropic distances (square + * of pixel/voxel spacing) * @param nThreads The number of threads for parallel computation. * Actually NOT used, since it's not easy to run a single row computation in * parallel. @@ -58,6 +64,7 @@ class DistanceTransform { inline static void distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, const bool squared = false, + std::vector alphas = std::vector(1, 1.0), const std::size_t nThreads = std::thread::hardware_concurrency()); /** @@ -70,6 +77,8 @@ class DistanceTransform { * local minimum for each sample. * @param squared Compute squared distances (L2)^2 - avoiding to * compute square roots - (true) or keep them normal (false - default). + * @param alphas Weighting factor for anisotropic distances (square + * of pixel/voxel spacing) * @param nThreads The number of threads for parallel computation. If * <= 1, the computation will be sequential. * @note Arrays f and D can also be the same. I should be first @@ -79,6 +88,7 @@ class DistanceTransform { inline static void distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, dope::DopeVector &I, const bool squared = false, + std::vector alphas = std::vector(DIM, 1.0), const std::size_t nThreads = std::thread::hardware_concurrency()); /** @@ -91,6 +101,8 @@ class DistanceTransform { * local minimum for each sample. * @param squared Compute squared distances (L2)^2 - avoiding to * compute square roots - (true) or keep them normal (false - default). + * @param alphas Weighting factor for anisotropic distances (square + * of pixel/voxel spacing) * @param nThreads The number of threads for parallel computation. * Actually NOT used, since it's not easy to run a single row computation in * parallel. @@ -100,6 +112,7 @@ class DistanceTransform { inline static void distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, dope::DopeVector &I, const bool squared = false, + std::vector alphas = std::vector(1, 1.0), const std::size_t nThreads = std::thread::hardware_concurrency()); /** @@ -126,12 +139,13 @@ class DistanceTransform { * window, in multi-threading). * @param D The resulting distance field of f (a window, in * multi-threading). - * @param d The dimension where to slice. - * @param order The order in which to permute the slices. + * @param alpha A multiplier stretching each parabola (L2 distance) + * vertically. */ template inline static void distanceL2Helper(const dope::DopeVector &f, - dope::DopeVector &D); + dope::DopeVector &D, + const Scalar alpha); /** * @brief The actual distance field computation is done by recursive calls @@ -141,7 +155,8 @@ class DistanceTransform { */ template inline static void distanceL2(const dope::DopeVector &f, - dope::DopeVector &D); + dope::DopeVector &D, + const Scalar alpha); /** * @brief The actual distance field computation as in the "Distance @@ -152,7 +167,8 @@ class DistanceTransform { */ template inline static void distanceL2(const dope::DopeVector &f, - dope::DopeVector &D); + dope::DopeVector &D, + const Scalar alpha); /** * @brief The loop iteration process that can be executed sequentially and @@ -170,7 +186,7 @@ class DistanceTransform { inline static void distanceL2Helper( const dope::DopeVector &f, dope::DopeVector &D, const dope::DopeVector &Ipre, - dope::DopeVector &Ipost); + dope::DopeVector &Ipost, const Scalar alpha); /** * @brief The actual distance field computation is done by recursive calls @@ -184,7 +200,7 @@ class DistanceTransform { inline static void distanceL2( const dope::DopeVector &f, dope::DopeVector &D, const dope::DopeVector &Ipre, - dope::DopeVector &Ipost); + dope::DopeVector &Ipost, const Scalar alpha); /** * @brief The actual distance field computation as in the "Distance @@ -199,7 +215,8 @@ class DistanceTransform { inline static void distanceL2(const dope::DopeVector &f, dope::DopeVector &D, const dope::DopeVector &Ipre, - dope::DopeVector &Ipost); + dope::DopeVector &Ipost, + const Scalar alpha); public: /** diff --git a/include/distance_transform/inlines/distance_transform.hpp b/include/distance_transform/inlines/distance_transform.hpp index c8b8a85..7f7d539 100644 --- a/include/distance_transform/inlines/distance_transform.hpp +++ b/include/distance_transform/inlines/distance_transform.hpp @@ -26,13 +26,20 @@ namespace dt { template inline void DistanceTransform::distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, - const bool squared, const std::size_t nThreads) { - dope::Index fSize, DSize; - fSize = f.allSizes(); - DSize = D.allSizes(); - if (DSize != fSize) - throw std::out_of_range("Matrixes do not have same size."); + const bool squared, std::vector alphas, + const std::size_t nThreads) { + // Check that the output array has the same size as the input array + const dope::Index fSize = f.allSizes(); + const dope::Index DSize = D.allSizes(); + if (DSize != fSize) { + throw std::out_of_range("Matrices do not have same size."); + } + const auto aSize = alphas.size(); + if (aSize != DIM) { + throw std::out_of_range("Spacing vector size is not equal to dimension."); + } + // Allocate working buffers dope::Grid fCopy(fSize); fCopy.import(f); dope::Grid DCopy(DSize); @@ -42,15 +49,23 @@ inline void DistanceTransform::distanceTransformL2( dope::Index order; + // A two-dimensional distance transform can be computed by first computing + // one-dimensional distance transforms along each column of the grid, and then + // computing one-dimensional distance transforms along each row of the result. + // This argument extends to arbitrary dimensions, resulting in the composition + // of transforms along each dimension of the underlying grid. for (dope::SizeType d = static_cast(0); d < DIM; ++d) { - // permute rotate - for (dope::SizeType o = static_cast(0); o < DIM; ++o) + // Rather than changing the direction in which we scan though the array, + // we permute and rotate the array and then operate on a fixed direction + for (dope::SizeType o = static_cast(0); o < DIM; ++o) { order[o] = (d + o) % DIM; + } dope::DopeVector tmpF_rotated = tmpF.permute(order); dope::DopeVector tmpD_rotated = tmpD.permute(order); - dope::Index winStart = dope::Index::Zero(), winSize; - tmpF_rotated.allSizes(winSize); + // Divide the image in various windows for multithreading purposes + dope::Index winStart = dope::Index::Zero(); + dope::Index winSize = tmpF_rotated.allSizes(); dope::SizeType range = winSize[0]; if (nThreads < range) { @@ -74,51 +89,68 @@ inline void DistanceTransform::distanceTransformL2( winSize[0] = tmpF_rotated.sizeAt(0); threads.at(i) = std::thread( static_cast &, - dope::DopeVector &)>( - &distanceL2Helper), - std::cref(tmpWindowsF.at(i)), std::ref(tmpWindowsD.at(i))); + dope::DopeVector &, + const Scalar)>(&distanceL2Helper), + std::cref(tmpWindowsF.at(i)), std::ref(tmpWindowsD.at(i)), + alphas[d]); + } + for (std::size_t i = 0; i < nWindows; ++i) { + threads.at(i).join(); } - for (std::size_t i = 0; i < nWindows; ++i) threads.at(i).join(); } else { - distanceL2Helper(tmpF_rotated, tmpD_rotated); + distanceL2Helper(tmpF_rotated, tmpD_rotated, alphas[d]); } std::swap(tmpD, tmpF); } - if (DIM % 2 == 0) DCopy = std::move(fCopy); + if (DIM % 2 == 0) { + DCopy = std::move(fCopy); + } D.import(DCopy); - if (!squared) element_wiseSquareRoot(D); + if (!squared) { + element_wiseSquareRoot(D); + } } template inline void DistanceTransform::distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, - const bool squared, const std::size_t) { - dope::Index1 fSize, DSize; - fSize = f.allSizes(); - DSize = D.allSizes(); - if (DSize != fSize) + const bool squared, std::vector alphas, const std::size_t) { + const dope::Index1 fSize = f.allSizes(); + const dope::Index1 DSize = D.allSizes(); + if (DSize != fSize) { throw std::out_of_range("Matrixes do not have same size."); + } + const auto aSize = alphas.size(); + if (aSize != 1) { + throw std::out_of_range("Spacing vector size is not equal to dimension."); + } - distanceL2(f, D); + distanceL2(f, D, alphas[0]); - if (!squared) element_wiseSquareRoot(D); + if (!squared) { + element_wiseSquareRoot(D); + } } template inline void DistanceTransform::distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, dope::DopeVector &I, const bool squared, - const std::size_t nThreads) { - dope::Index fSize, DSize, ISize; - fSize = f.allSizes(); - DSize = D.allSizes(); - ISize = I.allSizes(); - if (DSize != fSize || ISize != fSize) + std::vector alphas, const std::size_t nThreads) { + const dope::Index fSize = f.allSizes(); + const dope::Index DSize = D.allSizes(); + const dope::Index ISize = I.allSizes(); + if (DSize != fSize || ISize != fSize) { throw std::out_of_range("Matrixes do not have same size."); + } + const auto aSize = alphas.size(); + if (aSize != DIM) { + throw std::out_of_range("Spacing vector size is not equal to dimension."); + } dope::Grid fCopy(fSize); fCopy.import(f); @@ -135,8 +167,9 @@ inline void DistanceTransform::distanceTransformL2( for (dope::SizeType d = static_cast(0); d < DIM; ++d) { // permute rotate - for (dope::SizeType o = static_cast(0); o < DIM; ++o) + for (dope::SizeType o = static_cast(0); o < DIM; ++o) { order[o] = (d + o) % DIM; + } dope::DopeVector tmpF_rotated = tmpF.permute(order); dope::DopeVector tmpD_rotated = tmpD.permute(order); dope::DopeVector Ipre_rotated = Ipre.permute(order); @@ -175,14 +208,18 @@ inline void DistanceTransform::distanceTransformL2( static_cast &, dope::DopeVector &, const dope::DopeVector &, - dope::DopeVector &)>( - &distanceL2Helper), + dope::DopeVector &, + const Scalar)>(&distanceL2Helper), std::cref(tmpWindowsF.at(i)), std::ref(tmpWindowsD.at(i)), - std::cref(tmpWindowsIPre.at(i)), std::ref(tmpWindowsIPost.at(i))); + std::cref(tmpWindowsIPre.at(i)), std::ref(tmpWindowsIPost.at(i)), + alphas[d]); + } + for (std::size_t i = 0; i < nWindows; ++i) { + threads.at(i).join(); } - for (std::size_t i = 0; i < nWindows; ++i) threads.at(i).join(); } else { - distanceL2Helper(tmpF_rotated, tmpD_rotated, Ipre_rotated, Ipost_rotated); + distanceL2Helper(tmpF_rotated, tmpD_rotated, Ipre_rotated, Ipost_rotated, + alphas[d]); } std::swap(tmpD, tmpF); @@ -197,24 +234,32 @@ inline void DistanceTransform::distanceTransformL2( D.import(DCopy); I.import(ICopyPost); - if (!squared) element_wiseSquareRoot(D); + if (!squared) { + element_wiseSquareRoot(D); + } } template inline void DistanceTransform::distanceTransformL2( const dope::DopeVector &f, dope::DopeVector &D, dope::DopeVector &I, const bool squared, - const std::size_t) { - dope::Index1 fSize, DSize, ISize; - fSize = f.allSizes(); - DSize = D.allSizes(); - ISize = I.allSizes(); - if (DSize != fSize || ISize != fSize) + std::vector alphas, const std::size_t) { + const dope::Index1 fSize = f.allSizes(); + const dope::Index1 DSize = D.allSizes(); + const dope::Index1 ISize = I.allSizes(); + if (DSize != fSize || ISize != fSize) { throw std::out_of_range("Matrixes do not have same size."); + } + const auto aSize = alphas.size(); + if (aSize != 1) { + throw std::out_of_range("Spacing vector size is not equal to dimension."); + } - distanceL2(f, D, I); + distanceL2(f, D, I, Scalar(1.0)); - if (!squared) element_wiseSquareRoot(D); + if (!squared) { + element_wiseSquareRoot(D); + } } template @@ -230,13 +275,16 @@ inline void DistanceTransform::initializeIndices( inline void DistanceTransform::initializeIndices( dope::DopeVector &I) { - for (dope::SizeType q = static_cast(0); q < I.sizeAt(0); ++q) + for (dope::SizeType q = static_cast(0); q < I.sizeAt(0); + ++q) { I[q] = I.accumulatedOffset(q); + } } template inline void DistanceTransform::distanceL2Helper( - const dope::DopeVector &f, dope::DopeVector &D) { + const dope::DopeVector &f, dope::DopeVector &D, + const Scalar alpha) { dope::DopeVector f_dq; dope::DopeVector D_dq; @@ -244,29 +292,32 @@ inline void DistanceTransform::distanceL2Helper( ++q) { f.slice(0, q, f_dq); D.slice(0, q, D_dq); - distanceL2(f_dq, D_dq); + distanceL2(f_dq, D_dq, alpha); } } template inline void DistanceTransform::distanceL2( - const dope::DopeVector &f, dope::DopeVector &D) { + const dope::DopeVector &f, dope::DopeVector &D, + const Scalar alpha) { dope::DopeVector f_q, D_q; // compute distance at lower dimensions for each hyperplane for (dope::SizeType q = static_cast(0); q < f.sizeAt(0); ++q) { f.at(q, f_q); D.at(q, D_q); - distanceL2(f_q, D_q); + distanceL2(f_q, D_q, alpha); } } template inline void DistanceTransform::distanceL2(const dope::DopeVector &f, - dope::DopeVector &D) { + dope::DopeVector &D, + const Scalar alpha) { if (f.sizeAt(0) == static_cast(0) || - f.sizeAt(0) > D.sizeAt(0)) + f.sizeAt(0) > D.sizeAt(0)) { return; + } if (f.sizeAt(0) == static_cast(1)) { D[0] = f[0]; return; @@ -290,7 +341,8 @@ inline void DistanceTransform::distanceL2(const dope::DopeVector &f, --k; // compute horizontal position of intersection between the parabola from q // and the current lowest parabola - s = ((f[q] + q * q) - static_cast(f[v[k]] + v[k] * v[k])) / + s = ((f[q] / alpha + q * q) - + static_cast(f[v[k]] / alpha + v[k] * v[k])) / (static_cast(2 * q) - static_cast(2 * v[k])); } while (s <= z[k]); ++k; @@ -302,9 +354,11 @@ inline void DistanceTransform::distanceL2(const dope::DopeVector &f, k = static_cast(0); for (dope::SizeType q = static_cast(0); q < f.sizeAt(0); ++q) { - while (z[k + 1] < static_cast(q)) ++k; - D[q] = f[v[k]] + - (q - static_cast(v[k])) * (q - static_cast(v[k])); + while (z[k + 1] < static_cast(q)) { + ++k; + } + D[q] = f[v[k]] + alpha * (q - static_cast(v[k])) * + (q - static_cast(v[k])); } } @@ -312,7 +366,7 @@ template inline void DistanceTransform::distanceL2Helper( const dope::DopeVector &f, dope::DopeVector &D, const dope::DopeVector &Ipre, - dope::DopeVector &Ipost) { + dope::DopeVector &Ipost, const Scalar alpha) { dope::DopeVector f_dq; dope::DopeVector D_dq; dope::DopeVector Ipre_dq; @@ -324,7 +378,7 @@ inline void DistanceTransform::distanceL2Helper( D.slice(0, q, D_dq); Ipre.slice(0, q, Ipre_dq); Ipost.slice(0, q, Ipost_dq); - distanceL2(f_dq, D_dq, Ipre_dq, Ipost_dq); + distanceL2(f_dq, D_dq, Ipre_dq, Ipost_dq, alpha); } } @@ -332,7 +386,7 @@ template inline void DistanceTransform::distanceL2( const dope::DopeVector &f, dope::DopeVector &D, const dope::DopeVector &Ipre, - dope::DopeVector &Ipost) { + dope::DopeVector &Ipost, const Scalar alpha) { dope::DopeVector f_q, D_q; dope::DopeVector Ipre_q, Ipost_q; // compute distance at lower dimensions for each hyperplane @@ -342,7 +396,7 @@ inline void DistanceTransform::distanceL2( D.at(q, D_q); Ipre.at(q, Ipre_q); Ipost.at(q, Ipost_q); - distanceL2(f_q, D_q, Ipre_q, Ipost_q); + distanceL2(f_q, D_q, Ipre_q, Ipost_q, alpha); } } @@ -350,10 +404,11 @@ template inline void DistanceTransform::distanceL2( const dope::DopeVector &f, dope::DopeVector &D, const dope::DopeVector &Ipre, - dope::DopeVector &Ipost) { + dope::DopeVector &Ipost, const Scalar alpha) { if (f.sizeAt(0) == static_cast(0) || - f.sizeAt(0) > D.sizeAt(0)) + f.sizeAt(0) > D.sizeAt(0)) { return; + } if (f.sizeAt(0) == static_cast(1)) { D[0] = f[0]; Ipost[0] = Ipre[0]; @@ -378,7 +433,8 @@ inline void DistanceTransform::distanceL2( --k; // compute horizontal position of intersection between the parabola from q // and the current lowest parabola - s = ((f[q] + q * q) - static_cast(f[v[k]] + v[k] * v[k])) / + s = ((f[q] / alpha + q * q) - + static_cast(f[v[k]] / alpha + v[k] * v[k])) / (static_cast(2 * q) - static_cast(2 * v[k])); } while (s <= z[k]); ++k; @@ -390,9 +446,11 @@ inline void DistanceTransform::distanceL2( k = static_cast(0); for (dope::SizeType q = static_cast(0); q < f.sizeAt(0); ++q) { - while (z[k + 1] < static_cast(q)) ++k; - D[q] = f[v[k]] + - (q - static_cast(v[k])) * (q - static_cast(v[k])); + while (z[k + 1] < static_cast(q)) { + ++k; + } + D[q] = f[v[k]] + alpha * (q - static_cast(v[k])) * + (q - static_cast(v[k])); Ipost[q] = Ipre[v[k]]; } }