Skip to content

I
imagej-elphel
Commit 2eb14ba0 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

CLAUDE: Cluster energy normalization in recurrent buffer conditioning 

Per Andrey's continuity requirement: a sub-pixel target shift switching
between 1 and 2 surviving cells must not change the total injected energy
(filterConv5dROI does not conserve energy - it zeroes weaker components).

condition() now labels 4D-connected survivor clusters (3^4-1 connectivity,
BFS) and applies three regimes:
1. Single survivor: full 3x3x3x3 spread with fading peripherals (as before).
2. Compact cluster (2..curt_recur_max_clust near-equal survivors): spread
   limited along resolved dimensions; cluster total output = cluster max M,
   shared proportionally (survivor i contributes M*m_i/sum(m_j)) - exactly
   continuous across the 1<->2 survivor transition.
3. Cluster larger than curt_recur_max_clust (default 5): weak signal that
   passed the filter linearly - identity, no spreading/renormalization.

New parameter curt_recur_max_clust plumbed through dialog/properties/copy;
-RC<n> added to the -RECUR image titles.
Co-Authored-By: 's avatarClaude Fable 5 <noreply@anthropic.com>
parent 079d4237
Hide whitespace changes
Inline Side-by-side
Showing with 102 additions and 21 deletions
src/main/java/com/elphel/imagej/cuas/rt/CuasDetectRT.java
View file @ 2eb14ba0
......@@ -178,6 +178,7 @@ public class CuasDetectRT {
double curt_recur_leak = clt_parameters.imp.curt_recur_leak; // By Claude on 06/09/2026
double curt_recur_wthresh =clt_parameters.imp.curt_recur_wthresh; // By Claude on 06/09/2026
double curt_recur_spread = clt_parameters.imp.curt_recur_spread; // By Claude on 06/09/2026
int curt_recur_max_clust = clt_parameters.imp.curt_recur_max_clust; // By Claude on 06/09/2026
double[] curt_temp_weights = clt_parameters.imp.curt_temp_weights.clone();
......@@ -336,7 +337,8 @@ public class CuasDetectRT {
curt_vel_radius, // int velocity_radius, // By Claude on 06/09/2026
curt_vel_decimate, // int velocity_decimate, // By Claude on 06/09/2026
curt_recur_period, // int period, // By Claude on 06/09/2026
curt_recur_spread); // double spread) // By Claude on 06/09/2026
curt_recur_spread, // double spread, // By Claude on 06/09/2026
curt_recur_max_clust); // int max_clust) // By Claude on 06/09/2026
int max_recur = num_5d_lev0 / curt_recur_period + 1; // By Claude on 06/09/2026
recur_states = new double[max_recur][][][]; // By Claude on 06/09/2026
recur_ts = new String[max_recur]; // By Claude on 06/09/2026
......@@ -408,7 +410,8 @@ public class CuasDetectRT {
System.arraycopy(recur_ts, 0, recur_ts_show, 0, num_recur); // By Claude on 06/09/2026
String title_recur = title_conv5d+ // By Claude on 06/09/2026
"-RP"+curt_recur_period+"-RW"+d2s(curt_recur_w)+"-RL"+d2s(curt_recur_leak)+ // By Claude on 06/09/2026
"-RWT"+d2s(curt_recur_wthresh)+"-RS"+d2s(curt_recur_spread)+SUFFIX_RECUR; // By Claude on 06/09/2026
"-RWT"+d2s(curt_recur_wthresh)+"-RS"+d2s(curt_recur_spread)+ // By Claude on 06/09/2026
"-RC"+curt_recur_max_clust+SUFFIX_RECUR; // By Claude on 06/09/2026
ImagePlus imp_recur_rect = cuasRTUtils.showConvKernel5d( // By Claude on 06/09/2026
recur_show, // By Claude on 06/09/2026
curt_save_select, // By Claude on 06/09/2026
......
src/main/java/com/elphel/imagej/cuas/rt/CuasRecurrentBuffer.java
View file @ 2eb14ba0
......@@ -16,7 +16,13 @@ package com.elphel.imagej.cuas.rt;
* integer (1 px per minimum velocity step) when period == velocity_decimate.
* Points predicted outside the ROI are dropped (target left the watched area).
* - conditioned(measurement): adaptive 4D spreading of the filterConv5dROI survivors,
* see condition(). Each survivor contributes its value with total spread weight 1.
* see condition(). Survivors are grouped into 4D-connected clusters; each compact
* cluster (a single target) injects total energy equal to the cluster maximum,
* independent of the survivor count - this compensates the filter, which does not
* conserve energy (it zeroes weaker components), and keeps the output continuous
* when a sub-pixel position change switches between 1 and 2 surviving cells.
* Large clusters (weak signal passing the filter linearly) bypass the
* renormalization and spreading - identity.
* - w_eff: fixed `w` when wthresh <= 0; otherwise per-pixel adaptive
* w_eff = w * tanh(strength / wthresh), where strength is the strongest conditioned
* measurement at that pixel. Dark frames then give w_eff ~ 0 (pure prediction),
......@@ -36,6 +42,7 @@ public class CuasRecurrentBuffer {
private final int velocity_decimate;
private final int period; // time units between updates
private final double spread; // weight of an orthogonal 4D neighbor when spreading (diagonals get products)
private final int max_clust; // maximal 4D cluster size treated as a compact target; larger - linear mode
private double [][][] buffer; // [roi_width*roi_height][num_sub][vel_dim*vel_dim]
private int num_updates = 0;
......@@ -46,7 +53,8 @@ public class CuasRecurrentBuffer {
int velocity_radius,
int velocity_decimate,
int period,
double spread) {
double spread,
int max_clust) {
this.roi_width = roi_width;
this.roi_height = roi_height;
this.num_sub = num_sub;
......@@ -55,6 +63,7 @@ public class CuasRecurrentBuffer {
this.velocity_decimate = velocity_decimate;
this.period = period;
this.spread = spread;
this.max_clust = max_clust;
this.buffer = new double [roi_width*roi_height][num_sub][vel_dim*vel_dim];
}
......@@ -144,26 +153,86 @@ public class CuasRecurrentBuffer {
/**
* Adaptive 4D spreading of the filtered measurement ("conditioning").
* Each survivor (value > 0 after filterConv5dROI) is spread into its 4D
* (px, py, vx, vy) neighbourhood, except along "resolved" dimensions:
* a dimension is resolved when another survivor exists among the 3^4-1
* immediate 4D neighbors offset along that dimension - the survivors already
* span it, so expanding would only blur real resolution.
* - 1 isolated survivor: full 3x3x3x3 hypercube (orthogonal neighbors weight
* `spread`, diagonals `spread^2` etc.).
* - 2 survivors adjacent in velocity: that velocity axis not expanded (2x3x3x3).
* Total contributed weight per survivor is 1 regardless of geometry (weights are
* normalized over the in-bounds part of the adaptive kernel).
* Survivors (value > 0 after filterConv5dROI) are grouped into 4D-connected
* clusters (3^4-1 connectivity in (px, py, vx, vy)). Three regimes:
* - 1 isolated survivor: full 3x3x3x3 spread with fading peripherals
* (orthogonal neighbors weight `spread`, diagonals products of it).
* - 2..max_clust near-equal survivors (one target straddling cell boundaries;
* near-equality is guaranteed by the filter): each survivor's spread is
* limited along the "resolved" dimensions spanned by its neighbors
* (e.g. 2x3x3x3 for a velocity-adjacent pair - no fading added there).
* - more than max_clust survivors in one cluster: weak signal that passed the
* filter linearly - identity pass-through, no spreading, no renormalization.
* Energy normalization (compact clusters): the filter does not conserve energy -
* it zeroes weaker components - so a sub-pixel target shift switching between
* 1 and 2 surviving cells must not change the injected energy. Each cluster
* contributes a total equal to its maximum M, shared between survivors
* proportionally to their values: survivor i contributes M*m_i/sum(m_j).
* (1 survivor -> total m1; a near-equal pair -> ~M/2 each, total M = m1 -
* continuous across the transition.)
*/
private double [][][] condition(double [][][] meas) {
double [][][] out = new double [meas.length][num_sub][vel_dim*vel_dim];
for (int p = 0; p < meas.length; p++) {
int px = p % roi_width;
int py = p / roi_width;
for (int s = 0; s < num_sub; s++) {
final int nvel = vel_dim * vel_dim;
final int [] cell_stack = new int [meas.length * nvel];
final int [] label = new int [meas.length * nvel];
for (int s = 0; s < num_sub; s++) {
// Pass 1: label 4D-connected clusters of survivors, collect per-cluster size/max/sum
java.util.Arrays.fill(label, -1);
java.util.ArrayList<double[]> clusters = new java.util.ArrayList<double[]>(); // {size, max, sum}
for (int cell0 = 0; cell0 < label.length; cell0++) {
if ((label[cell0] >= 0) || (meas[cell0 / nvel][s][cell0 % nvel] <= 0)) continue;
int ncl = clusters.size();
double cl_size = 0, cl_max = 0, cl_sum = 0;
int sp = 0;
cell_stack[sp++] = cell0;
label[cell0] = ncl;
while (sp > 0) {
int cell = cell_stack[--sp];
int p = cell / nvel, v = cell % nvel;
double val = meas[p][s][v];
cl_size += 1;
cl_sum += val;
if (val > cl_max) cl_max = val;
int px = p % roi_width, py = p / roi_width;
int vx = v % vel_dim, vy = v / vel_dim;
for (int dpy = -1; dpy <= 1; dpy++) {
int ny = py + dpy;
if ((ny < 0) || (ny >= roi_height)) continue;
for (int dpx = -1; dpx <= 1; dpx++) {
int nx = px + dpx;
if ((nx < 0) || (nx >= roi_width)) continue;
int np = ny * roi_width + nx;
for (int dvy = -1; dvy <= 1; dvy++) {
int nvy = vy + dvy;
if ((nvy < 0) || (nvy >= vel_dim)) continue;
for (int dvx = -1; dvx <= 1; dvx++) {
int nvx = vx + dvx;
if ((nvx < 0) || (nvx >= vel_dim)) continue;
int ncell = np * nvel + nvy * vel_dim + nvx;
if ((label[ncell] < 0) && (meas[np][s][nvy * vel_dim + nvx] > 0)) {
label[ncell] = ncl;
cell_stack[sp++] = ncell;
}
}
}
}
}
}
clusters.add(new double[] {cl_size, cl_max, cl_sum});
}
// Pass 2: spread each survivor according to its cluster regime
for (int p = 0; p < meas.length; p++) {
int px = p % roi_width;
int py = p / roi_width;
for (int v = 0; v < meas[p][s].length; v++) {
double m = meas[p][s][v];
if (m <= 0) continue; // survivors only
double [] cl = clusters.get(label[p * nvel + v]);
if (cl[0] > max_clust) { // weak signal, linear mode - identity, no renormalization
out[p][s][v] += m;
continue;
}
int vx = v % vel_dim;
int vy = v / vel_dim;
// Find resolved dimensions: scan the 3^4-1 immediate 4D neighbors for other survivors
......@@ -196,7 +265,7 @@ public class CuasRecurrentBuffer {
int ry0 = res_y ? 0 : -1, ry1 = res_y ? 0 : 1;
int rvx0 = res_vx ? 0 : -1, rvx1 = res_vx ? 0 : 1;
int rvy0 = res_vy ? 0 : -1, rvy1 = res_vy ? 0 : 1;
// First pass: total in-bounds weight of the adaptive kernel
// Total in-bounds weight of the adaptive kernel
double sum_w = 0;
for (int dpy = ry0; dpy <= ry1; dpy++) {
int ny = py + dpy;
......@@ -216,8 +285,10 @@ public class CuasRecurrentBuffer {
}
}
if (sum_w <= 0) continue;
// Second pass: spread m with total weight 1
double scale = m / sum_w;
// Spread this survivor's share of the cluster energy: total = M*m_i/sum(m_j),
// so the whole cluster injects exactly its maximum M (continuous when a
// sub-pixel target shift changes the number of surviving cells)
double scale = (cl[1] * m / cl[2]) / sum_w;
for (int dpy = ry0; dpy <= ry1; dpy++) {
int ny = py + dpy;
if ((ny < 0) || (ny >= roi_height)) continue;
......
src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java
View file @ 2eb14ba0
......@@ -1156,6 +1156,7 @@ min_str_neib_fpn 0.35
public double curt_recur_leak = 1.0; // extra decay of the shifted prediction per update (gap decay = leak*(1-w) for fixed w) // By Claude on 06/09/2026
public double curt_recur_wthresh = 0.0; // if >0: adaptive per-pixel weight w_eff=w*tanh(strength/wthresh); 0 - fixed w // By Claude on 06/09/2026
public double curt_recur_spread = 0.5; // orthogonal 4D neighbor weight in adaptive spreading (diagonals get products) // By Claude on 06/09/2026
public int curt_recur_max_clust = 5; // max 4D-connected survivor cluster size treated as one target (energy-normalized spreading); larger - linear pass-through // By Claude on 06/09/2026
// debug/saving images
public boolean curt_save_c5full = false; // save fine velocities [direction][scene][subpixels]
public boolean curt_save_c5rect = true; // save fine velocities for selected rectangle only [scene][flattened image]
......@@ -3445,6 +3446,8 @@ min_str_neib_fpn 0.35
"If >0: per-pixel adaptive weight w_eff = w*tanh(strength/wthresh), dark frames give w_eff~0 (pure prediction). 0 - fixed w."); // By Claude on 06/09/2026
gd.addNumericField("4D spreading weight", this.curt_recur_spread, 6,8,"", // By Claude on 06/09/2026
"Weight of orthogonal 4D neighbors when spreading each surviving measurement point (diagonals get products of this)."); // By Claude on 06/09/2026
gd.addNumericField("Max target cluster size", this.curt_recur_max_clust, 0,3,"", // By Claude on 06/09/2026
"Survivor clusters up to this size are one target: spread with total energy = cluster max (compensates non-energy-preserving filtering). Larger clusters (weak linear signal) pass through unchanged."); // By Claude on 06/09/2026
gd.addMessage("=== Debug, saving images ===");
gd.addCheckbox ("Save full velocities hyperstack", this.curt_save_c5full,
......@@ -4957,6 +4960,7 @@ min_str_neib_fpn 0.35
this.curt_recur_leak = gd.getNextNumber(); // By Claude on 06/09/2026
this.curt_recur_wthresh = gd.getNextNumber(); // By Claude on 06/09/2026
this.curt_recur_spread = gd.getNextNumber(); // By Claude on 06/09/2026
this.curt_recur_max_clust=(int) gd.getNextNumber(); // By Claude on 06/09/2026
this.curt_save_c5full = gd.getNextBoolean();
this.curt_save_c5rect = gd.getNextBoolean();
......@@ -6294,6 +6298,7 @@ min_str_neib_fpn 0.35
properties.setProperty(prefix+"curt_recur_leak", this.curt_recur_leak+""); // double // By Claude on 06/09/2026
properties.setProperty(prefix+"curt_recur_wthresh", this.curt_recur_wthresh+""); // double // By Claude on 06/09/2026
properties.setProperty(prefix+"curt_recur_spread", this.curt_recur_spread+""); // double // By Claude on 06/09/2026
properties.setProperty(prefix+"curt_recur_max_clust", this.curt_recur_max_clust+"");// int // By Claude on 06/09/2026
properties.setProperty(prefix+"curt_save_c5full", this.curt_save_c5full+""); // boolean
properties.setProperty(prefix+"curt_save_c5rect", this.curt_save_c5rect+""); // boolean
......@@ -7618,6 +7623,7 @@ min_str_neib_fpn 0.35
if (properties.getProperty(prefix+"curt_recur_leak")!=null) this.curt_recur_leak=Double.parseDouble(properties.getProperty(prefix+"curt_recur_leak")); // By Claude on 06/09/2026
if (properties.getProperty(prefix+"curt_recur_wthresh")!=null) this.curt_recur_wthresh=Double.parseDouble(properties.getProperty(prefix+"curt_recur_wthresh")); // By Claude on 06/09/2026
if (properties.getProperty(prefix+"curt_recur_spread")!=null) this.curt_recur_spread=Double.parseDouble(properties.getProperty(prefix+"curt_recur_spread")); // By Claude on 06/09/2026
if (properties.getProperty(prefix+"curt_recur_max_clust")!=null) this.curt_recur_max_clust=Integer.parseInt(properties.getProperty(prefix+"curt_recur_max_clust")); // By Claude on 06/09/2026
if (properties.getProperty(prefix+"curt_save_c5full")!=null) this.curt_save_c5full=Boolean.parseBoolean(properties.getProperty(prefix+"curt_save_c5full"));
if (properties.getProperty(prefix+"curt_save_c5rect")!=null) this.curt_save_c5rect=Boolean.parseBoolean(properties.getProperty(prefix+"curt_save_c5rect"));
......@@ -8918,6 +8924,7 @@ min_str_neib_fpn 0.35
imp.curt_recur_leak = this.curt_recur_leak; // By Claude on 06/09/2026
imp.curt_recur_wthresh = this.curt_recur_wthresh; // By Claude on 06/09/2026
imp.curt_recur_spread = this.curt_recur_spread; // By Claude on 06/09/2026
imp.curt_recur_max_clust = this.curt_recur_max_clust; // By Claude on 06/09/2026
imp.curt_save_c5full = this.curt_save_c5full;
imp.curt_save_c5rect = this.curt_save_c5rect;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment