Skip to content

I
imagej-elphel
Commit 6a274e46 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

CLAUDE: robust average v2: NO NaN frame - all contributors + median/IQR rejection 

Andrey's correction: the eroded average left a wide NaN frame (and all
its descendants - EVAL diffs). The average virtual image must be defined
wherever ANY scene provides (a single bad first/last row/col at most):
combine ALL available contributors, let the count taper naturally near
the edges, and reject junk providers per pixel as OUTLIERS vs the MEDIAN
with an IQR-based robust sigma (immune to heavy contamination, unlike
mean/sigma). Fewer than 4 providers - plain mean; zero - NaN. Geometric
erosion (erodeDefined/erodedMasks/AVG_EDGE_ERODE) removed.

Next (Andrey): per-TILE outlier rejection in the 16-sensor TD
consolidation, anticipating the GPU migration (candidate indicator:
CC(0,0) flux, possibly + quadrant energy).
Co-Authored-By: 's avatarClaude Fable 5 <noreply@anthropic.com>
parent 8726cc75
Hide whitespace changes
Inline Side-by-side
Showing with 57 additions and 106 deletions
src/main/java/com/elphel/imagej/cuas/rt/CuasRender.java
View file @ 6a274e46
......@@ -606,9 +606,8 @@ public class CuasRender {
return rendered;
}
public static final int AVG_EDGE_ERODE = 24; // px: erode each scene's DEFINED region for the AVERAGE slices only (Andrey 07/07/2026) - partially covered virtual tiles at the warped-coverage edge emit tile-pitch junk (accumulated worst columns x=12,20,28 = 8-px pitch in ROWCOL-04); 3 tile pitches clears the partial ring. Individual scene slices keep everything (liberal). // By Claude on 07/07/2026
public static final double AVG_OUTLIER_SIGMA = 3.0; // reject contributions beyond this many sigmas from the first-pass mean... // By Claude on 07/07/2026
public static final int AVG_OUTLIER_MIN = 4; // ...when the pixel has at least this many (eroded) providers. // By Claude on 07/07/2026
public static final double AVG_OUTLIER_SIGMA = 3.0; // AVERAGE slices (Andrey's spec 07/07/2026 v2 - NO NaN frame, NO geometric erosion): combine ALL available contributors everywhere, contributor count tapers naturally near the edges; junk providers (partial-tile edge junk, x=12/20/28 tile-pitch class) are rejected per pixel as OUTLIERS vs the MEDIAN with this many robust sigmas (IQR/1.349 - junk-proof, unlike mean/sigma which contamination inflates)... // By Claude on 07/07/2026
public static final int AVG_OUTLIER_MIN = 4; // ...when the pixel has at least this many providers (fewer -> plain mean of what is there). NaN only where NO scene provides at all. // By Claude on 07/07/2026
/** Save a render-product hyperstack organized like -CUAS-INDIVIDUAL-CUAS-DBG so
* ImageJ can subtract the two directly (Andrey 07/06/2026):
......@@ -616,11 +615,10 @@ public class CuasRender {
* z (frames): "average" (all scenes), "average-f:l" (middle-half range,
* the renderSceneSequence default f=last/4, l=3*last/4),
* then the scenes.
* AVERAGE slices (Andrey 07/07/2026): per-scene contributions are edge-ERODED
* (AVG_EDGE_ERODE) and outlier-REJECTED (AVG_OUTLIER_SIGMA when >= AVG_OUTLIER_MIN
* providers) - the per-scene edge junk otherwise accumulates into stable margin
* artifacts; zero eroded providers -> honest NaN. Individual scene slices are saved
* as rendered (liberal - full coverage for inspection).
* AVERAGE slices (Andrey 07/07/2026 v2): ALL available contributors combined
* everywhere (count tapers naturally near the edges - no NaN frame), per-pixel
* MEDIAN/IQR outlier rejection kills the junk providers (see AVG_OUTLIER_*).
* Individual scene slices are saved as rendered (liberal - full coverage).
* By Claude on 07/06/2026 (extracted + reorganized from testRenderSequence). */
private static void saveRenderHyper(
final QuadCLT center_CLT,
......@@ -641,11 +639,9 @@ public class CuasRender {
// merged (rendered index num_comps-1) goes FIRST, then s00..s15
final int nr = (nc == 0) ? (num_comps - 1) : (nc - 1);
final String comp = (nc == 0) ? "merged" : String.format("s%02d", nc - 1);
// eroded per-scene validity for the AVERAGE slices (see class constants)
final boolean [][] eroded = erodedMasks(rendered, nr, img_width);
for (int navg = 0; navg < avg_ranges.length; navg++) {
stack.addSlice(comp+":"+avg_names[navg],
robustAverage(rendered, nr, avg_ranges[navg][0], avg_ranges[navg][1], eroded, npix));
robustAverage(rendered, nr, avg_ranges[navg][0], avg_ranges[navg][1], npix));
}
for (int ns = 0; ns < rendered.size(); ns++) {
final float [][] r = rendered.get(ns);
......@@ -661,7 +657,7 @@ public class CuasRender {
center_CLT.saveImagePlusInModelDirectory(null, imp); // title as filename
System.out.println("testRenderSequence(): saved "+suffix+" ("+num_comps+" components x "+
num_frames+" frames = 2 averages + "+rendered.size()+" scenes; merged first; "+
"averages: edge-eroded "+AVG_EDGE_ERODE+"px + "+AVG_OUTLIER_SIGMA+"-sigma outlier-rejected)");
"averages: all contributors, median/IQR "+AVG_OUTLIER_SIGMA+"-sigma outlier-rejected)");
}
/**
......@@ -694,8 +690,7 @@ public class CuasRender {
img_width, // int width,
img_height); // int height
final double [] kernel2d = cuasRTUtils.getKernel2d();
final boolean [][] eroded = erodedMasks(rendered, nr, img_width);
final float [] avg = robustAverage(rendered, nr, 0, rendered.size()-1, eroded, npix);
final float [] avg = robustAverage(rendered, nr, 0, rendered.size()-1, npix);
final int num_z = rendered.size() + 1; // average first
// per z-frame: [0]=log, [1]=log-diff, [2]=diff (computed once, added comp-major below)
final float [][][] comps = new float [3][num_z][];
......@@ -751,112 +746,68 @@ public class CuasRender {
return out;
}
/** Per-scene eroded validity masks for the robust averages (see class constants).
/** Robust temporal average of one component over [ns_from..ns_to] (Andrey's spec
* 07/07/2026 v2): ALL available contributors are combined everywhere - the count
* tapers naturally near the edges, NO geometric erosion and NO NaN frame. Junk
* providers are rejected per pixel as outliers vs the MEDIAN with an IQR-based
* robust sigma (immune to heavy contamination, unlike mean/sigma). Fewer than
* AVG_OUTLIER_MIN providers - plain mean; zero providers - NaN.
* By Claude on 07/07/2026. */
private static boolean [][] erodedMasks(
private static float [] robustAverage(
final ArrayList<float [][]> rendered,
final int nr,
final int img_width) {
final boolean [][] eroded = new boolean [rendered.size()][];
final int ns_from,
final int ns_to,
final int npix) {
final int nsc = ns_to - ns_from + 1;
final float [][] slices = new float [nsc][];
for (int ns = 0; ns < nsc; ns++) {
final float [][] r = rendered.get(ns_from + ns);
slices[ns] = ((r != null) && (r[nr] != null)) ? r[nr] : null;
}
final float [] avg = new float [npix];
final int STRIPE = 4096;
final int num_stripes = (npix + STRIPE - 1)/STRIPE;
final Thread[] threads = ImageDtt.newThreadArray();
final java.util.concurrent.atomic.AtomicInteger ai = new java.util.concurrent.atomic.AtomicInteger(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ns = ai.getAndIncrement(); ns < rendered.size(); ns = ai.getAndIncrement()) {
final float [][] r = rendered.get(ns);
if ((r != null) && (r[nr] != null)) {
eroded[ns] = erodeDefined(r[nr], img_width, AVG_EDGE_ERODE);
final float [] vals = new float [nsc];
for (int stripe = ai.getAndIncrement(); stripe < num_stripes; stripe = ai.getAndIncrement()) {
final int i1 = Math.min(npix, (stripe+1)*STRIPE);
for (int i = stripe*STRIPE; i < i1; i++) {
int n = 0;
for (int ns = 0; ns < nsc; ns++) if (slices[ns] != null) {
final float v = slices[ns][i];
if (!Float.isNaN(v)) vals[n++] = v;
}
if (n == 0) {
avg[i] = Float.NaN;
continue;
}
double s = 0.0;
if (n < AVG_OUTLIER_MIN) { // too few providers to judge outliers - plain mean
for (int j = 0; j < n; j++) s += vals[j];
avg[i] = (float) (s/n);
} else {
java.util.Arrays.sort(vals, 0, n);
final double med = 0.5*(vals[(n-1)/2] + vals[n/2]);
final double iqr = vals[(3*n)/4] - vals[n/4];
final double thr = AVG_OUTLIER_SIGMA * Math.max(iqr/1.349, 1.0E-12); // robust sigma; tiny floor keeps ==median survivors when IQR==0
int nk = 0;
for (int j = 0; j < n; j++) if (Math.abs(vals[j] - med) <= thr) {
s += vals[j];
nk++;
}
avg[i] = (nk > 0) ? ((float) (s/nk)) : ((float) med);
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return eroded;
}
/** Robust temporal average of one component over [ns_from..ns_to]: eroded per-scene
* contributions (partial-tile edge junk excluded by geometry), then outlier-rejected
* (AVG_OUTLIER_SIGMA vs the first-pass mean when >= AVG_OUTLIER_MIN providers);
* zero eroded providers -> NaN. By Claude on 07/07/2026. */
private static float [] robustAverage(
final ArrayList<float [][]> rendered,
final int nr,
final int ns_from,
final int ns_to,
final boolean [][] eroded,
final int npix) {
// pass 1: mean + sigma from the ERODED contributions
final double [] sum = new double [npix];
final double [] sum2 = new double [npix];
final int [] cnt = new int [npix];
for (int ns = ns_from; ns <= ns_to; ns++) {
final float [][] r = rendered.get(ns);
if ((r == null) || (r[nr] == null) || (eroded[ns] == null)) continue;
final float [] px = r[nr];
final boolean [] ok = eroded[ns];
for (int i = 0; i < npix; i++) if (ok[i]) {
final double d = px[i];
sum[i] += d;
sum2[i] += d*d;
cnt[i]++;
}
}
// pass 2: outlier-rejected mean (rejection only with enough providers)
final double [] sumk = new double [npix];
final int [] cntk = new int [npix];
final double [] mean = new double [npix];
final double [] rej = new double [npix]; // rejection threshold (|x-mean|), infinite when too few providers
for (int i = 0; i < npix; i++) {
if (cnt[i] > 0) {
mean[i] = sum[i]/cnt[i];
final double var = sum2[i]/cnt[i] - mean[i]*mean[i];
rej[i] = (cnt[i] >= AVG_OUTLIER_MIN) ?
(AVG_OUTLIER_SIGMA * Math.sqrt(Math.max(var, 0.0))) : Double.POSITIVE_INFINITY;
}
}
for (int ns = ns_from; ns <= ns_to; ns++) {
final float [][] r = rendered.get(ns);
if ((r == null) || (r[nr] == null) || (eroded[ns] == null)) continue;
final float [] px = r[nr];
final boolean [] ok = eroded[ns];
for (int i = 0; i < npix; i++) if (ok[i] && (Math.abs(px[i] - mean[i]) <= rej[i])) {
sumk[i] += px[i];
cntk[i]++;
}
}
final float [] avg = new float [npix];
for (int i = 0; i < npix; i++) {
avg[i] = (cntk[i] > 0) ? ((float) (sumk[i]/cntk[i])) : Float.NaN;
}
return avg;
}
/** Box erosion of the DEFINED (non-NaN) region: true only where the whole
* (2r+1)x(2r+1) neighborhood is defined AND inside the frame. Separable:
* per-row prefix counts of NaN, then per-column prefix counts of the row result.
* By Claude on 07/07/2026. */
private static boolean [] erodeDefined(final float [] px, final int width, final int r) {
final int npix = px.length;
final int height = npix/width;
final boolean [] okr = new boolean [npix]; // row-eroded
final int [] pre = new int [width+1];
for (int y = 0; y < height; y++) {
final int row0 = y*width;
for (int x = 0; x < width; x++) pre[x+1] = pre[x] + (Float.isNaN(px[row0+x]) ? 1 : 0);
for (int x = r; x < width-r; x++) {
okr[row0+x] = (pre[x+r+1] - pre[x-r]) == 0;
}
}
final boolean [] out = new boolean [npix];
final int [] prec = new int [height+1];
for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) prec[y+1] = prec[y] + (okr[y*width+x] ? 0 : 1);
for (int y = r; y < height-r; y++) {
out[y*width+x] = (prec[y+r+1] - prec[y-r]) == 0;
}
}
return out;
}
}
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