Saturated pixels not being masked

AGENTS.md

@@ -0,0 +1,52 @@
+# Repository instructions
+
+## Building and testing
+
+This project is CMake-based but uses the Elements build system. Never run CMake
+commands directly; use the project's Makefile through Elements instead.
+
+Run the following commands from the repository root as needed:
+
+- `make` — build the project.
+- `make install` — install the build products.
+- `make tests` — build and run the tests.
+
+Always run `make install` before integration runs such as executing `sourcextractor++`
+in `tests/sim12`, and use the installed executable and environment.
+
+## Coding guidelines
+
+- Follow `.editorconfig`: UTF-8, LF line endings, spaces, two-space indentation,
+  a maximum line length of 120 characters, no trailing whitespace, and a final
+  newline. Python uses four-space indentation; Makefiles use tabs.
+- Match the surrounding code when conventions vary. Keep changes focused and
+  avoid reformatting unrelated code.
+- In C++, use `PascalCase` for classes, `camelCase` for methods and functions,
+  `snake_case` for local variables and parameters, and `m_snake_case` for data
+  members. Match existing public API names.
+- Put opening braces on the same line as class, function, and control-flow
+  declarations. Use spaces after control-flow keywords (`if (...)`, `for (...)`).
+  Keep namespace contents and class access labels unindented, as in nearby code.
+- Always use braces `{}` around C++ `if`, `else if`, `else`, and `for` bodies,
+  even for a single statement. Do not use unbraced conditionals or `for` loops,
+  even when nearby code does so.
+- Follow the module layout: public headers in `<Module>/<Module>/`, implementation
+  files in `<Module>/src/lib/`, and tests in `<Module>/tests/src/`. Use `.h` and
+  `.cpp` files named after the class where applicable.
+- Use header guards and module-qualified project includes such as
+  `"SEFramework/Task/SourceTask.h"`. Follow nearby include grouping and guard names.
+  Keep code in its existing namespace (`SourceXtractor`, `ModelFitting`, etc.).
+- Use constructor initializer lists for members, `explicit` for single-argument
+  constructors where implicit conversion is not intended, and `override` for
+  overridden virtual methods. Mark read-only member functions `const` and use
+  const references for read-only access to nontrivial objects.
+- Follow existing ownership patterns with `std::unique_ptr` and `std::shared_ptr`;
+  use RAII for resource management. Use `auto` where the type is clear from context.
+- Preserve numerical types and precision choices in the surrounding code,
+  including the aliases in `SEUtils/Types.h` where used.
+- Preserve existing license headers. Use Doxygen comments for public C++ API
+  documentation and explain non-obvious algorithms, assumptions, and units.
+- Follow the existing Boost.Test style for C++ tests: `*_test.cpp` files,
+  `BOOST_AUTO_TEST_SUITE`, and test cases or fixtures matching nearby tests.
+- In Python, use `snake_case` functions and variables, `PascalCase` classes, and
+  docstrings consistent with the surrounding module.

SEImplementation/CMakeLists.txt

@@ -187,6 +187,10 @@ elements_add_unit_test(JacobianSource_test tests/src/Plugin/Jacobian/JacobianSou
 		             LINK_LIBRARIES SEImplementation
 		             TYPE Boost)
 if (NOT WITHOUT_MODELFITTING)
+elements_add_unit_test(FlexibleModelFittingWeight_test
+                     tests/src/Plugin/FlexibleModelFitting/FlexibleModelFittingWeight_test.cpp
+                     LINK_LIBRARIES SEImplementation
+                     TYPE Boost)
 elements_add_unit_test(MoffatModelFitting_test tests/src/Plugin/MoffatModelFitting/MoffatModelFitting_test.cpp
                      LINK_LIBRARIES SEImplementation
                      TYPE Boost)

SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h

@@ -73,6 +73,7 @@ class FlexibleModelFittingIterativeTask : public GroupTask {
   void computeProperties(SourceGroupInterface& group) const override;
 
 private:
+  friend struct FlexibleModelFittingWeightFixture;
   struct SourceState {
     std::unordered_map<int, double> parameters_initial_values;
     std::unordered_map<int, double> parameters_values;

SEImplementation/SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingTask.h

@@ -54,6 +54,7 @@ class FlexibleModelFittingTask : public GroupTask {
   void computeProperties(SourceGroupInterface& group) const override;
 
 private:
+  friend struct FlexibleModelFittingWeightFixture;
 
   bool isFrameValid(SourceGroupInterface& group, int frame_index) const;
 

SEImplementation/SEImplementation/Plugin/MoffatModelFitting/MoffatModelFittingTask.h

@@ -25,6 +25,7 @@
 #define _SEIMPLEMENTATION_PLUGIN_MOFFATMODELFITTING_MOFFATMODELFITTINGTASK_H_
 
 #include "SEFramework/Task/SourceTask.h"
+#include "SEFramework/Image/VectorImage.h"
 
 #include "SEImplementation/Image/ImagePsf.h"
 
@@ -41,6 +42,9 @@ class MoffatModelFittingTask : public SourceTask {
   void computeProperties(SourceInterface& source) const override;
 
 private:
+  friend struct MoffatModelFittingWeightFixture;
+
+  std::shared_ptr<VectorImage<SeFloat>> createWeightImage(SourceInterface& source) const;
 
   std::string m_least_squares_engine;
   unsigned int m_max_iterations;

SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.cpp

@@ -308,7 +308,7 @@ std::shared_ptr<VectorImage<SeFloat>> FlexibleModelFittingIterativeTask::createW
   auto ref_coordinates = source.getProperty<ReferenceCoordinates>().getCoordinateSystem();
 
   auto frame_image = frame_images.getLockedImage(LayerSubtractedImage);
-  auto frame_image_thresholded = frame_images.getLockedImage(LayerThresholdedImage);
+  auto original_image = frame_images.getLockedImage(LayerOriginalImage);
   auto variance_map = frame_images.getLockedImage(LayerVarianceMap);
 
   const auto& frame_info = source.getProperty<MeasurementFrameInfo>(frame_index);
@@ -356,7 +356,9 @@ std::shared_ptr<VectorImage<SeFloat>> FlexibleModelFittingIterativeTask::createW
       auto dx = x + rect.getTopLeft().m_x - ellipse.m_x;
       auto dy = y + rect.getTopLeft().m_y - ellipse.m_y;
 
-      if (saturation > 0 && pixel_val > saturation) {
+      // Saturation is defined before background subtraction, in the original image units.
+      if (saturation > 0 && original_image->getValue(
+          rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y) >= saturation) {
         weight->at(x, y) = 0;
       } else if (gain > 0.0 && pixel_val > 0.0) {
         weight->at(x, y) = sqrt(1.0 / (back_var + pixel_val / gain));

SEImplementation/src/lib/Plugin/FlexibleModelFitting/FlexibleModelFittingTask.cpp

@@ -95,7 +95,7 @@ std::shared_ptr<VectorImage<SeFloat>> FlexibleModelFittingTask::createWeightImag
   const auto& frame_images = group.begin()->getProperty<MeasurementFrameImages>(frame_index);
 
   auto frame_image = frame_images.getLockedImage(LayerSubtractedImage);
-  auto frame_image_thresholded = frame_images.getLockedImage(LayerThresholdedImage);
+  auto original_image = frame_images.getLockedImage(LayerOriginalImage);
   auto variance_map = frame_images.getLockedImage(LayerVarianceMap);
 
   const auto& frame_info = group.begin()->getProperty<MeasurementFrameInfo>(frame_index);
@@ -109,7 +109,9 @@ std::shared_ptr<VectorImage<SeFloat>> FlexibleModelFittingTask::createWeightImag
     for (int x = 0; x < rect.getWidth(); x++) {
       auto back_var = variance_map->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y);
       auto pixel_val = frame_image->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y);
-      if (saturation > 0 && pixel_val > saturation) {
+      // Saturation is defined before background subtraction, in the original image units.
+      if (saturation > 0 && original_image->getValue(
+          rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y) >= saturation) {
         weight->at(x, y) = 0;
       }
       else if (gain > 0.0 && pixel_val > 0.0) {

SEImplementation/src/lib/Plugin/MoffatModelFitting/MoffatModelFittingTask.cpp

@@ -63,6 +63,7 @@
 #include "ModelFitting/Engine/DataVsModelResiduals.h"
 
 #include "SEImplementation/Plugin/DetectionFrameSourceStamp/DetectionFrameSourceStamp.h"
+#include "SEImplementation/Plugin/DetectionFrameImages/DetectionFrameImages.h"
 
 #include "SEImplementation/Plugin/PixelCentroid/PixelCentroid.h"
 #include "SEImplementation/Property/PixelCoordinateList.h"
@@ -150,10 +151,60 @@ struct SourceModel {
 }
 
 
+std::shared_ptr<VectorImage<SeFloat>> MoffatModelFittingTask::createWeightImage(SourceInterface& source) const {
+  const auto& stamp = source.getProperty<DetectionFrameSourceStamp>();
+  const auto& source_stamp = stamp.getStamp();
+  const auto& variance_stamp = stamp.getVarianceStamp();
+  const auto& thresholded_stamp = stamp.getThresholdedStamp();
+  auto stamp_top_left = stamp.getTopLeft();
+  const auto& pixel_coordinates = source.getProperty<PixelCoordinateList>().getCoordinateList();
+
+  auto weight = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
+  std::fill(weight->getData().begin(), weight->getData().end(), 1);
+
+  for (int y = 0; y < source_stamp.getHeight(); y++) {
+    for (int x = 0; x < source_stamp.getWidth(); x++) {
+      weight->at(x, y) = (thresholded_stamp.getValue(x, y) >= 0) ? 0 : 1;
+    }
+  }
+
+  for (auto pixel : pixel_coordinates) {
+    pixel -= stamp_top_left;
+    weight->at(pixel.m_x, pixel.m_y) = 1;
+  }
+
+  const auto& detection_frame_info = source.getProperty<DetectionFrameInfo>();
+  SeFloat gain = detection_frame_info.getGain();
+  SeFloat saturation = detection_frame_info.getSaturation();
+
+  std::shared_ptr<ImageChunk<SeFloat>> original_stamp;
+  if (saturation > 0) {
+    const auto& frame_images = source.getProperty<DetectionFrameImages>();
+    original_stamp = frame_images.getImageChunk(LayerOriginalImage,
+        stamp_top_left.m_x, stamp_top_left.m_y, source_stamp.getWidth(), source_stamp.getHeight());
+  }
+
+  for (int y = 0; y < source_stamp.getHeight(); y++) {
+    for (int x = 0; x < source_stamp.getWidth(); x++) {
+      auto back_var = variance_stamp.getValue(x, y);
+      // The threshold is defined in original-image units, before background subtraction.
+      if (saturation > 0 && original_stamp->getValue(x, y) >= saturation) {
+        weight->at(x, y) = 0;
+      } else if (weight->at(x, y)>0) {
+        if (gain > 0.0) {
+          weight->at(x, y) = sqrt(1.0 / (back_var + source_stamp.getValue(x, y) / gain));
+        } else {
+          weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
+        }
+      }
+    }
+  }
+
+  return weight;
+}
+
 void MoffatModelFittingTask::computeProperties(SourceInterface& source) const {
   auto& source_stamp = source.getProperty<DetectionFrameSourceStamp>().getStamp();
-  auto& variance_stamp = source.getProperty<DetectionFrameSourceStamp>().getVarianceStamp();
-  auto& thresholded_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdedStamp();
   auto& threshold_map_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdMapStamp();
   PixelCoordinate stamp_top_left = source.getProperty<DetectionFrameSourceStamp>().getTopLeft();
 
@@ -207,38 +258,7 @@ void MoffatModelFittingTask::computeProperties(SourceInterface& source) const {
     std::move(extended_models)
   };
 
-  auto weight = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
-  std::fill(weight->getData().begin(), weight->getData().end(), 1);
-
-  for (int y=0; y < source_stamp.getHeight(); y++) {
-    for (int x=0; x < source_stamp.getWidth(); x++) {
-      weight->at(x, y) = (thresholded_stamp.getValue(x, y) >= 0) ? 0 : 1;
-    }
-  }
-
-  for (auto pixel : pixel_coordinates) {
-    pixel -= stamp_top_left;
-    weight->at(pixel.m_x, pixel.m_y) = 1;
-  }
-
-  const auto& detection_frame_info = source.getProperty<DetectionFrameInfo>();
-  SeFloat gain = detection_frame_info.getGain();
-  SeFloat saturation = detection_frame_info.getSaturation();
-
-  for (int y=0; y < source_stamp.getHeight(); y++) {
-    for (int x=0; x < source_stamp.getWidth(); x++) {
-      auto back_var = variance_stamp.getValue(x, y);
-      if (saturation > 0 && source_stamp.getValue(x, y) >= saturation) {
-        weight->at(x, y) = 0;
-      } else if (weight->at(x, y)>0) {
-        if (gain > 0.0) {
-          weight->at(x, y) = sqrt(1.0 / (back_var + source_stamp.getValue(x, y) / gain));
-        } else {
-          weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
-        }
-      }
-    }
-  }
+  auto weight = createWeightImage(source);
 
   // FIXME we should be able to use the source_stamp Image interface directly
   auto image = VectorImage<SeFloat>::create(source_stamp);

SEImplementation/tests/src/Plugin/FlexibleModelFitting/FlexibleModelFittingWeight_test.cpp

@@ -0,0 +1,132 @@
+/** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université
+ *
+ * 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
+ */
+#include <boost/test/unit_test.hpp>
+
+#include <cmath>
+
+#include "SEFramework/Source/SimpleSource.h"
+#include "SEFramework/Source/SimpleSourceGroup.h"
+#include "SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingIterativeTask.h"
+#include "SEImplementation/Plugin/FlexibleModelFitting/FlexibleModelFittingTask.h"
+#include "SEImplementation/Plugin/MeasurementFrameCoordinates/MeasurementFrameCoordinates.h"
+#include "SEImplementation/Plugin/MeasurementFrameGroupRectangle/MeasurementFrameGroupRectangle.h"
+#include "SEImplementation/Plugin/MeasurementFrameImages/MeasurementFrameImages.h"
+#include "SEImplementation/Plugin/MeasurementFrameInfo/MeasurementFrameInfo.h"
+#include "SEImplementation/Plugin/MeasurementFrameRectangle/MeasurementFrameRectangle.h"
+#include "SEImplementation/Plugin/ReferenceCoordinates/ReferenceCoordinates.h"
+#include "SEImplementation/Plugin/WorldCentroid/WorldCentroid.h"
+
+namespace SourceXtractor {
+namespace {
+
+class WeightTestCoordinates : public CoordinateSystem {
+public:
+  WorldCoordinate imageToWorld(ImageCoordinate coord) const override {
+    return WorldCoordinate(coord.m_x, coord.m_y);
+  }
+
+  ImageCoordinate worldToImage(WorldCoordinate) const override {
+    return ImageCoordinate(2, 2);
+  }
+};
+
+}
+
+struct FlexibleModelFittingWeightFixture {
+  void checkWeights(bool iterative, SeFloat saturation, SeFloat scale, bool varying_background, SeFloat gain, SeFloat background_level = 30) {
+    const int size = 5;
+    auto image = VectorImage<SeFloat>::create(size, size);
+    auto background = VectorImage<SeFloat>::create(size, size);
+    auto variance = VectorImage<SeFloat>::create(size, size);
+    variance->fillValue(4 * scale * scale);
+
+    // Below, exactly at, and above saturation; also exercise nonpositive signal.
+    const SeFloat pixels[] = {99, 100, 101, 20, 10};
+    for (int y = 0; y < size; ++y) {
+      for (int x = 0; x < size; ++x) {
+        image->at(x, y) = pixels[x] * scale;
+        background->at(x, y) = (background_level + (varying_background ? y * 5 : 0)) * scale;
+      }
+    }
+
+    auto coordinates = std::make_shared<WeightTestCoordinates>();
+    auto frame = std::make_shared<MeasurementImageFrame>(
+        image, variance, 1e6, coordinates, gain, saturation, 0);
+    frame->setBackgroundLevel(background, 2 * scale);
+    auto source = std::make_unique<SimpleSource>();
+    source->setIndexedProperty<MeasurementFrameImages>(0, frame, size, size);
+    source->setIndexedProperty<MeasurementFrameInfo>(0, size, size, gain, saturation, 1e6, 2 * scale);
+    source->setIndexedProperty<MeasurementFrameCoordinates>(0, coordinates);
+    source->setProperty<ReferenceCoordinates>(coordinates);
+    source->setProperty<WorldCentroid>(2, 2);
+    source->setIndexedProperty<MeasurementFrameRectangle>(0, ImageCoordinate(1, 1), ImageCoordinate(3, 3));
+
+    std::shared_ptr<VectorImage<SeFloat>> weights;
+    if (iterative) {
+      FlexibleModelFittingIterativeTask task("", 1, 1, {}, {}, {}, {});
+      weights = task.createWeightImage(*source, 0);
+    } else {
+      SimpleSourceGroup group;
+      group.addSource(std::move(source));
+      group.setIndexedProperty<MeasurementFrameGroupRectangle>(0, PixelCoordinate(0, 0), PixelCoordinate(4, 4));
+      FlexibleModelFittingTask task("", 1, 1, {}, {}, {});
+      weights = task.createWeightImage(group, 0);
+    }
+
+    BOOST_REQUIRE_EQUAL(weights->getWidth(), size);
+    BOOST_REQUIRE_EQUAL(weights->getHeight(), size);
+    for (int y = 0; y < size; ++y) {
+      for (int x = 0; x < size; ++x) {
+        if (saturation > 0 && x >= 1 && x <= 2) {
+          BOOST_CHECK_EQUAL(weights->getValue(x, y), 0);
+        } else {
+          // Shot noise continues to use the background-subtracted signal.
+          double signal = image->getValue(x, y) - background->getValue(x, y);
+          double expected_variance = 4 * scale * scale;
+          if (gain > 0 && signal > 0) {
+            expected_variance += signal / gain;
+          }
+          BOOST_CHECK_CLOSE(weights->getValue(x, y), 1 / std::sqrt(expected_variance), 0.001);
+        }
+      }
+    }
+  }
+};
+
+BOOST_FIXTURE_TEST_SUITE(FlexibleModelFittingWeight_test, FlexibleModelFittingWeightFixture)
+
+BOOST_AUTO_TEST_CASE(saturation_before_background_subtraction) {
+  for (bool iterative : {true, false}) {
+    checkWeights(iterative, 100, 1, false, 2);
+    checkWeights(iterative, 100, 1, false, 2, 0);
+    checkWeights(iterative, 100, 1, false, 2, -30);
+    checkWeights(iterative, 100, 1, true, 2);
+    checkWeights(iterative, 100, 1, true, 0);
+    checkWeights(iterative, 200, 2, true, 1);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(disabled_saturation) {
+  for (bool iterative : {true, false}) {
+    checkWeights(iterative, 0, 1, true, 2);
+    checkWeights(iterative, -1, 1, true, 2);
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+}