Add remap functions using coordinates maps

Author: Fabien Servant

src/aliceVision/image/CMakeLists.txt

@@ -19,6 +19,7 @@ set(image_files_headers
     pixelTypes.hpp
     Rgb.hpp
     Sampler.hpp
+    remap.hpp
     cache.hpp
 )
 
@@ -31,6 +32,7 @@ set(image_files_sources
     io.cpp
     imageAlgo.cpp
     jetColorMap.cpp
+    remap.cpp
     cache.cpp
 )
 

src/aliceVision/image/remap.cpp

@@ -0,0 +1,15 @@
+// This file is part of the AliceVision project.
+// Copyright (c) 2025 AliceVision contributors.
+// This Source Code Form is subject to the terms of the Mozilla Public License,
+// v. 2.0. If a copy of the MPL was not distributed with this file,
+// You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include <aliceVision/image/remap.hpp>
+
+namespace aliceVision {
+namespace image {
+
+
+
+}  // namespace image
+}  // namespace aliceVision

src/aliceVision/image/remap.hpp

@@ -0,0 +1,200 @@
+// This file is part of the AliceVision project.
+// Copyright (c) 2025 AliceVision contributors.
+// This Source Code Form is subject to the terms of the Mozilla Public License,
+// v. 2.0. If a copy of the MPL was not distributed with this file,
+// You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <aliceVision/image/Image.hpp>
+#include <aliceVision/image/Sampler.hpp>
+#include <aliceVision/system/Logger.hpp>
+#include <type_traits>
+
+namespace aliceVision {
+namespace image {
+
+/**
+ * @brief Fill pixels of an image such that dest(i, j) = input(map(i, j).y, map(i, j).x)
+ * @param[in] src source image (of size H1,W1) to remap
+ * @param[in] map image (of size H2, W2) containing coordinates in the source image
+ * @param[in] fillColor the default pixel value if not possible to remap
+ * @param[out] output image to store the result (will be resized to size H2, W2)
+ */
+template <typename P>
+void remap(const Image<P> & source, const Image<Vec2> & map, const P & fillColor, Image<P> & output)
+{
+    if (output.width() != map.width() || output.height() != map.height())
+    {
+        output.resize(map.width(), map.height());
+    }
+
+    const image::Sampler2d<image::SamplerLinear> sampler;
+
+    for (int i = 0; i < map.height(); i++)
+    {
+        for (int j = 0; j < map.width(); j++)
+        {
+            const Vec2 & coord = map(i, j);
+
+            const double & x = coord.x();
+            const double & y = coord.y();
+
+            // pick pixel if it is in the image domain
+            if (source.contains(y, x))
+            {
+                output(i, j) = sampler(source, y, x);
+            }
+            else 
+            {
+                output(i, j) = fillColor;
+            }
+        }
+    }
+}
+
+/**
+ * @brief Fill pixels of an image such that dest(i, j) = input(map(i, j).y, map(i, j).x)
+ * @param[in] src source image (of size H1,W1) to remap
+ * @param[in] map image (of size H2, W2) containing coordinates in the source image
+ * @param[in] fillColor the default pixel value if not possible to remap
+ * @param[out] output image to store the result (will be resized to size H2, W2)
+ */
+template <typename P>
+void remapInter(const Image<P> & source, const Image<Vec2> & map, const P& fillColor, Image<P> & output)
+{
+    if (output.width() != map.width() || output.height() != map.height())
+    {
+        output.resize(map.width(), map.height());
+    }
+
+    int srcWidth = source.width();
+    int srcHeight = source.height();
+
+    const image::Sampler2d<image::SamplerLinear> sampler;
+
+    for (int i = 0; i < map.height(); i++)
+    {
+        for (int j = 0; j < map.width(); j++)
+        {
+            const Vec2 & coord = map(i, j);
+
+            const double & x = coord.x();
+            const double & y = coord.y();
+            
+            double scaleX = 1.0;
+            if (j < map.width() - 1 && j > 0) 
+            {
+                scaleX = std::abs(map(i, j + 1).x() - map(i, j - 1).x()) * 0.5;
+            } 
+            else if (j < map.width() - 1) 
+            {
+                scaleX = std::abs(map(i, j + 1).x() - x);
+            } 
+            else if (j > 0) 
+            {
+                scaleX = std::abs(x - map(i, j - 1).x());
+            }
+            
+            double scaleY = 1.0;
+            if (i < map.height() - 1 && i > 0) 
+            {
+                scaleY = std::abs(map(i + 1, j).y() - map(i - 1, j).y()) * 0.5;
+            } 
+            else if (i < map.height() - 1) 
+            {
+                scaleY = std::abs(map(i + 1, j).y() - y);
+            } 
+            else if (i > 0) 
+            {
+                scaleY = std::abs(y - map(i - 1, j).y());
+            }
+        
+            const double max_scale = std::max(scaleX, scaleY);
+            if (scaleX > 1.0 && scaleY > 1.0)
+            {
+                //We are doing downsampling here.
+                //Use area interpolation to avoid aliasing
+
+                double halfWidth = scaleX * 0.5;
+                double halfHeight = scaleY * 0.5;
+
+                double x1 = coord.x() - halfWidth;
+                double x2 = coord.x() + halfWidth;
+                double y1 = coord.y() - halfHeight;
+                double y2 = coord.y() + halfHeight;
+
+                int ix1 = static_cast<int>(x1);
+                int ix2 = static_cast<int>(std::ceil(x2));
+                int iy1 = static_cast<int>(y1);
+                int iy2 = static_cast<int>(std::ceil(y2));
+
+                //Not using sampler as the neigborhood size is dynamic
+
+                P sum(0.0);
+                double wsum = 0.0;
+
+                for (int by = iy1; by < iy2; by++)
+                {
+                    if (by < 0 || by >= srcHeight)
+                    {
+                        continue;
+                    }
+
+                    //Compute vertical occupancy
+                    // (by + 1) - y1
+                    // or y2 - by
+                    double yocc = std::max(0.0, std::min(double(by + 1), y2) - std::max(double(by), y1));
+
+                    for (int bx = ix1; bx < ix2; bx++)
+                    {
+                        if (bx < 0 || bx >= srcWidth)
+                        {
+                            continue;
+                        }
+
+                        //Compute horizontal occupancy
+                        // (bx + 1) - x1
+                        // or x2 - bx
+                        double xocc = std::max(0.0, std::min(double(bx + 1), x2) - std::max(double(bx), x1));
+                        double weight = yocc * xocc;
+
+                        sum += source(by, bx) * weight;
+                        wsum += weight;
+                    }
+                }
+
+                //Output pixel is the weighted mean of the 
+                //Source pixels in the area covered by the output pixel.
+                if (wsum > 0.0)
+                {
+                    sum /= wsum;
+
+                    if constexpr (std::is_same<P, RGBAfColor>())
+                    {
+                        sum.a() = 1.0;
+                    }
+                }
+
+                output(i, j) = sum;
+            }
+            else 
+            {
+                // Not downsampling, do bilinear interpolation
+
+                // pick pixel if it is in the image domain
+                if (source.contains(y, x))
+                {
+                    output(i, j) = sampler(source, y, x);
+                }
+                else 
+                {
+                    output(i, j) = fillColor;
+                }
+            }
+        }
+    }
+}
+
+}  // namespace image
+}  // namespace aliceVision