Improve debug log

.github/workflows/dependencies.txt

@@ -17,3 +17,4 @@ levmar-devel
 log4cpp-devel
 onnxruntime-devel
 wcslib-devel
+opencv-devel

SEBenchmarks/src/program/BenchConvolution.cpp

@@ -106,7 +106,7 @@ class BenchConvolution : public Elements::Program {
 
 #ifdef WITH_OPENCV
         logger.info() << "Timing OpenCV implementation";
-        auto opencv_result = benchmark<OpenCVConvolution>(image, kernel, repeat, measures);
+        auto opencv_result = benchmark<OpenCVConvolution<SeFloat>>(image, kernel, repeat, measures);
 #endif
 
         if (krn_size <= 10 || img_size <= 20) {

SEFramework/CMakeLists.txt

@@ -33,6 +33,7 @@ find_package(CCfits)
 find_package(BoostDLL)
 find_package(FFTW COMPONENTS single double)
 find_package(WCSLIB REQUIRED)
+find_package(OpenCV REQUIRED)
 
 
 #===============================================================================
@@ -59,8 +60,8 @@ elements_add_library(SEFramework
                      src/lib/CoordinateSystem/*.cpp
                      LINK_LIBRARIES
                         ElementsKernel SEUtils Table Configuration MathUtils FilePool
-                        WCSLIB ${CCFITS_LIBRARIES} ${FFTW_LIBRARIES}
-                     INCLUDE_DIRS SEUtils ${CCFITS_INCLUDE_DIRS} ${BoostDLL_INCLUDE_DIRS} ${FFTW_INCLUDE_DIRS}
+			WCSLIB ${CCFITS_LIBRARIES} ${FFTW_LIBRARIES} ${OpenCV_LIBRARIES}
+			INCLUDE_DIRS SEUtils ${CCFITS_INCLUDE_DIRS} ${BoostDLL_INCLUDE_DIRS} ${FFTW_INCLUDE_DIRS} ${OpenCV_INCLUDE_DIRS}
                      PUBLIC_HEADERS SEFramework)
 
 #===============================================================================

SEFramework/SEFramework/Convolution/OpenCVConvolution.h

@@ -28,24 +28,38 @@
 
 namespace SourceXtractor {
 
+template <typename T = SeFloat>
 class OpenCVConvolution {
 public:
-  explicit OpenCVConvolution(std::shared_ptr<const VectorImage<SeFloat>> img)
-    : m_kernel(img->getWidth(), img->getHeight(), CV_32F) {
+  explicit OpenCVConvolution(std::shared_ptr<const VectorImage<T>> img)
+    : m_kernel(img->getWidth(), img->getHeight(), CV_32F) , m_original_kernel{VectorImage<T>::create(*MirrorImage<T>::create(img))} {
     cv::Mat aux(img->getWidth(), img->getHeight(), CV_32F);
     std::copy(img->getData().begin(), img->getData().end(), aux.begin<float>());
     cv::flip(aux, m_kernel, -1);
   }
 
   virtual ~OpenCVConvolution() = default;
 
-  void convolve(std::shared_ptr<VectorImage<SeFloat>> image) const {
+  void convolve(std::shared_ptr<WriteableImage<T>> image) const {
+    auto vector_image = VectorImage<T>::create(image);
+
     cv::Mat image_cv (image->getHeight(), image->getWidth(), CV_32F);
-    std::copy(image->getData().begin(), image->getData().end(), image_cv.begin<float>());
+    // for (int y = 0; y < image->getHeight(); ++y) {
+    //   for (int x = 0; x < image->getWidth(); ++x) {
+    //     image_cv.at<float>(y, x) = image->getValue(x, y);
+    //   }
+    // }
+    std::copy(vector_image->getData().begin(), vector_image->getData().end(), image_cv.begin<float>());
 
     cv::filter2D(image_cv, image_cv, -1, m_kernel);
 
-    std::copy(image_cv.begin<float>(), image_cv.end<float>(), image->getData().begin());
+    for (int y = 0; y < image->getHeight(); ++y) {
+      for (int x = 0; x < image->getWidth(); ++x) {
+        image->setValue(x, y, image_cv.at<float>(y, x));
+      }
+    }
+
+    //std::copy(image_cv.begin<float>(), image_cv.end<float>(), image->getData().begin());
   }
 
   std::size_t getWidth() const {
@@ -56,8 +70,13 @@ class OpenCVConvolution {
     return m_kernel.size[1];
   }
 
+  std::shared_ptr<const Image<T>> getKernel() const {
+    return m_original_kernel;
+  }
+
 private:
   cv::Mat m_kernel;// (size, size, CV_32F);
+  std::shared_ptr<const VectorImage<T>> m_original_kernel;
 };
 
 } // end of SourceXtractor

SEFramework/SEFramework/Pipeline/SourceGrouping.h

@@ -121,6 +121,7 @@ class SourceGrouping : public SourceGroupingInterface {
   std::shared_ptr<SourceGroupFactory> m_group_factory;
   std::list<std::unique_ptr<SourceGroupInterface>> m_source_groups;
   unsigned int m_hard_limit;
+  unsigned int m_total_sources_waiting = 0;
 
 }; /* End of SourceGrouping class */
 

SEFramework/src/lib/FFT/FFT.cpp

@@ -26,6 +26,8 @@
 #include <fftw3.h>
 #include <map>
 
+#include <iostream>
+
 namespace SourceXtractor {
 
 /**
@@ -130,6 +132,9 @@ auto FFT<T>::createForwardPlan(int width, int height, std::vector<T>& inout) ->
   }
 
   // No available plan yet, so get one from FFTW
+
+  std::cout << "!!!!!!!!!!Creating new FFTW forward plan for size " << width << "x" << height << std::endl;
+
   boost::upgrade_to_unique_lock<boost::shared_mutex> write_lock{read_lock};
   boost::lock_guard<boost::mutex>                    lock_planner{fftw_global_plan_mutex};
 
@@ -162,14 +167,21 @@ auto FFT<T>::createInversePlan(int width, int height, std::vector<T>& inout) ->
   static boost::shared_mutex                        mutex;
   static std::map<std::tuple<int, int>, plan_ptr_t> plan_cache;
 
+
+  auto before_time = std::chrono::high_resolution_clock::now();
   boost::upgrade_lock<boost::shared_mutex> read_lock{mutex};
+  auto after_time =  std::chrono::high_resolution_clock::now();
+  std::cout << "!!!!Time to get read lock: " << std::chrono::duration_cast<std::chrono::microseconds>(after_time - before_time).count() << " us" << std::endl;
 
   auto pi = plan_cache.find(std::make_tuple(width, height));
   if (pi != plan_cache.end()) {
     return pi->second;
   }
 
   // No available plan yet, so get one from FFTW
+
+  std::cout << "!!!!!!!!!!Creating new FFTW inverse plan for size " << width << "x" << height << std::endl;
+
   boost::upgrade_to_unique_lock<boost::shared_mutex> write_lock{read_lock};
   boost::lock_guard<boost::mutex>                    lock_planner{fftw_global_plan_mutex};
 

SEFramework/src/lib/Pipeline/SourceGrouping.cpp

@@ -21,18 +21,25 @@
  */
 
 #include "SEFramework/Pipeline/SourceGrouping.h"
+#include <ElementsKernel/Logging.h>
 #include <vector>
 
 
 namespace SourceXtractor {
 
+static Elements::Logging logger = Elements::Logging::getLogger("SourceGrouping");
+
+
 SourceGrouping::SourceGrouping(std::shared_ptr<GroupingCriteria> grouping_criteria,
                                std::shared_ptr<SourceGroupFactory> group_factory,
                                unsigned int hard_limit)
         : m_grouping_criteria(grouping_criteria), m_group_factory(group_factory), m_hard_limit(hard_limit) {
 }
 
 void SourceGrouping::receiveSource(std::unique_ptr<SourceInterface> source) {
+  m_total_sources_waiting++;
+  logger.debug() << "Receiving source, sources in grouping: " << m_total_sources_waiting;
+
   // Pointer which points to the group of the source
   SourceGroupInterface* matched_group = nullptr;
 
@@ -87,6 +94,8 @@ void SourceGrouping::receiveSource(std::unique_ptr<SourceInterface> source) {
 void SourceGrouping::receiveProcessSignal(const ProcessSourcesEvent& process_event) {
   std::vector<std::list<std::unique_ptr<SourceGroupInterface>>::iterator> groups_to_process;
 
+  logger.debug() << "Received processing signal, total sources waiting in grouping: " << m_total_sources_waiting;
+
   // We iterate through all the SourceGroups we have
   for (auto group_it = m_source_groups.begin(); group_it != m_source_groups.end(); ++group_it) {
     // We look at its Sources and if we find at least one that needs to be processed we put it in groups_to_process
@@ -101,9 +110,13 @@ void SourceGrouping::receiveProcessSignal(const ProcessSourcesEvent& process_eve
   // For each SourceGroup that we put in groups_to_process,
   for (auto& group : groups_to_process) {
     // we remove it from our list of stored SourceGroups and notify our observers
+    m_total_sources_waiting -= (*group)->size();
+    logger.debug() << "Sending group size " << (*group)->size() << ", sources remaining in grouping: " << m_total_sources_waiting;
     sendSource(std::move(*group));
     m_source_groups.erase(group);
   }
+
+  logger.debug() << "Processing signal handled, total sources remaining in grouping: " << m_total_sources_waiting;
 }
 
 std::set<PropertyId> SourceGrouping::requiredProperties() const {

SEFramework/src/lib/Source/EntangledSource.cpp

@@ -21,9 +21,12 @@
 
 #include "SEFramework/Source/SourceGroupWithOnDemandProperties.h"
 #include "SEFramework/Task/GroupTask.h"
+#include "ElementsKernel/Logging.h"
 
 namespace SourceXtractor {
 
+static Elements::Logging logger = Elements::Logging::getLogger("EntangledSource");
+
 SourceGroupWithOnDemandProperties::EntangledSource::EntangledSource(std::shared_ptr<SourceInterface> source, SourceGroupWithOnDemandProperties& group)
         : m_source(source), m_group(group) {
   // Normally, it should not be possible that the given source is of type
@@ -38,7 +41,6 @@ SourceGroupWithOnDemandProperties::EntangledSource::EntangledSource(std::shared_
 }
 
 const Property& SourceGroupWithOnDemandProperties::EntangledSource::getProperty(const PropertyId& property_id) const {
-
   // If we already have the property stored in this object, returns it
   if (m_property_holder.isPropertySet(property_id)) {
     return m_property_holder.getProperty(property_id);
@@ -63,7 +65,8 @@ const Property& SourceGroupWithOnDemandProperties::EntangledSource::getProperty(
       throw PropertyNotFoundException(property_id);
     }
 
-  // Use the task to make the property
+    //logger.debug() << "Computing property " << property_id.getString();
+    // Use the task to make the property
     group_task->computeProperties(m_group);
 
     // The property should now be available either in this object or in the group object

SEFramework/src/lib/Source/SourceWithOnDemandProperties.cpp

@@ -46,6 +46,7 @@ const Property& SourceWithOnDemandProperties::getProperty(const PropertyId& prop
     // if not, get the task that makes it and execute, we should have it then
     auto task = m_task_provider->getTask<SourceTask>(property_id);
     if (task) {
+      //logger.debug() << "Computing property " << property_id.getString();
       task->computeProperties(const_cast<SourceWithOnDemandProperties&>(*this));
       return m_property_holder.getProperty(property_id);
     }

SEImplementation/CMakeLists.txt

@@ -37,6 +37,7 @@ find_package(BoostPython ${PYTHON_EXPLICIT_VERSION})
 find_package(Boost 1.53 REQUIRED)
 find_package(OnnxRuntime)
 find_package(Log4CPP REQUIRED)
+find_package(OpenCV REQUIRED)
 
 if (${Boost_VERSION} LESS "106700")
     message(WARNING "
@@ -104,7 +105,7 @@ elements_add_library(SEImplementation
               BoostPython PythonLibs PythonInterp
               Pyston
               ${OPT_LIBRARIES}
-            INCLUDE_DIRS ${BoostPython_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS} ${OPT_INCLUDES}
+	      INCLUDE_DIRS ${BoostPython_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS} ${OPT_INCLUDES} ${OpenCV_INCLUDE_DIRS}
             PUBLIC_HEADERS SEImplementation)
 
 #===============================================================================

SEImplementation/SEImplementation/Grouping/AssocGrouping.h

@@ -42,6 +42,7 @@ class AssocGrouping : public SourceGroupingInterface {
   std::shared_ptr<SourceGroupFactory> m_group_factory;
   std::map<unsigned int, std::unique_ptr<SourceGroupInterface>> m_source_groups;
   unsigned int m_hard_limit;
+  unsigned int m_total_sources_waiting = 0;
 
 };
 

SEImplementation/SEImplementation/Grouping/MoffatGrouping.h

@@ -60,6 +60,8 @@ class MoffatGrouping : public SourceGroupingInterface {
   size_t m_group_counter;
   std::map<unsigned int, std::shared_ptr<Group>> m_groups;
   QuadTree<std::shared_ptr<SourceInfo>> m_tree;
+
+  unsigned int m_total_sources_waiting = 0;
 };
 
 }

SEImplementation/SEImplementation/Grouping/SplitSourcesGrouping.h

@@ -42,7 +42,7 @@ class SplitSourcesGrouping : public SourceGroupingInterface {
   std::shared_ptr<SourceGroupFactory> m_group_factory;
   std::map<unsigned int, std::unique_ptr<SourceGroupInterface>> m_source_groups;
   unsigned int m_hard_limit;
-
+  unsigned int m_total_sources_waiting = 0;
 };
 
 }

SEImplementation/SEImplementation/Image/DownSampledImagePsf.h

@@ -52,9 +52,12 @@ class DownSampledImagePsf {
   std::shared_ptr<VectorImage<SourceXtractor::SeFloat>> getScaledKernel(SeFloat scale) const;
   void convolve(std::shared_ptr<WriteableImage<float>> image) const;
 
-  // For ConvolutionContext optimization
+
+#ifdef USE_DFT_CONVOLUTION_FOR_PSF
+  // // For ConvolutionContext optimization
   std::unique_ptr<DFTConvolution<SeFloat>::ConvolutionContext> prepare(const std::shared_ptr<const Image<SeFloat>>& model_ptr) const;
   void convolve(std::shared_ptr<WriteableImage<float>> image, std::unique_ptr<DFTConvolution<SeFloat>::ConvolutionContext>& context) const;
+#endif
 
 private:
   double m_down_scaling;
@@ -66,15 +69,19 @@ class DownSampledImagePsf {
 
 namespace ModelFitting {
 
+
 /**
  * Specialization of PsfTraits, as DFTConvolution has the concept of context
  */
-template<>
+#ifdef USE_DFT_CONVOLUTION_FOR_PSF
+ template<>
 struct PsfTraits<SourceXtractor::DownSampledImagePsf> {
   using context_t = typename std::unique_ptr<SourceXtractor::ImagePsf::ConvolutionContext>;
   static constexpr bool has_context = true;
 };
 
+#endif
+
 } // end of ModelFitting
 
 #endif /* _SEIMPLEMENTATION_IMAGE_DOWNSAMPLEDIMAGEPSF_H_ */

SEImplementation/SEImplementation/Image/ImagePsf.h

@@ -31,12 +31,35 @@
 #include "SEFramework/Image/ProcessedImage.h"
 #include "SEFramework/Convolution/Convolution.h"
 
+//#define USE_DFT_CONVOLUTION_FOR_PSF
+//#define USE_DIRECT_CONVOLUTION_FOR_PSF
+#define USE_OPENCV_CONVOLUTION_FOR_PSF
+
+#ifdef USE_DIRECT_CONVOLUTION_FOR_PSF
+#include "SEFramework/Convolution/DirectConvolution.h"
+#endif
+
+#ifdef USE_OPENCV_CONVOLUTION_FOR_PSF
+#include "SEFramework/Convolution/OpenCVConvolution.h"
+#endif
 
 namespace SourceXtractor {
 
+#ifdef USE_DIRECT_CONVOLUTION_FOR_PSF
+class ImagePsf: public DirectConvolution<SeFloat, PaddedImage<SeFloat, Reflect101Coordinates>> {
+private:
+  typedef DirectConvolution<SeFloat, PaddedImage<SeFloat, Reflect101Coordinates>> base_t;
+#endif
+#ifdef USE_DFT_CONVOLUTION_FOR_PSF
 class ImagePsf: public DFTConvolution<SeFloat, PaddedImage<SeFloat, Reflect101Coordinates>> {
 private:
-  typedef DFTConvolution<SeFloat, PaddedImage<SeFloat, Reflect101Coordinates>> base_t;
+typedef DFTConvolution<SeFloat, PaddedImage<SeFloat, Reflect101Coordinates>> base_t;
+#endif
+#ifdef USE_OPENCV_CONVOLUTION_FOR_PSF
+class ImagePsf: public OpenCVConvolution<SeFloat> {
+private:
+  typedef OpenCVConvolution<SeFloat> base_t;
+#endif
 
 public:
 
@@ -78,12 +101,15 @@ namespace ModelFitting {
 /**
  * Specialization of PsfTraits, as DFTConvolution has the concept of context
  */
+#ifdef USE_DFT_CONVOLUTION_FOR_PSF
 template<>
 struct PsfTraits<SourceXtractor::ImagePsf> {
   using context_t = typename std::unique_ptr<SourceXtractor::ImagePsf::ConvolutionContext>;
   static constexpr bool has_context = true;
 };
 
+#endif
+
 } // end of ModelFitting
 
 

SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h

@@ -85,6 +85,10 @@ class FlexibleModelFittingIterativeTask : public GroupTask {
     std::vector<int> iterations_per_meta;
     std::vector<SeFloat> fitting_areas_x;
     std::vector<SeFloat> fitting_areas_y;
+
+    // Cached images
+    std::vector<std::shared_ptr<VectorImage<SeFloat>>> cached_image_copies;
+    std::vector<std::shared_ptr<VectorImage<SeFloat>>> cached_weight_images;
   };
 
   struct FittingState {
@@ -111,8 +115,8 @@ class FlexibleModelFittingIterativeTask : public GroupTask {
   void updateCheckImages(SourceGroupInterface& group, double pixel_scale, FittingState& state) const;
   SeFloat computeChiSquared(SourceGroupInterface& group, SourceInterface& source, int index,
       double pixel_scale, FlexibleModelFittingParameterManager& manager, int& total_data_points, FittingState& state) const;
-  SeFloat computeChiSquaredForFrame(std::shared_ptr<const Image<SeFloat>> image,
-      std::shared_ptr<const Image<SeFloat>> model, std::shared_ptr<const Image<SeFloat>> weights, int& data_points) const;
+  SeFloat computeChiSquaredForFrame(std::shared_ptr<const VectorImage<SeFloat>> image,
+      std::shared_ptr<const VectorImage<SeFloat>> model, std::shared_ptr<const VectorImage<SeFloat>> weights, int& data_points) const;
   int fitSourcePrepareParameters(FlexibleModelFittingParameterManager& parameter_manager,
                                  ModelFitting::EngineParameterManager& engine_parameter_manager,
                                  SourceInterface& source, int index, FittingState& state) const;