find_spots: Exclude masked pixels from erosion check (dials#2790)

In order to prevent erroneous spot erasure, as seen in dials#2785.

Set masked pixels to 'true' in the dispersion threshold mask so that they
do not affect erosion and erode from masked regions as if they were
background.

Closes dials#2785.

Author: dimitrivlachos

AUTHORS

@@ -10,6 +10,7 @@ Daniel Paley
 Daniel Tchon
 David Waterman
 Derek Mendez
+Dimitrios Vlachos
 Graeme Winter
 Huw Jenkins
 Hans-Christian Stadler

newsfragments/2790.bugfix

@@ -0,0 +1 @@
+``dials.find_spots``: Prevent masked pixels from causing a "halo" of excluded pixels around them when using ``threashold.algorithm=extended_dispersion`` (the default). This previously caused strong pixels to be excluded from module edges and incorrectly classed strong pixels as background if close to masked pixels.

src/dials/algorithms/image/threshold/local.h

@@ -901,7 +901,13 @@ namespace dials { namespace algorithms {
       final_mask_ = af::versa<bool, af::c_grid<2> >(image.accessor(), false);
       global_mask_ = af::versa<bool, af::c_grid<2> >(image.accessor(), false);
       for (std::size_t i = 0; i < image.size(); ++i) {
-        if (temp[i]) {
+        if (!temp[i]) {
+          /*
+           * Masked pixels are set as non-background in order to
+           * prevent them from affecting the erosion calculation.
+           */
+          cv_mask_[i] = true;
+        } else {
           double bnd_b = gain[i] + nsig_b * gain[i] * std::sqrt(2.0 / (count[i] - 1));
           cv_mask_[i] = cv_[i] > bnd_b;
           global_mask_[i] = image[i] > threshold;
@@ -1094,7 +1100,13 @@ namespace dials { namespace algorithms {
 
           // Compute the thresholds
           dst[k] = false;
-          if (mask[k] && m >= min_count_ && x >= 0) {
+          if (!mask[k]) {
+            /*
+             * Masked pixels are set as non-background in order to
+             * prevent them from affecting the erosion calculation.
+             */
+            dst[k] = true;
+          } else if (m >= min_count_ && x >= 0) {
             double a = m * y - x * x - x * (m - 1);
             double c = x * nsig_b_ * std::sqrt(2 * (m - 1));
             dst[k] = (a > c);
@@ -1165,7 +1177,13 @@ namespace dials { namespace algorithms {
 
           // Compute the thresholds
           dst[k] = false;
-          if (mask[k] && m >= min_count_ && x >= 0) {
+          if (!mask[k]) {
+            /*
+             * Masked pixels are set as non-background in order to
+             * prevent them from affecting the erosion calculation.
+             */
+            dst[k] = true;
+          } else if (m >= min_count_ && x >= 0) {
             double a = m * y - x * x;
             double c = gain[k] * x * (m - 1 + nsig_b_ * std::sqrt(2 * (m - 1)));
             dst[k] = (a > c);