Adding Gaussian hat average mode to density map

This commit is contained in:
Matthias Koefferlein 2026-05-24 19:17:23 +02:00
parent 0648ee1bd0
commit 472a7aa7fd
3 changed files with 185 additions and 93 deletions

View File

@ -57,85 +57,7 @@
</property>
</widget>
</item>
<item row="0" column="0">
<widget class="QLabel" name="label_7">
<property name="text">
<string>Pixel size</string>
</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">
<string>µm</string>
</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">
<enum>Qt::Horizontal</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>353</width>
<height>20</height>
</size>
</property>
</spacer>
</item>
<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">
<item row="3" column="1" colspan="3">
<widget class="QComboBox" name="cb_boundary_mode">
<property name="sizePolicy">
<sizepolicy hsizetype="Fixed" vsizetype="Fixed">
@ -168,6 +90,114 @@
</item>
</widget>
</item>
<item row="3" column="0">
<widget class="QLabel" name="label_12">
<property name="text">
<string>Boundary mode</string>
</property>
</widget>
</item>
<item row="4" 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="0" column="3">
<spacer name="horizontalSpacer_2">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>353</width>
<height>20</height>
</size>
</property>
</spacer>
</item>
<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="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">
<string>µm</string>
</property>
</widget>
</item>
<item row="4" column="0">
<widget class="QLabel" name="label_9">
<property name="text">
<string>Threads</string>
</property>
</widget>
</item>
<item row="0" column="0">
<widget class="QLabel" name="label_7">
<property name="text">
<string>Pixel size</string>
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QLabel" name="label_13">
<property name="text">
<string>Averaging mode</string>
</property>
</widget>
</item>
<item row="2" column="1" colspan="3">
<widget class="QComboBox" name="cb_averaging_mode">
<property name="sizePolicy">
<sizepolicy hsizetype="Fixed" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="sizeAdjustPolicy">
<enum>QComboBox::AdjustToContentsOnFirstShow</enum>
</property>
<item>
<property name="text">
<string>Square hat</string>
</property>
</item>
<item>
<property name="text">
<string>Gaussian hat (σ = window size / 2)</string>
</property>
</item>
</widget>
</item>
</layout>
</widget>
</item>

View File

@ -66,10 +66,17 @@ 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 averaging_mode_square ("averaging-mode-square");
static const std::string averaging_mode_gaussian ("averaging-mode-gaussian");
// items in average_mode_cb
static const std::vector<std::string> averaging_modes = { averaging_mode_square, averaging_mode_gaussian };
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_averaging_mode ("density-map-averaging-mode");
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");
@ -90,6 +97,7 @@ public:
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_averaging_mode, averaging_mode_square));
options.push_back (std::make_pair (cfg_density_map_threads, "1"));
}
@ -307,6 +315,18 @@ DensityMapDialog::configure (const std::string &name, const std::string &value)
cb_boundary_mode->setCurrentIndex (mode);
return true;
} else if (name == cfg_density_map_averaging_mode) {
int mode = 0;
for (size_t i = 0; i < averaging_modes.size (); ++i) {
if (averaging_modes[i] == value) {
mode = int (i);
}
}
cb_averaging_mode->setCurrentIndex (mode);
return true;
} else if (name == cfg_density_map_threads) {
int thr = 1;
@ -390,6 +410,12 @@ DensityMapDialog::make_density_map ()
}
par.boundary_mode = boundary_modes [boundary_mode];
int averaging_mode = cb_averaging_mode->currentIndex ();
if (averaging_mode < 0 || averaging_mode >= int (averaging_modes.size ())) {
return;
}
par.averaging_mode = averaging_modes [averaging_mode];
int region_mode = region_cb->currentIndex ();
if (region_mode < 0 || region_mode >= int (region_modes.size ())) {
return;
@ -514,6 +540,7 @@ DensityMapDialog::make_density_map ()
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_averaging_mode, par.averaging_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));
@ -644,7 +671,39 @@ DensityMapDialog::compute_density_map (const DensityMapParameters &par)
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);
if (par.averaging_mode == averaging_mode_square) {
std::vector<double> weights;
int wh = nw / 2;
weights.resize (wh * 2 + 1, 1.0);
if (nw % 2 == 0) {
weights.front () = weights.back () = 0.5;
}
average_window (*img_object, par.boundary_mode, wh, weights);
} else if (par.averaging_mode == averaging_mode_gaussian) {
std::vector<double> weights;
int wh = nw * 3 / 2;
// Note: the averaging kernel (f = exp(-r^2) = exp(-(x^2+y^2)) can be decomposed
// into f = f(x)*f(y); f(x) = exp(-x^2); f(y) = exp(-y^2).
// Hence we can use a 1d-kernel here and average vertically first, then horizontally.
weights.reserve (wh * 2 + 1);
for (int i = -wh; i <= wh; ++i) {
double x = double (i) / (nw * 0.5);
weights.push_back (exp (-x * x));
}
average_window (*img_object, par.boundary_mode, wh, weights);
}
}
}
@ -658,8 +717,15 @@ inline int safe_mod (int a, int b)
}
void
DensityMapDialog::average_window (img::Object &img_object, const std::string boundary_mode, int nw)
DensityMapDialog::average_window (img::Object &img_object, const std::string boundary_mode, int wh, const std::vector<double> &weights)
{
tl_assert (weights.size () == size_t (wh * 2 + 1));
double ws = 0.0;
for (auto i = weights.begin (); i != weights.end (); ++i) {
ws += *i;
}
bool periodic = (boundary_mode == boundary_mode_periodic);
int nx = img_object.width ();
@ -673,7 +739,6 @@ DensityMapDialog::average_window (img::Object &img_object, const std::string bou
} 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++;
@ -689,13 +754,10 @@ DensityMapDialog::average_window (img::Object &img_object, const std::string bou
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;
double f = weights [dy + wh];
std::vector<double>::iterator d = vavg_data.begin () + y * nx;
if (periodic || (y + dy >= 0 && y + dy < ny)) {
@ -715,7 +777,7 @@ DensityMapDialog::average_window (img::Object &img_object, const std::string bou
// horizontal sum
outside *= nw; // because we do normalization later
outside *= ws; // because we do normalization later
// TODO: transposing the image would make things more efficient
@ -724,13 +786,10 @@ DensityMapDialog::average_window (img::Object &img_object, const std::string bou
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;
double f = weights [dx + wh];
std::vector<double>::iterator d = havg_data.begin () + x;
@ -753,7 +812,7 @@ DensityMapDialog::average_window (img::Object &img_object, const std::string bou
}
// take the average
double s = 1.0 / (double (nw) * double (nw));
double s = 1.0 / (ws * ws);
for (auto i = havg_data.begin (); i != havg_data.end (); ++i) {
*i *= s;
}

View File

@ -76,6 +76,9 @@ private:
// The boundary mode
std::string boundary_mode;
// The averaging mode
std::string averaging_mode;
// The number of threads to use
int threads;
};
@ -88,7 +91,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);
void average_window (img::Object &img_object, const std::string boundary_mode, int wh, const std::vector<double> &weights);
};