From 1f4fda2f96d396856c01161911d4700cc8808872 Mon Sep 17 00:00:00 2001 From: fabien servant Date: Fri, 24 Jul 2020 20:35:20 +0200 Subject: [PATCH 1/3] [software] panoramaCompositing: testing ordered depths --- .../pipeline/main_panoramaCompositing.cpp | 393 ++++++++++++------ 1 file changed, 260 insertions(+), 133 deletions(-) diff --git a/src/software/pipeline/main_panoramaCompositing.cpp b/src/software/pipeline/main_panoramaCompositing.cpp index 2e8ea82c87..c7cf8642f5 100644 --- a/src/software/pipeline/main_panoramaCompositing.cpp +++ b/src/software/pipeline/main_panoramaCompositing.cpp @@ -38,7 +38,6 @@ namespace po = boost::program_options; namespace bpt = boost::property_tree; namespace fs = boost::filesystem; - typedef struct { size_t offset_x; size_t offset_y; @@ -47,6 +46,111 @@ typedef struct { std::string weights_path; } ConfigView; +bool feathering(aliceVision::image::Image & output, const aliceVision::image::Image & color, const aliceVision::image::Image & inputMask) { + + std::vector> feathering; + std::vector> feathering_mask; + feathering.push_back(color); + feathering_mask.push_back(inputMask); + + int lvl = 0; + int width = color.Width(); + int height = color.Height(); + + while (1) { + const image::Image & src = feathering[lvl]; + const image::Image & src_mask = feathering_mask[lvl]; + + image::Image half(width / 2, height / 2); + image::Image half_mask(width / 2, height / 2); + + for (int i = 0; i < half.Height(); i++) { + + int di = i * 2; + for (int j = 0; j < half.Width(); j++) { + int dj = j * 2; + + int count = 0; + half(i, j) = image::RGBfColor(0.0,0.0,0.0); + + if (src_mask(di, dj)) { + half(i, j) += src(di, dj); + count++; + } + + if (src_mask(di, dj + 1)) { + half(i, j) += src(di, dj + 1); + count++; + } + + if (src_mask(di + 1, dj)) { + half(i, j) += src(di + 1, dj); + count++; + } + + if (src_mask(di + 1, dj + 1)) { + half(i, j) += src(di + 1, dj + 1); + count++; + } + + if (count > 0) { + half(i, j) /= float(count); + half_mask(i, j) = 1; + } + else { + half_mask(i, j) = 0; + } + } + + + } + + feathering.push_back(half); + feathering_mask.push_back(half_mask); + + + width = half.Width(); + height = half.Height(); + + if (width < 2 || height < 2) break; + + lvl++; + } + + + for (int lvl = feathering.size() - 2; lvl >= 0; lvl--) { + + image::Image & src = feathering[lvl]; + image::Image & src_mask = feathering_mask[lvl]; + image::Image & ref = feathering[lvl + 1]; + image::Image & ref_mask = feathering_mask[lvl + 1]; + + for (int i = 0; i < src_mask.Height(); i++) { + for (int j = 0; j < src_mask.Width(); j++) { + if (!src_mask(i, j)) { + int mi = i / 2; + int mj = j / 2; + + if (mi >= ref_mask.Height()) { + mi = ref_mask.Height() - 1; + } + + if (mj >= ref_mask.Width()) { + mj = ref_mask.Width() - 1; + } + + src_mask(i, j) = ref_mask(mi, mj); + src(i, j) = ref(mi, mj); + } + } + } + } + + output = feathering[0]; + + return true; +} + void drawBorders(aliceVision::image::Image & inout, aliceVision::image::Image & mask, size_t offset_x, size_t offset_y) { @@ -129,14 +233,14 @@ void drawBorders(aliceVision::image::Image & inout, aliceVisi } } -void drawSeams(aliceVision::image::Image & inout, aliceVision::image::Image & labels) { +void drawSeams(aliceVision::image::Image & inout, aliceVision::image::Image & labels) { for (int i = 1; i < labels.Height() - 1; i++) { for (int j = 1; j < labels.Width() - 1; j++) { - unsigned char label = labels(i, j); - unsigned char same = true; + IndexT label = labels(i, j); + IndexT same = true; same &= (labels(i - 1, j - 1) == label); same &= (labels(i - 1, j + 1) == label); @@ -154,7 +258,7 @@ void drawSeams(aliceVision::image::Image & inout, aliceVision } } -void getMaskFromLabels(aliceVision::image::Image & mask, aliceVision::image::Image & labels, unsigned char index, size_t offset_x, size_t offset_y) { +void getMaskFromLabels(aliceVision::image::Image & mask, aliceVision::image::Image & labels, IndexT index, size_t offset_x, size_t offset_y) { for (int i = 0; i < mask.Height(); i++) { @@ -592,8 +696,8 @@ class LaplacianPyramid { /*Make sure pyramid size can be divided by 2 on each levels*/ double max_scale = 1.0 / pow(2.0, max_levels - 1); - width = size_t(ceil(double(width) * max_scale) / max_scale); - height = size_t(ceil(double(height) * max_scale) / max_scale); + //width = size_t(ceil(double(width) * max_scale) / max_scale); + //height = size_t(ceil(double(height) * max_scale) / max_scale); /*Prepare pyramid*/ for (int lvl = 0; lvl < max_levels; lvl++) { @@ -606,6 +710,92 @@ class LaplacianPyramid { } } + bool augment(size_t new_max_levels) { + + if (new_max_levels <= _levels.size()) { + return false; + } + + size_t old_max_level = _levels.size(); + + image::Image current_color = _levels[_levels.size() - 1]; + image::Image current_weights = _weights[_weights.size() - 1]; + + _levels[_levels.size() - 1].fill(image::RGBfColor(0.0f, 0.0f, 0.0f)); + _weights[_weights.size() - 1].fill(0.0f); + + image::Image current_mask(current_color.Width(), current_color.Height(), true, 0); + image::Image current_color_feathered(current_color.Width(), current_color.Height()); + + for (int i = 0; i < current_color.Height(); i++) { + for (int j = 0; j < current_color.Width(); j++) { + if (current_weights(i, j) < 1e-6) { + current_color(i, j) = image::RGBfColor(0.0); + continue; + } + + current_color(i, j).r() = current_color(i, j).r() / current_weights(i, j); + current_color(i, j).g() = current_color(i, j).g() / current_weights(i, j); + current_color(i, j).b() = current_color(i, j).b() / current_weights(i, j); + current_mask(i ,j) = 255; + } + } + + + feathering(current_color_feathered, current_color, current_mask); + current_color = current_color_feathered; + + + for (int l = old_max_level; l < new_max_levels; l++) { + + _levels.emplace_back(_levels[l - 1].Width() / 2, _levels[l - 1].Height() / 2, true, image::RGBfColor(0.0f, 0.0f, 0.0f)); + _weights.emplace_back(_weights[l - 1].Width() / 2, _weights[l - 1].Height() / 2, true, 0.0f); + } + + int width = current_color.Width(); + int height = current_color.Height(); + image::Image next_color; + image::Image next_weights; + + for (int l = old_max_level - 1; l < new_max_levels - 1; l++) + { + aliceVision::image::Image buf(width, height); + aliceVision::image::Image buf2(width, height); + aliceVision::image::Image bufw(width, height); + + next_color = aliceVision::image::Image(width / 2, height / 2); + next_weights = aliceVision::image::Image(width / 2, height / 2); + + convolveGaussian5x5(buf, current_color); + downscale(next_color, buf); + + convolveGaussian5x5(bufw, current_weights); + downscale(next_weights, bufw); + + upscale(buf, next_color); + convolveGaussian5x5(buf2, buf); + + for (int i = 0; i < buf2.Height(); i++) { + for (int j = 0; j < buf2.Width(); j++) { + buf2(i,j) *= 4.0f; + } + } + + substract(current_color, current_color, buf2); + + merge(current_color, current_weights, l, 0, 0); + + current_color = next_color; + current_weights = next_weights; + width /= 2; + height /= 2; + } + + merge(current_color, current_weights, _levels.size() - 1, 0, 0); + + return true; + } + bool apply(const aliceVision::image::Image & source, const aliceVision::image::Image & weights, size_t offset_x, size_t offset_y) { int width = source.Width(); @@ -749,13 +939,12 @@ class DistanceSeams { public: DistanceSeams(size_t outputWidth, size_t outputHeight) : _weights(outputWidth, outputHeight, true, 0.0f), - _labels(outputWidth, outputHeight, true, 255), - currentIndex(0) + _labels(outputWidth, outputHeight, true, 255) { } virtual ~DistanceSeams() = default; - virtual bool append(const aliceVision::image::Image & inputMask, const aliceVision::image::Image & inputWeights, size_t offset_x, size_t offset_y) + virtual bool append(const aliceVision::image::Image & inputMask, const aliceVision::image::Image & inputWeights, IndexT currentIndex, size_t offset_x, size_t offset_y) { if (inputMask.size() != inputWeights.size()) { return false; @@ -783,19 +972,16 @@ class DistanceSeams { } } - currentIndex++; - return true; } - const image::Image & getLabels() { + const image::Image & getLabels() { return _labels; } private: image::Image _weights; - image::Image _labels; - unsigned char currentIndex; + image::Image _labels; }; class LaplacianCompositer : public Compositer @@ -809,115 +995,27 @@ class LaplacianCompositer : public Compositer } - bool feathering(aliceVision::image::Image & output, const aliceVision::image::Image & color, const aliceVision::image::Image & inputMask) { - - std::vector> feathering; - std::vector> feathering_mask; - feathering.push_back(color); - feathering_mask.push_back(inputMask); - - int lvl = 0; - int width = color.Width(); - int height = color.Height(); - - while (1) { - const image::Image & src = feathering[lvl]; - const image::Image & src_mask = feathering_mask[lvl]; - - image::Image half(width / 2, height / 2); - image::Image half_mask(width / 2, height / 2); - - for (int i = 0; i < half.Height(); i++) { - - int di = i * 2; - for (int j = 0; j < half.Width(); j++) { - int dj = j * 2; - - int count = 0; - half(i, j) = image::RGBfColor(0.0,0.0,0.0); - - if (src_mask(di, dj)) { - half(i, j) += src(di, dj); - count++; - } - - if (src_mask(di, dj + 1)) { - half(i, j) += src(di, dj + 1); - count++; - } - - if (src_mask(di + 1, dj)) { - half(i, j) += src(di + 1, dj); - count++; - } - - if (src_mask(di + 1, dj + 1)) { - half(i, j) += src(di + 1, dj + 1); - count++; - } - - if (count > 0) { - half(i, j) /= float(count); - half_mask(i, j) = 1; - } - else { - half_mask(i, j) = 0; - } - } - - - } - - feathering.push_back(half); - feathering_mask.push_back(half_mask); - + virtual bool append(const aliceVision::image::Image & color, const aliceVision::image::Image & inputMask, const aliceVision::image::Image & inputWeights, size_t offset_x, size_t offset_y) + { + /*Get smalles size*/ + size_t minsize = std::min(color.Height(), color.Width()); - width = half.Width(); - height = half.Height(); - - if (width < 2 || height < 2) break; - - lvl++; - } - + /* + minsize / 2^x = 8 + minsize / 8 = 2^x + x = log2(minsize/8) + */ + size_t optimal_scale = size_t(floor(std::log2(double(minsize) / 8.0))); + if (optimal_scale < _bands) { + ALICEVISION_LOG_ERROR("Decreasing scale !"); + return false; + } - for (int lvl = feathering.size() - 2; lvl >= 0; lvl--) { - - image::Image & src = feathering[lvl]; - image::Image & src_mask = feathering_mask[lvl]; - image::Image & ref = feathering[lvl + 1]; - image::Image & ref_mask = feathering_mask[lvl + 1]; - - for (int i = 0; i < src_mask.Height(); i++) { - for (int j = 0; j < src_mask.Width(); j++) { - if (!src_mask(i, j)) { - int mi = i / 2; - int mj = j / 2; - - if (mi >= ref_mask.Height()) { - mi = ref_mask.Height() - 1; - } - - if (mj >= ref_mask.Width()) { - mj = ref_mask.Width() - 1; - } - - src_mask(i, j) = ref_mask(mi, mj); - src(i, j) = ref(mi, mj); - } - } - } + if (optimal_scale > _bands) { + _bands = optimal_scale; + _pyramid_panorama.augment(_bands); } - output = feathering[0]; - - return true; - } - - virtual bool append(const aliceVision::image::Image & color, const aliceVision::image::Image & inputMask, const aliceVision::image::Image & inputWeights, size_t offset_x, size_t offset_y) - { - aliceVision::image::Image color_big_feathered; - size_t new_offset_x, new_offset_y; aliceVision::image::Image color_pot; aliceVision::image::Image mask_pot; @@ -925,7 +1023,6 @@ class LaplacianCompositer : public Compositer makeImagePyramidCompatible(color_pot, new_offset_x, new_offset_y, color, offset_x, offset_y, _bands); makeImagePyramidCompatible(mask_pot, new_offset_x, new_offset_y, inputMask, offset_x, offset_y, _bands); makeImagePyramidCompatible(weights_pot, new_offset_x, new_offset_y, inputWeights, offset_x, offset_y, _bands); - aliceVision::image::Image feathered; @@ -942,7 +1039,7 @@ class LaplacianCompositer : public Compositer } _pyramid_panorama.apply(feathered, weights_pot, new_offset_x, new_offset_y); - + return true; } @@ -1069,7 +1166,7 @@ int aliceVision_main(int argc, char **argv) bool isMultiBand = false; if (compositerType == "multiband") { - compositer = std::unique_ptr(new LaplacianCompositer(panoramaSize.first, panoramaSize.second, 8)); + compositer = std::unique_ptr(new LaplacianCompositer(panoramaSize.first, panoramaSize.second, 1)); isMultiBand = true; } else if (compositerType == "alpha") @@ -1081,10 +1178,15 @@ int aliceVision_main(int argc, char **argv) compositer = std::unique_ptr(new Compositer(panoramaSize.first, panoramaSize.second)); } + + // Compute seams + std::vector> viewsToDraw; + std::unique_ptr distanceseams(new DistanceSeams(panoramaSize.first, panoramaSize.second)); if (isMultiBand) { + std::map>> indexed_by_scale; for (const auto& viewIt : sfmData.getViews()) { // Load mask @@ -1103,26 +1205,52 @@ int aliceVision_main(int argc, char **argv) image::Image weights; image::readImage(weightsPath, weights, image::EImageColorSpace::NO_CONVERSION); - distanceseams->append(mask, weights, offsetX, offsetY); + distanceseams->append(mask, weights, viewIt.first, offsetX, offsetY); + + /*Get smalles size*/ + size_t minsize = std::min(mask.Height(), mask.Width()); + + /* + minsize / 2^x = 8 + minsize / 8 = 2^x + x = log2(minsize/8) + */ + size_t optimal_scale = size_t(floor(std::log2(double(minsize) / 8.0))); + indexed_by_scale[optimal_scale].push_back(viewIt.second); + } + + for (auto item : indexed_by_scale) { + for (auto view : item.second) { + viewsToDraw.push_back(view); + } + } + } + else { + for (auto& viewIt : sfmData.getViews()) + { + viewsToDraw.push_back(viewIt.second); } } - image::Image labels = distanceseams->getLabels(); + + image::Image labels = distanceseams->getLabels(); distanceseams.reset(); distanceseams = nullptr; oiio::ParamValueList outputMetadata; // Do compositing - int pos = 0; - for (const auto& viewIt : sfmData.getViews()) + for (const auto & view : viewsToDraw) { - if(!sfmData.isPoseAndIntrinsicDefined(viewIt.second.get())) + IndexT viewId = view->getViewId(); + + if(!sfmData.isPoseAndIntrinsicDefined(view.get())) { // skip unreconstructed views continue; } + // Load image and convert it to linear colorspace - const std::string imagePath = (fs::path(warpingFolder) / (std::to_string(viewIt.first) + ".exr")).string(); + const std::string imagePath = (fs::path(warpingFolder) / (std::to_string(viewId) + ".exr")).string(); ALICEVISION_LOG_INFO("Load image with path " << imagePath); image::Image source; image::readImage(imagePath, source, image::EImageColorSpace::NO_CONVERSION); @@ -1137,13 +1265,13 @@ int aliceVision_main(int argc, char **argv) const std::size_t offsetY = metadata.find("AliceVision:offsetY")->get_int(); // Load mask - const std::string maskPath = (fs::path(warpingFolder) / (std::to_string(viewIt.first) + "_mask.exr")).string(); + const std::string maskPath = (fs::path(warpingFolder) / (std::to_string(viewId) + "_mask.exr")).string(); ALICEVISION_LOG_INFO("Load mask with path " << maskPath); image::Image mask; image::readImage(maskPath, mask, image::EImageColorSpace::NO_CONVERSION); // Load Weights - const std::string weightsPath = (fs::path(warpingFolder) / (std::to_string(viewIt.first) + "_weight.exr")).string(); + const std::string weightsPath = (fs::path(warpingFolder) / (std::to_string(viewId) + "_weight.exr")).string(); ALICEVISION_LOG_INFO("Load weights with path " << weightsPath); image::Image weights; image::readImage(weightsPath, weights, image::EImageColorSpace::NO_CONVERSION); @@ -1152,7 +1280,7 @@ int aliceVision_main(int argc, char **argv) if (isMultiBand) { image::Image seams(weights.Width(), weights.Height()); - getMaskFromLabels(seams, labels, pos, offsetX, offsetY); + getMaskFromLabels(seams, labels, viewId, offsetX, offsetY); // Composite image into panorama compositer->append(source, mask, seams, offsetX, offsetY); @@ -1161,7 +1289,6 @@ int aliceVision_main(int argc, char **argv) { compositer->append(source, mask, weights, offsetX, offsetY); } - ++pos; } // Build image From d12031eb559d2f810d6785c40d2c5e86a9a18f80 Mon Sep 17 00:00:00 2001 From: fabien servant Date: Mon, 27 Jul 2020 10:32:43 +0200 Subject: [PATCH 2/3] [software] panoramaCompositing: add some checks to the images validity --- .../pipeline/main_panoramaCompositing.cpp | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/software/pipeline/main_panoramaCompositing.cpp b/src/software/pipeline/main_panoramaCompositing.cpp index c7cf8642f5..2c90f15ed5 100644 --- a/src/software/pipeline/main_panoramaCompositing.cpp +++ b/src/software/pipeline/main_panoramaCompositing.cpp @@ -1137,7 +1137,6 @@ int aliceVision_main(int argc, char **argv) // load input scene sfmData::SfMData sfmData; - std::cout << sfmData.getViews().size() << std::endl; if(!sfmDataIO::Load(sfmData, sfmDataFilepath, sfmDataIO::ESfMData(sfmDataIO::VIEWS|sfmDataIO::EXTRINSICS|sfmDataIO::INTRINSICS))) { ALICEVISION_LOG_ERROR("The input file '" + sfmDataFilepath + "' cannot be read"); @@ -1189,6 +1188,12 @@ int aliceVision_main(int argc, char **argv) std::map>> indexed_by_scale; for (const auto& viewIt : sfmData.getViews()) { + if(!sfmData.isPoseAndIntrinsicDefined(viewIt.second.get())) + { + // skip unreconstructed views + continue; + } + // Load mask const std::string maskPath = (fs::path(warpingFolder) / (std::to_string(viewIt.first) + "_mask.exr")).string(); ALICEVISION_LOG_INFO("Load mask with path " << maskPath); @@ -1228,6 +1233,12 @@ int aliceVision_main(int argc, char **argv) else { for (auto& viewIt : sfmData.getViews()) { + if(!sfmData.isPoseAndIntrinsicDefined(viewIt.second.get())) + { + // skip unreconstructed views + continue; + } + viewsToDraw.push_back(viewIt.second); } } @@ -1298,6 +1309,12 @@ int aliceVision_main(int argc, char **argv) { for (const auto& viewIt : sfmData.getViews()) { + if(!sfmData.isPoseAndIntrinsicDefined(viewIt.second.get())) + { + // skip unreconstructed views + continue; + } + // Load mask const std::string maskPath = (fs::path(warpingFolder) / (std::to_string(viewIt.first) + "_mask.exr")).string(); ALICEVISION_LOG_INFO("Load mask with path " << maskPath); From 18aaa07367c797d03171d069c86ded15d2bb1fa1 Mon Sep 17 00:00:00 2001 From: fabien servant Date: Mon, 27 Jul 2020 12:46:45 +0200 Subject: [PATCH 3/3] [software] panoramaWarping: remove useless code --- .../pipeline/main_panoramaWarping.cpp | 253 ------------------ 1 file changed, 253 deletions(-) diff --git a/src/software/pipeline/main_panoramaWarping.cpp b/src/software/pipeline/main_panoramaWarping.cpp index 96a37759aa..c3c3d96683 100644 --- a/src/software/pipeline/main_panoramaWarping.cpp +++ b/src/software/pipeline/main_panoramaWarping.cpp @@ -38,259 +38,6 @@ namespace po = boost::program_options; namespace bpt = boost::property_tree; namespace fs = boost::filesystem; - -Eigen::VectorXf gaussian_kernel_vector(size_t kernel_length, float sigma) { - - Eigen::VectorXd x; - x.setLinSpaced(kernel_length + 1, -sigma, +sigma); - - Eigen::VectorXd cdf(kernel_length + 1); - for (int i = 0; i < kernel_length + 1; i++) { - cdf(i) = 0.5 * (1 + std::erf(x(i)/sqrt(2.0))); - } - - Eigen::VectorXd k1d(kernel_length); - for (int i = 0; i < kernel_length; i++) { - k1d(i) = cdf(i + 1) - cdf(i); - } - - double sum = k1d.sum(); - k1d = k1d / sum; - - return k1d.cast(); -} - -Eigen::MatrixXf gaussian_kernel(size_t kernel_length, float sigma) { - - Eigen::VectorXf k1d = gaussian_kernel_vector(kernel_length, sigma); - Eigen::MatrixXf K = k1d * k1d.transpose(); - - - double sum = K.sum(); - K = K / sum; - - return K; -} - - - -bool convolve(image::Image & output, const image::Image & input, const image::Image & mask, const Eigen::MatrixXf & kernel) { - - if (output.size() != input.size()) { - return false; - } - - if (output.size() != mask.size()) { - return false; - } - - if (kernel.size() % 2 == 0) { - return false; - } - - if (kernel.rows() != kernel.cols()) { - return false; - } - - int radius = kernel.rows() / 2; - - Eigen::MatrixXf kernel_scaled = kernel; - - for (int i = 0; i < output.Height(); i++) { - for (int j = 0; j < output.Width(); j++) { - - float sum = 0.0f; - float sum_mask = 0.0f; - - if (!mask(i, j)) { - output(i, j) = 0.0f; - continue; - } - - for (int k = 0; k < kernel.rows(); k++) { - float ni = i + k - radius; - if (ni < 0 || ni >= output.Height()) { - continue; - } - - for (int l = 0; l < kernel.cols(); l++) { - float nj = j + l - radius; - if (nj < 0 || nj >= output.Width()) { - continue; - } - - if (!mask(ni, nj)) { - continue; - } - - float val = kernel(k, l) * input(ni, nj); - sum += val; - sum_mask += kernel(k, l); - } - } - - output(i, j) = sum / sum_mask; - } - } - - return true; -} - -bool convolve(image::Image & output, const image::Image & input, const image::Image & mask, const Eigen::MatrixXf & kernel) { - - if (output.size() != input.size()) { - return false; - } - - if (output.size() != mask.size()) { - return false; - } - - if (kernel.size() % 2 == 0) { - return false; - } - - if (kernel.rows() != kernel.cols()) { - return false; - } - - int radius = kernel.rows() / 2; - - Eigen::MatrixXf kernel_scaled = kernel; - - for (int i = 0; i < output.Height(); i++) { - for (int j = 0; j < output.Width(); j++) { - - image::RGBfColor sum(0.0f); - float sum_mask = 0.0f; - - if (!mask(i, j)) { - output(i, j) = sum; - continue; - } - - for (int k = 0; k < kernel.rows(); k++) { - float ni = i + k - radius; - if (ni < 0 || ni >= output.Height()) { - continue; - } - - for (int l = 0; l < kernel.cols(); l++) { - float nj = j + l - radius; - if (nj < 0 || nj >= output.Width()) { - continue; - } - - if (!mask(ni, nj)) { - continue; - } - - sum.r() += kernel(k, l) * input(ni, nj).r(); - sum.g() += kernel(k, l) * input(ni, nj).g(); - sum.b() += kernel(k, l) * input(ni, nj).b(); - sum_mask += kernel(k, l); - } - } - - output(i, j) = sum / sum_mask; - } - } - - return true; -} - -bool computeDistanceMap(image::Image & distance, const image::Image & mask) { - - int m = mask.Height(); - int n = mask.Width(); - - int maxval = m * n; - - distance = image::Image (n, m, false); - for(int x = 0; x < n; ++x) { - - //A corner is when mask becomes 0 - bool b = !mask(0, x); - if (b) { - distance(0, x) = 0; - } - else { - distance(0, x) = maxval * maxval; - } - - for (int y = 1; y < m; y++) { - bool b = !mask(y, x); - if (b) { - distance(y, x) = 0; - } - else { - distance(y, x) = 1 + distance(y - 1, x); - } - } - - for (int y = m - 2; y >= 0; y--) { - if (distance(y + 1, x) < distance(y, x)) { - distance(y, x) = 1 + distance(y + 1, x); - } - } - } - - for (int y = 0; y < m; y++) { - int q; - std::map s; - std::map t; - - q = 0; - s[0] = 0; - t[0] = 0; - - std::function f = [distance, y](int x, int i) { - int gi = distance(y, i); - return (x - i)*(x - i) + gi * gi; - }; - - std::function sep = [distance, y](int i, int u) { - int gu = distance(y, u); - int gi = distance(y, i); - - int nom = (u * u) - (i * i) + (gu * gu) - (gi * gi); - int denom = 2 * (u - i); - - return nom / denom; - }; - - for (int u = 1; u < n; u++) { - - while (q >= 0 && (f(t[q], s[q]) > f(t[q], u))) { - q = q - 1; - } - - if (q < 0) { - q = 0; - s[0] = u; - } - else { - int w = 1 + sep(s[q], u); - if (w < n) { - q = q + 1; - s[q] = u; - t[q] = w; - } - } - } - - for (int u = n - 1; u >= 0; u--) { - distance(y, u) = f(u, s[q]); - if (u == t[q]) { - q = q - 1; - } - } - } - - return true; -} - - namespace SphericalMapping { /**