Add peak min / max pixels threshold, add a bunch of descriptions

src/calng/correction_addons/peakfinder9.py

@@ -17,6 +17,17 @@ class Peakfinder9(BaseCorrectionAddon):
         (
             UINT32_ELEMENT(schema)
             .key(f"{prefix}.windowRadius")
+            .description(
+                "The algorithm works by first deciding (in parallel, looking at each "
+                "pixel) if 'current pixel' is a candidate as start of a peak. Then, it "
+                "moves out in rings, expanding peak up until at most windowRadius. "
+                "The border of the window is used for background estimation. The other "
+                "parameters adjust rejection criteria for the 'current pixel' "
+                "(minPeakValueOverNeighbors, minSnrMaxPixel), criteria for adding "
+                "pixels to the peak (minSnrPeakPixels), and rejection criteria for the "
+                "full peak after expansion (minSnrWholePeak, minPeakPixels, "
+                "maxPeakPixels)."
+            )
             .tags("managed")
             .assignmentOptional()
             .defaultValue(2)
@@ -25,6 +36,10 @@ class Peakfinder9(BaseCorrectionAddon):
 
             UINT32_ELEMENT(schema)
             .key(f"{prefix}.maxPeaks")
+            .description(
+                "For GPU buffer reasons, we need an upper limit on how many peaks "
+                "there will be per frame. Just leave this reasonably high."
+            )
             .tags("managed")
             .assignmentOptional()
             .defaultValue(500)
@@ -33,6 +48,11 @@ class Peakfinder9(BaseCorrectionAddon):
 
             FLOAT_ELEMENT(schema)
             .key(f"{prefix}.minPeakValueOverNeighbors")
+            .description(
+                "Current pixel is rejected as start of peak if it does not exceed the "
+                "maximum pixel value on the window border by at least this amount; "
+                "note, this is in absolute terms unlike other parameters."
+            )
             .tags("managed")
             .assignmentOptional()
             .defaultValue(10)
@@ -41,6 +61,11 @@ class Peakfinder9(BaseCorrectionAddon):
 
             FLOAT_ELEMENT(schema)
             .key(f"{prefix}.minSnrMaxPixel")
+            .description(
+                "Current pixel is rejected as start of peak if it does not exceed mean "
+                "by at least mean by this amount times sigma. Both mean and sigma are "
+                "based on window border pixels."
+            )
             .tags("managed")
             .assignmentOptional()
             .defaultValue(5)
@@ -49,6 +74,12 @@ class Peakfinder9(BaseCorrectionAddon):
 
             FLOAT_ELEMENT(schema)
             .key(f"{prefix}.minSnrPeakPixels")
+            .description(
+                "If current pixel was not rejected until now, we start expanding peak "
+                "layer by layer. For a pixel to be added to the peak, it must (be in "
+                "the next layer, and) exceed the mean by at least this number times "
+                "sigma."
+            )
             .tags("managed")
             .assignmentOptional()
             .defaultValue(4)
@@ -57,14 +88,47 @@ class Peakfinder9(BaseCorrectionAddon):
 
             FLOAT_ELEMENT(schema)
             .key(f"{prefix}.minSnrWholePeak")
+            .description(
+                "After expanding peak with additional peak pixels, the whole peak is "
+                "rejected if its total mass does not exceed the mean by at least this "
+                "number times sigma."
+            )
             .tags("managed")
             .assignmentOptional()
             .defaultValue(6)
             .reconfigurable()
             .commit(),
 
+            UINT32_ELEMENT(schema)
+            .key(f"{prefix}.minPeakPixels")
+            .description(
+                "After expanding peak, it is rejected if the number of pixels is not "
+                "at least this."
+            )
+            .tags("managed")
+            .assignmentOptional()
+            .defaultValue(1)
+            .reconfigurable()
+            .commit(),
+
+            UINT32_ELEMENT(schema)
+            .key(f"{prefix}.maxPeakPixels")
+            .description(
+                "After expanding peak, it is rejected if the number of pixels exceeds "
+                "this."
+            )
+            .tags("managed")
+            .assignmentOptional()
+            .defaultValue(100)
+            .reconfigurable()
+            .commit(),
+
             FLOAT_ELEMENT(schema)
             .key(f"{prefix}.minSigma")
+            .description(
+                "When background estimate for stdev is very low, this is used to set "
+                "minimum sigma value used for the thresholds above."
+            )
             .tags("managed")
             .assignmentOptional()
             .defaultValue(5)
@@ -73,6 +137,7 @@ class Peakfinder9(BaseCorrectionAddon):
 
             UINT32_ELEMENT(schema)
             .key(f"{prefix}.blockX")
+            .description("GPU execution tuning parameter you probably want to ignore.")
             .tags("managed")
             .assignmentOptional()
             .defaultValue(1)
@@ -81,6 +146,7 @@ class Peakfinder9(BaseCorrectionAddon):
 
             UINT32_ELEMENT(schema)
             .key(f"{prefix}.blockY")
+            .description("GPU execution tuning parameter you probably want to ignore.")
             .tags("managed")
             .assignmentOptional()
             .defaultValue(1)
@@ -89,6 +155,7 @@ class Peakfinder9(BaseCorrectionAddon):
 
             UINT32_ELEMENT(schema)
             .key(f"{prefix}.blockZ")
+            .description("GPU execution tuning parameter you probably want to ignore.")
             .tags("managed")
             .assignmentOptional()
             .defaultValue(64)
@@ -214,6 +281,8 @@ class Peakfinder9(BaseCorrectionAddon):
             cupy.float32(self._config["minSnrPeakPixels"]),
             cupy.float32(self._config["minSnrWholePeak"]),
             cupy.float32(self._config["minSigma"]),
+            cupy.uint32(self._config["minPeakPixels"]),
+            cupy.uint32(self._config["maxPeakPixels"]),
             cupy.uint32(self._config["maxPeaks"]),
         )
 

src/calng/kernels/peakfinder9_gpu.cu

@@ -58,14 +58,19 @@ public:
 extern "C" __global__ void pf9(const unsigned short num_frames,
                                const unsigned short ss_dim,
                                const unsigned short fs_dim,
-                               const unsigned short window_radius,
-                               const float min_peak_over_border,
-                               const float min_snr_biggest_pixel,
-                               const float min_snr_peak_pixels,
-                               const float min_snr_whole_peak,
-                               const float min_sigma,
-                               const unsigned int max_peaks,
-                               const float* image,
+                               const unsigned short window_radius,  // radius to border: border is for background estimate and bounds peak size
+                               const float min_peak_over_border,  // for "current pixel" to be candidate, it must be this much (absolute value) over max value on border
+                               // for the next parameters, snr is based on stdev of border pixels
+                               // note: hardcoded limit says there must be at least 10 non-masked border pixels for this
+                               const float min_snr_biggest_pixel,  // "current pixel" is compared against this
+                               // we "expand" peak in rings
+                               const float min_snr_peak_pixels,  // additional peak pixels get compared against this
+                               const float min_snr_whole_peak,  // total mass of peak after expansion is compared against this
+                               const float min_sigma,  // when background has suspiciously low stdev, this sigma kicks in
+                               const unsigned min_peak_pixels,  // if less than this number of pixels are considered part of peak, reject
+                               const unsigned max_peak_pixels,  // ditto but upper bound (count must be <=)
+                               const unsigned max_peaks,  // for buffer reasons
+                               const float* image,  // input data
                                unsigned int* output_counts,
                                float* output_x,
                                float* output_y,
@@ -150,29 +155,28 @@ extern "C" __global__ void pf9(const unsigned short num_frames,
 	float peak_weighted_ss;
 	float peak_weighted_fs;
 	float peak_total_mass = 0;
+	unsigned num_peak_pixels = 1;  // will expand from center pixel, count it already
 	{
-		/* TODO: more compact form */
 		float peak_weighted_ss_nom = 0;
 		float peak_weighted_fs_nom = 0;
 		// expand peak in rings
-        {
-            const float peak_pixel_threshold = mean + min_snr_peak_pixels * sigma;
-            for (unsigned short layer=1; layer<=window_radius; ++layer) {
-                float total_mass_before = peak_total_mass;
-                masked_frame.fun_ring(ss_index, fs_index, layer, [&] (unsigned short i, unsigned short j, float val) {
-                    if (val > peak_pixel_threshold) {
-                        float val_over_mean = val - mean;
-                        peak_total_mass += val_over_mean;
-                        peak_weighted_ss_nom += val_over_mean * static_cast<float>(i);
-                        peak_weighted_fs_nom += val_over_mean * static_cast<float>(j);
-                    }
-                });
-                // in case nothing was added, stop expanding
-                if (peak_total_mass == total_mass_before) {
-                    break;
-                }
-            }
-        }
+		const float peak_pixel_threshold = mean + min_snr_peak_pixels * sigma;
+		for (unsigned short layer=1; layer<=window_radius; ++layer) {
+			unsigned num_peak_pixels_before = num_peak_pixels;
+			masked_frame.fun_ring(ss_index, fs_index, layer, [&] (unsigned short i, unsigned short j, float val) {
+				if (val > peak_pixel_threshold) {
+					float val_over_mean = val - mean;
+					peak_total_mass += val_over_mean;
+					peak_weighted_ss_nom += val_over_mean * static_cast<float>(i);
+					peak_weighted_fs_nom += val_over_mean * static_cast<float>(j);
+					++num_peak_pixels;
+				}
+			});
+			// in case nothing was added, stop expanding
+			if (num_peak_pixels == num_peak_pixels_before) {
+				break;
+			}
+		}
 		if (peak_total_mass <= mean + min_snr_whole_peak * sigma) {
 			return;
 		}