Implementing averaging window for density map

src/plugins/tools/density_map/lay_plugin/DensityMapDialog.ui

@@ -7,7 +7,7 @@
     <x>0</x>
     <y>0</y>
     <width>657</width>
-    <height>425</height>
+    <height>480</height>
    </rect>
   </property>
   <property name="windowTitle">
@@ -64,6 +64,16 @@
         </property>
        </widget>
       </item>
+      <item row="1" column="1">
+       <widget class="QLineEdit" name="le_window_size">
+        <property name="sizePolicy">
+         <sizepolicy hsizetype="Fixed" vsizetype="Fixed">
+          <horstretch>0</horstretch>
+          <verstretch>0</verstretch>
+         </sizepolicy>
+        </property>
+       </widget>
+      </item>
       <item row="0" column="2">
        <widget class="QLabel" name="label_8">
         <property name="text">
@@ -71,6 +81,33 @@
         </property>
        </widget>
       </item>
+      <item row="3" column="1">
+       <widget class="QSpinBox" name="sb_threads">
+        <property name="sizePolicy">
+         <sizepolicy hsizetype="Fixed" vsizetype="Fixed">
+          <horstretch>0</horstretch>
+          <verstretch>0</verstretch>
+         </sizepolicy>
+        </property>
+        <property name="minimum">
+         <number>1</number>
+        </property>
+       </widget>
+      </item>
+      <item row="1" column="0">
+       <widget class="QLabel" name="label_11">
+        <property name="text">
+         <string>Averaging window size</string>
+        </property>
+       </widget>
+      </item>
+      <item row="3" column="0">
+       <widget class="QLabel" name="label_9">
+        <property name="text">
+         <string>Threads</string>
+        </property>
+       </widget>
+      </item>
       <item row="0" column="3">
        <spacer name="horizontalSpacer_2">
         <property name="orientation">
@@ -84,24 +121,51 @@
         </property>
        </spacer>
       </item>
-      <item row="1" column="1">
-       <widget class="QSpinBox" name="sb_threads">
+      <item row="1" column="2" colspan="2">
+       <widget class="QLabel" name="label_10">
+        <property name="text">
+         <string>µm (a multiple of the pixel size or empty)</string>
+        </property>
+       </widget>
+      </item>
+      <item row="2" column="0">
+       <widget class="QLabel" name="label_12">
+        <property name="text">
+         <string>Boundary mode</string>
+        </property>
+       </widget>
+      </item>
+      <item row="2" column="1" colspan="3">
+       <widget class="QComboBox" name="cb_boundary_mode">
         <property name="sizePolicy">
          <sizepolicy hsizetype="Fixed" vsizetype="Fixed">
           <horstretch>0</horstretch>
           <verstretch>0</verstretch>
          </sizepolicy>
         </property>
-        <property name="minimum">
-         <number>1</number>
-        </property>
-       </widget>
-      </item>
-      <item row="1" column="0">
-       <widget class="QLabel" name="label_9">
-        <property name="text">
-         <string>Threads</string>
+        <property name="sizeAdjustPolicy">
+         <enum>QComboBox::AdjustToContents</enum>
         </property>
+        <item>
+         <property name="text">
+          <string>Periodic (repeat region seamlessly)</string>
+         </property>
+        </item>
+        <item>
+         <property name="text">
+          <string>Use zero density outside</string>
+         </property>
+        </item>
+        <item>
+         <property name="text">
+          <string>Use 100% density outside</string>
+         </property>
+        </item>
+        <item>
+         <property name="text">
+          <string>Use average density outside</string>
+         </property>
+        </item>
        </widget>
       </item>
      </layout>
@@ -473,7 +537,7 @@
       <enum>Qt::Horizontal</enum>
      </property>
      <property name="standardButtons">
-      <set>QDialogButtonBox::Apply|QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
+      <set>QDialogButtonBox::Apply|QDialogButtonBox::Close|QDialogButtonBox::Ok</set>
      </property>
     </widget>
    </item>

src/plugins/tools/density_map/lay_plugin/layDensityMapDialog.cc

@@ -58,9 +58,19 @@ static const std::vector<std::string> region_modes = {
   region_mode_by_rulers
 };
 
+static const std::string boundary_mode_periodic ("boundary-mode-periodic");
+static const std::string boundary_mode_zero ("boundary-mode-zero");
+static const std::string boundary_mode_one ("boundary-mode-one");
+static const std::string boundary_mode_average ("boundary-mode-average");
+
+//  items in boundary_mode_cb
+static const std::vector<std::string> boundary_modes = { boundary_mode_periodic, boundary_mode_zero, boundary_mode_one, boundary_mode_average };
+
 static const std::string cfg_density_map_region_mode ("density-map-region-mode");
 static const std::string cfg_density_map_layer_mode ("density-map-layer-mode");
 static const std::string cfg_density_map_pixel_size ("density-map-pixel-size");
+static const std::string cfg_density_map_window_size ("density-map-window-size");
+static const std::string cfg_density_map_boundary_mode ("density-map-boundary-mode");
 static const std::string cfg_density_map_threads ("density-map-threads");
 static const std::string cfg_density_map_source_layer ("density-map-source-layer");
 static const std::string cfg_density_map_box_layer ("density-map-box-layer");
@@ -78,6 +88,8 @@ public:
     options.push_back (std::make_pair (cfg_density_map_layer_mode, layer_mode_specific));
     options.push_back (std::make_pair (cfg_density_map_region_mode, region_mode_global_bbox));
     options.push_back (std::make_pair (cfg_density_map_pixel_size, "100"));
+    options.push_back (std::make_pair (cfg_density_map_window_size, "0"));
+    options.push_back (std::make_pair (cfg_density_map_boundary_mode, boundary_mode_periodic));
     options.push_back (std::make_pair (cfg_density_map_threads, "1"));
   }
 
@@ -269,6 +281,32 @@ DensityMapDialog::configure (const std::string &name, const std::string &value)
     }
     return true;
 
+  } else if (name == cfg_density_map_window_size) {
+
+    double ws = 0.0;
+    try {
+      tl::from_string (value, ws);
+      if (ws > 0) {
+        le_window_size->setText (tl::to_qstring (tl::to_string (ws)));
+      } else {
+        le_window_size->setText (QString ());
+      }
+    } catch (...) {
+    }
+    return true;
+
+  } else if (name == cfg_density_map_boundary_mode) {
+
+    int mode = 0;
+    for (size_t i = 0; i < boundary_modes.size (); ++i) {
+      if (boundary_modes[i] == value) {
+        mode = int (i);
+      }
+    }
+
+    cb_boundary_mode->setCurrentIndex (mode);
+    return true;
+
   } else if (name == cfg_density_map_threads) {
 
     int thr = 1;
@@ -335,9 +373,23 @@ DensityMapDialog::make_density_map ()
   tl::from_string_ext (tl::to_string (le_pixel_size->text ()), par.pixel_size);
 
   if (par.pixel_size < 1e-6) {
-    throw tl::Exception (tl::to_string (QObject::tr ("Pixel size must be positive and not zero")));
+    throw tl::Exception (tl::to_string (QObject::tr ("The pixel size must be positive and not zero")));
   }
 
+  par.window_size = 0.0;
+  if (! le_window_size->text ().simplified ().isEmpty ()) {
+    tl::from_string_ext (tl::to_string (le_window_size->text ()), par.window_size);
+    if (par.window_size < 1e-6) {
+      throw tl::Exception (tl::to_string (QObject::tr ("The window size must be positive and not zero or empty")));
+    }
+  }
+
+  int boundary_mode = cb_boundary_mode->currentIndex ();
+  if (boundary_mode < 0 || boundary_mode >= int (boundary_modes.size ())) {
+    return;
+  }
+  par.boundary_mode = boundary_modes [boundary_mode];
+
   int region_mode = region_cb->currentIndex ();
   if (region_mode < 0 || region_mode >= int (region_modes.size ())) {
     return;
@@ -461,8 +513,10 @@ DensityMapDialog::make_density_map ()
   //  Commit the parameters
   dispatcher ()->config_set (cfg_density_map_layer_mode, lm);
   dispatcher ()->config_set (cfg_density_map_region_mode, rm);
+  dispatcher ()->config_set (cfg_density_map_boundary_mode, par.boundary_mode);
   dispatcher ()->config_set (cfg_density_map_threads, tl::to_string (par.threads));
   dispatcher ()->config_set (cfg_density_map_pixel_size, tl::to_string (par.pixel_size));
+  dispatcher ()->config_set (cfg_density_map_window_size, tl::to_string (par.window_size));
 
   if (lm == layer_mode_specific) {
     dispatcher ()->config_set (cfg_density_map_source_layer, source_layer.to_string ());
@@ -584,6 +638,127 @@ DensityMapDialog::compute_density_map (const DensityMapParameters &par)
   //  Execute the tiling processor
 
   tp.execute (tl::to_string (tr ("Computing density map")));
+
+  //  Do the averaging if requested
+
+  unsigned int nw = (unsigned int) (std::max (0.0, std::min (1000.0, floor (par.window_size / par.pixel_size + 0.5))));
+
+  if (nw > 1) {
+    average_window (*img_object, par.boundary_mode, nw);
+  }
+}
+
+inline int safe_mod (int a, int b)
+{
+  if (a < 0) {
+    return b - (-a % b);
+  } else {
+    return a % b;
+  }
+}
+
+void
+DensityMapDialog::average_window (img::Object &img_object, const std::string boundary_mode, int nw)
+{
+  bool periodic = (boundary_mode == boundary_mode_periodic);
+
+  int nx = img_object.width ();
+  int ny = img_object.height ();
+
+  //  compute the outside value if not periodic
+  double outside = 0.0;
+  if (boundary_mode == boundary_mode_one) {
+
+    outside = 1.0;
+
+  } else if (boundary_mode == boundary_mode_average) {
+
+    double outside = 0.0;
+    const float *d = img_object.float_data ();
+    for (size_t i = size_t (nx) * size_t (ny); i > 0; --i) {
+      outside += *d++;
+    }
+    outside /= double (nx) * double (ny);
+
+  }
+
+  std::vector<double> vavg_data;
+  vavg_data.resize (size_t (nx) * size_t (ny), 0.0);
+
+  //  vertical sum
+
+  for (int y = 0; y < ny; ++y) {
+
+    int wh = nw / 2;
+    double fb = (nw % 2 == 0) ? 0.5 : 1.0;
+
+    for (int dy = -wh; dy <= wh; ++dy) {
+
+      //  top and bottom row count half in case of even nw
+      double f = (dy == -wh || dy == wh) ? fb : 1.0;
+
+      std::vector<double>::iterator d = vavg_data.begin () + y * nx;
+      if (periodic || (y + dy >= 0 && y + dy < ny)) {
+        const float *s = img_object.float_data () + safe_mod (y + dy, ny) * nx;
+        for (int ix = 0; ix < nx; ++ix) {
+          *d++ += f * *s++;
+        }
+      } else {
+        for (int ix = 0; ix < nx; ++ix) {
+          *d++ += f * outside;
+        }
+      }
+
+    }
+
+  }
+
+  //  horizontal sum
+
+  outside *= nw;   //  because we do normalization later
+
+  //  TODO: transposing the image would make things more efficient
+
+  std::vector<double> havg_data;
+  havg_data.resize (size_t (nx) * size_t (ny), 0.0);
+
+  for (int x = 0; x < nx; ++x) {
+
+    int wh = nw / 2;
+    double fb = (nw % 2 == 0) ? 0.5 : 1.0;
+
+    for (int dx = -wh; dx <= wh; ++dx) {
+
+      //  top and bottom row count half in case of even nw
+      double f = (dx == -wh || dx == wh) ? fb : 1.0;
+
+      std::vector<double>::iterator d = havg_data.begin () + x;
+
+      if (periodic || (x + dx >= 0 && x + dx < nx)) {
+        std::vector<double>::const_iterator s = vavg_data.begin () + safe_mod (x + dx, nx);
+        for (int iy = 0; iy < ny; ++iy) {
+          *d += f * *s;
+          d += nx;
+          s += nx;
+        }
+      } else {
+        for (int iy = 0; iy < ny; ++iy) {
+          *d += f * outside;
+          d += nx;
+        }
+      }
+
+    }
+
+  }
+
+  //  take the average
+  double s = 1.0 / (double (nw) * double (nw));
+  for (auto i = havg_data.begin (); i != havg_data.end (); ++i) {
+    *i *= s;
+  }
+
+  img_object.set_data (nx, ny, havg_data);
 }
 
 }

src/plugins/tools/density_map/lay_plugin/layDensityMapDialog.h

@@ -30,6 +30,11 @@
 #include "layBrowser.h"
 #include "layMarker.h"
 
+namespace img
+{
+  class Object;
+}
+
 namespace lay
 {
 
@@ -65,6 +70,12 @@ private:
     //  The pixel size
     double pixel_size;
 
+    //  The window size or zero for "no window"
+    double window_size;
+
+    //  The boundary mode
+    std::string boundary_mode;
+
     //  The number of threads to use
     int threads;
   };
@@ -77,6 +88,7 @@ private:
 
   void make_density_map ();
   void compute_density_map (const DensityMapParameters &par);
+  void average_window (img::Object &img_object, const std::string boundary_mode, int nw);
 
 };