Merge pull request #41 from ebertin/feature/dft_background

Use DFT convolution for kernel > 5x5

Author: marcschefer

SEBenchmarks/CMakeLists.txt

@@ -15,6 +15,7 @@ elements_subdir(SEBenchmarks)
 #===============================================================================
 elements_depends_on_subdirs(ModelFitting)
 elements_depends_on_subdirs(SEFramework)
+elements_depends_on_subdirs(SEImplementation)
 
 #===============================================================================
 # Add the find_package macro (a pure CMake command) here to locate the
@@ -48,6 +49,8 @@ endif ()
 #===============================================================================
 elements_add_executable(BenchConvolution src/program/BenchConvolution.cpp
         LINK_LIBRARIES ModelFitting SEFramework ${Boost_LIBRARIES} ${OpenCV_LIBS})
+elements_add_executable(BenchBackgroundConvolution src/program/BenchBackgroundConvolution.cpp
+        LINK_LIBRARIES SEFramework SEImplementation ${Boost_LIBRARIES})
 
 #===============================================================================
 # Declare the Boost tests here

SEBenchmarks/src/program/BenchBackgroundConvolution.cpp

@@ -0,0 +1,170 @@
+/*
+ * Copyright (C) 2012-2020 Euclid Science Ground Segment
+ *
+ * This library is free software; you can redistribute it and/or modify it under
+ * the terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation; either version 3.0 of the License, or (at your option)
+ * any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this library; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/**
+ * @file src/program/BenchBackgroundConvolution.cpp
+ * @date 27/03/19
+ * @author aalvarez
+ */
+
+#include <map>
+#include <string>
+
+#include <boost/program_options.hpp>
+#include <boost/timer/timer.hpp>
+#include <random>
+#include "ElementsKernel/ProgramHeaders.h"
+#include "ElementsKernel/Real.h"
+#include "SEImplementation/Segmentation/BackgroundConvolution.h"
+#include "SEFramework/Image/VectorImage.h"
+#include "SEUtils/IsClose.h"
+
+#if BOOST_VERSION < 105600
+#include <boost/units/detail/utility.hpp>
+using boost::units::detail::demangle;
+#else
+using boost::core::demangle;
+#endif
+
+namespace po = boost::program_options;
+namespace timer = boost::timer;
+using namespace SExtractor;
+
+static Elements::Logging logger = Elements::Logging::getLogger("BenchBackgroundConvolution");
+
+class BenchBackgroundConvolution : public Elements::Program {
+private:
+  std::default_random_engine random_generator;
+  std::uniform_real_distribution<SeFloat> random_dist{0, 1};
+
+public:
+
+  po::options_description defineSpecificProgramOptions() override {
+    po::options_description options{};
+    options.add_options()
+      ("image-start", po::value<int>()->default_value(100), "Image start size")
+      ("image-step-size", po::value<int>()->default_value(3), "Image step size")
+      ("image-nsteps", po::value<int>()->default_value(1), "Number of steps for the image")
+      ("kernel-start", po::value<int>()->default_value(3), "Kernel start size")
+      ("kernel-step-size", po::value<int>()->default_value(4), "Kernel step size")
+      ("kernel-nsteps", po::value<int>()->default_value(2), "Number of steps for the kernel")
+      ("repeat", po::value<int>()->default_value(5), "Repeat")
+      ("measures", po::value<int>()->default_value(10), "Number of measures");
+    return options;
+  }
+
+  std::shared_ptr<VectorImage<SeFloat>> generateImage(int size) {
+    auto img = VectorImage<SeFloat>::create(size, size);
+    for (int x = 0; x < size; ++x) {
+      for (int y = 0; y < size; ++y) {
+        img->setValue(x, y, random_dist(random_generator));
+      }
+    }
+    return img;
+  }
+
+  Elements::ExitCode mainMethod(std::map<std::string, po::variable_value>& args) override {
+
+    auto img_start = args["image-start"].as<int>();
+    auto img_step_size = args["image-step-size"].as<int>();
+    auto img_nsteps = args["image-nsteps"].as<int>();
+    auto krn_start = args["kernel-start"].as<int>();
+    auto krn_step_size = args["kernel-step-size"].as<int>();
+    auto krn_nsteps = args["kernel-nsteps"].as<int>();
+    auto repeat = args["repeat"].as<int>();
+    auto measures = args["measures"].as<int>();
+
+    std::cout << "Image,Kernel,Implementation,Time" << std::endl;
+
+    for (int img_step = 0; img_step < img_nsteps; ++img_step) {
+      auto img_size = img_start + img_step * img_step_size;
+      auto image = generateImage(img_size);
+      auto variance = generateImage(img_size);
+
+      for (int krn_step = 0; krn_step < krn_nsteps; ++krn_step) {
+        auto krn_size = krn_start + krn_step * krn_step_size;
+
+        logger.info() << "Using an image of " << img_size << "x" << img_size;
+        logger.info() << "Using a kernel of " << krn_size << "x" << krn_size;
+
+        auto kernel = generateImage(krn_size);
+
+        logger.info() << "Timing Direct implementation";
+        auto direct_result = benchmark<BgConvolutionImageSource>(image, variance, kernel, repeat, measures);
+
+        logger.info() << "Timing DFT implementation";
+        auto dft_result = benchmark<BgDFTConvolutionImageSource>(image, variance, kernel, repeat, measures);
+
+        logger.info() << "Comparing results";
+        verifyResults(direct_result, dft_result);
+      }
+    }
+
+    return Elements::ExitCode::OK;
+  }
+
+  template<typename BackgroundConvolution>
+  std::shared_ptr<VectorImage<SeFloat>>
+  benchmark(std::shared_ptr<VectorImage<SeFloat>>& image, std::shared_ptr<VectorImage<SeFloat>>& variance,
+            std::shared_ptr<VectorImage<SeFloat>>& kernel, int repeat, int measures) {
+    auto conv_name = demangle(typeid(BackgroundConvolution).name());
+
+    auto bg_convolution = std::make_shared<BackgroundConvolution>(image, variance, 0.5, kernel);
+
+    std::shared_ptr<VectorImage<SeFloat>> result;
+
+    for (int m = 0; m < measures; ++m) {
+      logger.info() << conv_name << " " << m + 1 << "/" << measures;
+      timer::cpu_timer timer;
+      timer.stop();
+
+      for (int r = 0; r < repeat; ++r) {
+        timer.start();
+        result = bg_convolution->getImageTile(0, 0, image->getWidth(), image->getHeight())->getImage();
+        timer.stop();
+      }
+
+      std::cout << image->getWidth() << ',' << kernel->getWidth() << ",\"" << conv_name << "\"," << timer.elapsed().wall
+                << std::endl;
+    }
+
+    return result;
+  }
+
+  void verifyResults(std::shared_ptr<VectorImage<SeFloat>> a, std::shared_ptr<VectorImage<SeFloat>> b) {
+    bool all_equal = true;
+    for (int x = 0; x < a->getWidth(); ++x) {
+      for (int y = 0; y < a->getHeight(); ++y) {
+        auto av = a->getValue(x, y);
+        auto bv = b->getValue(x, y);
+        if (!isClose(av, bv, 1e-3, 1e-3)) {
+          logger.info() << "Mismatch at " << x << ',' << y << ": "
+                        << av << " != " << bv;
+          all_equal = false;
+        }
+      }
+    }
+    if (all_equal) {
+      logger.info() << "All elements are equal!";
+    } else {
+      logger.warn() << "Convoluted images are not equal!";
+    }
+  }
+};
+
+MAIN_FOR(BenchBackgroundConvolution)

SEBenchmarks/src/program/BenchConvolution.cpp

@@ -30,6 +30,7 @@
 #include <random>
 #include "ElementsKernel/ProgramHeaders.h"
 #include "ElementsKernel/Real.h"
+#include "SEUtils/IsClose.h"
 #include "SEFramework/Image/VectorImage.h"
 #include "SEFramework/Convolution/DirectConvolution.h"
 #include "SEFramework/Convolution/DFT.h"
@@ -180,11 +181,6 @@ class BenchConvolution : public Elements::Program {
       logger.warn() << "Convoluted images are not equal!";
     }
   }
-
-  bool isClose(SeFloat a, SeFloat b, SeFloat atol=1e-8, SeFloat rtol=1e-5) const {
-    return std::abs(a - b) <= (atol + rtol * std::abs(b));
-  }
-
 };
 
 MAIN_FOR(BenchConvolution)

SEFramework/CMakeLists.txt

@@ -133,6 +133,9 @@ elements_add_unit_test(RecenterImage_test tests/src/Image/RecenterImage_test.cpp
 elements_add_unit_test(MirrorImage_test tests/src/Image/MirrorImage_test.cpp
                      LINK_LIBRARIES SEFramework
                      TYPE Boost)
+elements_add_unit_test(FunctionalImage_test tests/src/Image/FunctionalImage_test.cpp
+                     LINK_LIBRARIES SEFramework
+                     TYPE Boost)
 elements_add_unit_test(FFT_test tests/src/FFT/FFT_test.cpp
                      LINK_LIBRARIES SEFramework
                      TYPE Boost)

SEFramework/SEFramework/Convolution/DFT.h

@@ -39,7 +39,8 @@ class DFTConvolution {
     return m_kernel->getHeight();
   }
 
-  void convolve(std::shared_ptr<WriteableImage<T>> image_ptr) const {
+  template <typename ...Args>
+  void convolve(std::shared_ptr<WriteableImage<T>> image_ptr, Args... padding_args) const {
     // Dimension of the working padded images
     auto padded_width = image_ptr->getWidth() + m_kernel->getWidth() - 1;
     auto padded_height = image_ptr->getHeight() + m_kernel->getHeight() - 1;
@@ -51,7 +52,7 @@ class DFTConvolution {
     auto total_size = padded_height * padded_width;
 
     // Padded image
-    auto padded = TPadding::create(image_ptr, padded_width, padded_height);
+    auto padded = TPadding::create(image_ptr, padded_width, padded_height, std::forward<Args>(padding_args)...);
 
     // Buffers
     std::vector<real_t> real_buffer(total_size * 2);

SEFramework/SEFramework/Convolution/DirectConvolution.h

@@ -21,13 +21,16 @@ class DirectConvolution {
 
   virtual ~DirectConvolution() = default;
 
-  void convolve(std::shared_ptr<WriteableImage<T>> image) const {
+  template <typename ...Args>
+  void convolve(std::shared_ptr<WriteableImage<T>> image, Args... padding_args) const {
     auto padded_width = image->getWidth() + m_kernel->getWidth() - 1;
     auto padded_height = image->getHeight() + m_kernel->getHeight() - 1;
     auto lpad = m_kernel->getWidth() / 2;
     auto tpad = m_kernel->getHeight() / 2;
 
-    auto padded = VectorImage<T>::create(TPadding::create(image, padded_width, padded_height));
+    auto padded = VectorImage<T>::create(
+      TPadding::create(image, padded_width, padded_height, std::forward<Args>(padding_args)...)
+    );
 
     for (int iy = tpad; iy < padded->getHeight() - tpad; ++iy) {
       for (int ix = lpad; ix < padded->getWidth() - lpad; ++ix) {

SEFramework/SEFramework/Image/BufferedImage.h

@@ -39,16 +39,13 @@ class BufferedImage : public ImageBase<T> {
 
   /// Returns the value of the pixel with the coordinates (x,y)
   virtual T getValue(int x, int y) const override {
-    //std::cout << "BufferedImage::getValue() " << x << " " << y << "\n";
     assert(x >= 0 && y >=0 && x < m_source->getWidth() && y < m_source->getHeight());
 
-    auto& current_tile = getEntryForThisThread();
-
-    if (current_tile == nullptr || !current_tile->isPixelInTile(x, y)) {
-      current_tile = m_tile_manager->getTileForPixel(x, y, m_source);
+    if (m_current_tile == nullptr || !m_current_tile->isPixelInTile(x, y)) {
+      m_current_tile = m_tile_manager->getTileForPixel(x, y, m_source);
     }
 
-    return current_tile->getValue(x, y);
+    return m_current_tile->getValue(x, y);
   }
 
   /// Returns the width of the image in pixels
@@ -80,16 +77,7 @@ class BufferedImage : public ImageBase<T> {
 protected:
   std::shared_ptr<const ImageSource<T>> m_source;
   std::shared_ptr<TileManager> m_tile_manager;
-
-  std::shared_ptr<ImageTile<T>> &getEntryForThisThread() const {
-    if (!m_current_tile.get()) {
-      m_current_tile.reset(new std::shared_ptr<ImageTile<T>>{nullptr});
-    }
-    return *m_current_tile;
-  }
-
-private:
-  mutable boost::thread_specific_ptr<std::shared_ptr<ImageTile<T>>> m_current_tile;
+  mutable std::shared_ptr<ImageTile<T>> m_current_tile;
 };
 
 }

SEFramework/SEFramework/Image/FunctionalImage.h

@@ -0,0 +1,51 @@
+/*
+ * @file SEFramework/Image/FunctionalImage.h
+ * @date 27/03/19
+ * @author Alejandro Alvarez Ayllon
+ */
+
+#ifndef _SEFRAMEWORK_IMAGE_FUNCTIONALIMAGE_H
+#define _SEFRAMEWORK_IMAGE_FUNCTIONALIMAGE_H
+
+#include "SEFramework/Image/ImageBase.h"
+
+namespace SExtractor {
+
+template<typename T>
+class FunctionalImage : public ImageBase<T> {
+public:
+  using FunctorType = std::function<T(int x, int y)>;
+
+protected:
+  FunctionalImage(int width, int height, FunctorType functor)
+    : m_width{width}, m_height{height}, m_functor{functor} {
+  }
+
+public:
+  virtual ~FunctionalImage() = default;
+
+  template<typename ...Args>
+  static std::shared_ptr<ImageBase<T>> create(Args&& ... args) {
+    return std::shared_ptr<FunctionalImage<T>>(new FunctionalImage{std::forward<Args>(args)...});
+  }
+
+  T getValue(int x, int y) const override {
+    return m_functor(x, y);
+  }
+
+  int getWidth() const override {
+    return m_width;
+  }
+
+  int getHeight() const override {
+    return m_height;
+  }
+
+private:
+  int m_width, m_height;
+  FunctorType m_functor;
+};
+
+} // end SExtractor
+
+#endif // _SEFRAMEWORK_IMAGE_FUNCTIONALIMAGE_H

SEFramework/SEFramework/Image/WriteableBufferedImage.h

@@ -21,7 +21,7 @@ class WriteableBufferedImage : public BufferedImage<T>, public WriteableImage<T>
   WriteableBufferedImage(std::shared_ptr<const ImageSource<T>> source, std::shared_ptr<TileManager> tile_manager)
       : BufferedImage<T>(source, tile_manager) {}
 
-  using BufferedImage<T>::getEntryForThisThread;
+  using BufferedImage<T>::m_current_tile;
 
 public:
 
@@ -35,15 +35,12 @@ class WriteableBufferedImage : public BufferedImage<T>, public WriteableImage<T>
   virtual void setValue(int x, int y, T value) override {
     assert(x >= 0 && y >=0 && x < BufferedImage<T>::m_source->getWidth() && y < BufferedImage<T>::m_source->getHeight());
 
-    auto& current_tile = getEntryForThisThread();
-
-    if (current_tile == nullptr || !current_tile->isPixelInTile(x, y)) {
-      current_tile =
-          BufferedImage<T>::m_tile_manager->getTileForPixel(x, y, BufferedImage<T>::m_source);
+    if (m_current_tile == nullptr || !m_current_tile->isPixelInTile(x, y)) {
+      m_current_tile = BufferedImage<T>::m_tile_manager->getTileForPixel(x, y, BufferedImage<T>::m_source);
     }
 
-    current_tile->setModified(true);
-    current_tile->setValue(x, y, value);
+    m_current_tile->setModified(true);
+    m_current_tile->setValue(x, y, value);
   }
 
 };