/**
 **
 ** CuasRender.java - dedicated CUAS RT renderer: conditioned raw-jp4 scenes
 ** rendered on the virtual-camera uniform grid at borrowed poses.
 **
 ** Copyright (C) 2026 Elphel, Inc.
 **
 ** -----------------------------------------------------------------------------**
 **
 **  CuasRender.java is free software: you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License as published by
 **  the Free Software Foundation, either version 3 of the License, or
 **  (at your option) any later version.
 **
 **  This program is distributed in the hope that it will be useful,
 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 **  GNU General Public License for more details.
 **
 **  You should have received a copy of the GNU General Public License
 **  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ** -----------------------------------------------------------------------------**
 **
 */
package com.elphel.imagej.cuas.rt;

import java.util.ArrayList;

import com.elphel.imagej.cameras.CLTParameters;
import com.elphel.imagej.cuas.CuasMotion;
import com.elphel.imagej.gpu.GpuQuad;
import com.elphel.imagej.gpu.TpTask;
import com.elphel.imagej.tileprocessor.ErsCorrection;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.OpticalFlow;
import com.elphel.imagej.tileprocessor.QuadCLT;
import com.elphel.imagej.tileprocessor.TileProcessor; // By Claude on 07/07/2026 (fillNaNs - FG-cluster Laplacian fill)
import com.elphel.imagej.tileprocessor.TwoQuadCLT;

import ij.ImagePlus;
import ij.ImageStack;

/**
 * The dedicated CUAS RT renderer (design ruling 2026-07-05, spec in
 * imagej-elphel-internal handoffs/2026-07-05_cuas_rt_dataflow_and_grids.md par.6):
 * replaces the legacy multi-purpose renderSceneSequence() for CUAS RT.
 *  - Grid: ALWAYS the virtual-camera uniform grid (no transformFromVirtual /
 *    transformCameraVew / conditionInitialDS in this path).
 *  - Disparity: the raw center DSI (the -INTER-INTRA-LMA vintage, same source as
 *    reliable_tiles) - ONE disparity source for all of RT.
 *  - Pixels: raw /jp4/*.tiff + CuasConditioning (conditionSceneToGpu) - never
 *    QuadCLT image products.
 *  - QuadCLT instances: BORROW-only (poses + ERS rates from scenes_poses, geometry,
 *    tile dims). Any new implied QuadCLT dependency must be raised with Andrey.
 *  - All methods STATIC.
 *
 * First delivered piece: the MINIMAL ingest+render sequence test (testRenderSequence)
 * that certifies conditionSceneToGpu + the virtual-grid render chain by generating
 * per-sensor + merged per-scene stacks directly comparable to the oracle
 * -CUAS-INDIVIDUAL-CUAS-DBG / -CUAS-MERGED-CUAS-DBG products. Today's lesson: this
 * step was jumped over (conditionSceneToGpu had a single caller inside the pose test
 * and was never image-certified standalone) - test and debug HERE first, only then
 * return to the pose test.
 *
 * By Claude on 07/05/2026, from Andrey's design.
 */
public class CuasRender {
	public static final double [] ZERO3 = {0.0, 0.0, 0.0};

	/**
	 * CUAS disparity conditioning for RENDERING (static, QuadCLT-free - the legacy
	 * conditionInitialDS takes a QuadCLT instance and is out per the RT rules).
	 * Requirements (Andrey 2026-07-05): disparity must be DEFINED FOR ALL TILES - the
	 * sky has nothing to match but is the background for drones, so it renders at
	 * infinity (disparity 0); NaN is acceptable only outside the image. The raw plain
	 * DSI slice has no NaN but carries garbage (crazy tiles < -1 in the sky, negatives
	 * near high-contrast posts) - all non-physical: clamp to >= 0.
	 * SMALL FOREGROUND clusters (Andrey 07/07/2026): where a small foreground object (the
	 * roadside "bomb" sign - 3 tiles at disp ~22, later lost) stands in front of a distant
	 * background, BOTH disparities are legitimate - but for CUAS the far background is what
	 * matters: in some sequences the drone flies exactly behind the sign and must stay
	 * detectable. Connected clusters of local NEAR-outliers (disparity above the local
	 * 7x7 median by fg_over) of at most max_fg_tiles tiles are re-filled with zero-curvature
	 * (Laplacian, TileProcessor.fillNaNs) interpolation of the surrounding background.
	 * @param disparity plain disparity (slice 0 of -INTER-INTRA-LMA), NOT the LMA slice
	 *                  (which is NaN wherever the LMA did not fit - the whole sky)
	 * @param tilesX       tile grid width
	 * @param max_fg_tiles largest connected near-outlier cluster to fill; <=0 - disable
	 * @param fg_over      disparity above the local median that makes a tile a near-outlier (px)
	 * @return conditioned copy: NaN -> 0.0 (infinity), negatives -> 0.0, small FG filled
	 */
	public static double [] conditionDisparityCuas(
			final double [] disparity,
			final int       tilesX,
			final int       max_fg_tiles,
			final double    fg_over) {
		final double [] cond = new double [disparity.length];
		for (int i = 0; i < disparity.length; i++) {
			final double d = disparity[i];
			cond[i] = (Double.isNaN(d) || (d < 0.0)) ? 0.0 : d;
		}
		if ((max_fg_tiles <= 0) || (tilesX <= 0)) {
			return cond;
		}
		final int tilesY = cond.length / tilesX;
		// 1. local near-outlier candidates: disparity > (7x7 neighborhood median) + fg_over
		//    (a <=9-tile cluster stays a minority of the 49-tile window)
		final int med_rad = 3;
		final boolean [] cand = new boolean [cond.length];
		final double [] wnd = new double [(2*med_rad+1)*(2*med_rad+1)];
		for (int ty = 0; ty < tilesY; ty++) {
			for (int tx = 0; tx < tilesX; tx++) {
				int nw = 0;
				for (int dy = -med_rad; dy <= med_rad; dy++) {
					final int yy = ty + dy;
					if ((yy < 0) || (yy >= tilesY)) continue;
					for (int dx = -med_rad; dx <= med_rad; dx++) {
						final int xx = tx + dx;
						if ((xx < 0) || (xx >= tilesX)) continue;
						wnd[nw++] = cond[yy*tilesX + xx];
					}
				}
				java.util.Arrays.sort(wnd, 0, nw);
				cand[ty*tilesX + tx] = cond[ty*tilesX + tx] > (wnd[nw/2] + fg_over);
			}
		}
		// 2. connected (4-conn) candidate clusters of at most max_fg_tiles -> NaN holes
		final boolean [] visited = new boolean [cond.length];
		final int [] stack = new int [cond.length];
		final int [] comp =  new int [cond.length];
		final double [] holed = cond.clone();
		int num_clusters = 0, num_tiles = 0;
		for (int seed = 0; seed < cond.length; seed++) if (cand[seed] && !visited[seed]) {
			int sp = 0, nc = 0;
			stack[sp++] = seed;
			visited[seed] = true;
			while (sp > 0) {
				final int t = stack[--sp];
				comp[nc++] = t;
				final int ty = t/tilesX, tx = t%tilesX;
				if ((tx > 0)        && cand[t-1]      && !visited[t-1])      { visited[t-1]=true;      stack[sp++] = t-1; }
				if ((tx < tilesX-1) && cand[t+1]      && !visited[t+1])      { visited[t+1]=true;      stack[sp++] = t+1; }
				if ((ty > 0)        && cand[t-tilesX] && !visited[t-tilesX]) { visited[t-tilesX]=true; stack[sp++] = t-tilesX; }
				if ((ty < tilesY-1) && cand[t+tilesX] && !visited[t+tilesX]) { visited[t+tilesX]=true; stack[sp++] = t+tilesX; }
			}
			if (nc <= max_fg_tiles) {
				for (int i = 0; i < nc; i++) holed[comp[i]] = Double.NaN;
				num_clusters++;
				num_tiles += nc;
			}
		}
		if (num_clusters == 0) {
			return cond;
		}
		// 3. zero-curvature (Laplacian) fill of the holes from the surrounding background
		final double [] filled = TileProcessor.fillNaNs(
				holed,     // final double []  data,
				null,      // final boolean [] prohibit,
				tilesX,    // int              width,
				2*tilesX,  // final int        grow (large - fill everything)
				0.7,       // double           diagonal_weight, // relative to ortho
				100,       // int              num_passes,
				0.03);     // final double     max_rchange
		System.out.println("conditionDisparityCuas(): filled "+num_clusters+" small foreground cluster(s), "+
				num_tiles+" tile(s) total, with zero-curvature background disparity (max_fg_tiles="+
				max_fg_tiles+", fg_over="+fg_over+")");
		return filled;
	}

	/**
	 * Render ONE scene (already ingested to the GPU by CuasConditioning.conditionSceneToGpu)
	 * on the virtual-camera uniform grid at the given (borrowed or fitted) pose.
	 * No correlation - tasks + convert (erased CLT: NaN outside task tiles) + imclt.
	 * @param clt_parameters processing parameters (LPF sigmas, margin)
	 * @param image_dtt      ImageDtt bound to the center GPU (created once per sequence)
	 * @param center_CLT     virtual center (grid dims, ErsCorrection reference side)
	 * @param scene          scene instance (BORROW: geometry + its ErsCorrection with
	 *                       rates already set; image data must already be on the GPU)
	 * @param center_disparity raw center DSI disparity (virtual grid)
	 * @param scene_xyz      scene position (ZERO3 for CUAS rotation-only)
	 * @param scene_atr      scene orientation (borrowed/fitted)
	 * @param margin         tile margin
	 * @param mb_vectors     null (no motion-blur compensation) or [2][ntiles] per-tile blur
	 *                       vectors (for CUAS uniform over the scene: dpx/dt, dpy/dt from
	 *                       the wobble radius + RPM). With MB the convert runs TWICE (the
	 *                       shifted, negative-scale second pass removes the bolometer
	 *                       exponential tail; scale pair sums to 1, ratio from
	 *                       (1-pix time)/mb_tau, offset grows if scales would exceed
	 *                       mb_max_gain - all inside setInterTasksMotionBlur).
	 * @param mb_tau         bolometer time constant, s (imp.mb_tau, 0.008)
	 * @param mb_max_gain    max scale before the offset is stretched (imp.mb_max_gain_inter)
	 * @param debugLevel     debug level
	 * @return [17][width*height]: renders of the 16 sensors (grid-transformed to the
	 *         virtual grid) + [16] = the consolidated (NaN-aware weighted average)
	 *         merged render; null if no valid tasks
	 */
	public static float [][] renderSceneVirtual(
			final CLTParameters clt_parameters,
			final ImageDtt      image_dtt,
			final QuadCLT       center_CLT,
			final QuadCLT       scene,
			final double []     center_disparity,
			final double []     scene_xyz,
			final double []     scene_atr,
			final int           margin,
			final double [][]   mb_vectors,  // By Claude on 07/05/2026
			final double        mb_tau,      // By Claude on 07/05/2026
			final double        mb_max_gain, // By Claude on 07/05/2026
			final int           debugLevel) {
		final GpuQuad gpuQuad = center_CLT.getGPUQuad();
		// 1. virtual uniform grid -> scene grid at the pose (ERS per the instances' rates)
		final double [][] scene_pXpYD = OpticalFlow.transformToScenePxPyD(
				null, center_disparity, scene_xyz, scene_atr, scene, center_CLT);
		final double disparity_corr = clt_parameters.imp.disparity_corr + center_CLT.getDispInfinityRef();
		final int img_width = scene.getGeometryCorrection().getSensorWH()[0];
		if (mb_vectors != null) {
			// Motion-blur compensation: TWO task sets (original positive scale + shifted
			// negative scale), double convert_direct - the oracle machinery
			// (setInterTasksMotionBlur / interCorrTDMotionBlur). The FIRST convert erases
			// the CLT to NaN, so tiles absent from the task set reach the NaN-aware channel
			// consolidation (and the renders) as NaN - never stale data (Andrey 07/05/2026).
			final TpTask [][] tp_tasks2 = GpuQuad.setInterTasksMotionBlur(
					scene.getNumSensors(),
					img_width,
					false,                      // GPU calculates port coordinates
					scene_pXpYD,
					null,                       // selection: ALL tiles with valid disparity
					mb_tau,
					mb_max_gain,
					mb_vectors,
					scene.getGeometryCorrection(),
					disparity_corr,
					margin,
					null,                       // valid_tiles
					ImageDtt.THREADS_MAX);
			if ((tp_tasks2 == null) || (tp_tasks2[0] == null) || (tp_tasks2[0].length == 0)) {
				System.out.println("renderSceneVirtual(): no MB tasks for scene "+scene.getImageName());
				return null;
			}
			image_dtt.interCorrTDMotionBlur(
					clt_parameters.img_dtt,
					1,                          // erase_clt_first: NaN outside the task set // By Claude on 07/05/2026
					tp_tasks2,
					null,
					clt_parameters.gpu_sigma_r,
					clt_parameters.gpu_sigma_b,
					clt_parameters.gpu_sigma_g,
					clt_parameters.gpu_sigma_m,
					scene.isMonochrome() ? 1.0 : clt_parameters.gpu_sigma_rb_corr,
					clt_parameters.getGpuCorrSigma(scene.isMonochrome()),
					clt_parameters.getGpuCorrLoGSigma(scene.isMonochrome()),
					clt_parameters.corr_red,
					clt_parameters.corr_blue,
					0,                          // sensor_mask_inter = 0: convert only
					ImageDtt.THREADS_MAX,
					debugLevel);
		} else {
			final TpTask [] tp_tasks = GpuQuad.setInterTasks(
					scene.getNumSensors(),
					img_width,
					false,                          // GPU calculates port coordinates
					scene_pXpYD,                    // per-tile pX,pY,disparity
					null,                           // selection: ALL tiles with valid disparity
					scene.getGeometryCorrection(),
					disparity_corr,
					margin,
					null,                           // valid_tiles
					ImageDtt.THREADS_MAX);
			if ((tp_tasks == null) || (tp_tasks.length == 0)) {
				System.out.println("renderSceneVirtual(): no tasks for scene "+scene.getImageName());
				return null;
			}
			// 2. LPFs + tasks + offsets + convert (mask=0: no correlation)
			image_dtt.interCorrTD(
					clt_parameters.img_dtt,
					tp_tasks,
					null,
					clt_parameters.gpu_sigma_r,
					clt_parameters.gpu_sigma_b,
					clt_parameters.gpu_sigma_g,
					clt_parameters.gpu_sigma_m,
					scene.isMonochrome() ? 1.0 : clt_parameters.gpu_sigma_rb_corr,
					clt_parameters.getGpuCorrSigma(scene.isMonochrome()),
					clt_parameters.getGpuCorrLoGSigma(scene.isMonochrome()),
					clt_parameters.corr_red,
					clt_parameters.corr_blue,
					0,                              // sensor_mask_inter = 0: convert only
					ImageDtt.THREADS_MAX,
					debugLevel);
			// re-convert with erase_clt=1: NaN outside task tiles (interCorrTD converts with
			// erase=-1, leaving the previous scene's tiles as ghosts in a render)
			gpuQuad.execConvertDirect(false, null, 1);
		}
		// 3. per-sensor renders (virtual grid)
		final float [][] per_sensor = CuasMotion.perSensorImagesFromTD(gpuQuad, false);
		// 4. merged = the NaN-aware weighted average consolidated in TD, rendered via slot 0
		final float [][] fclt16 = gpuQuad.getCltData(false);
		final float [] avg_td = CuasTD.consolidateSensorsTD(fclt16, null);
		gpuQuad.setCltData(0, avg_td, false);
		final float [] merged = CuasMotion.perSensorImagesFromTD(gpuQuad, false)[0];
		final float [][] rslt = new float [per_sensor.length + 1][];
		System.arraycopy(per_sensor, 0, rslt, 0, per_sensor.length);
		rslt[per_sensor.length] = merged;
		return rslt;
	}

	/**
	 * MINIMAL standalone test of the RT ingest+render chain (runs INSTEAD of detection
	 * and INSTEAD of the pose test): for every scene, ingest raw /jp4/ via
	 * CuasConditioning.conditionSceneToGpu, borrow the STORED pose + STORED ERS rates
	 * (exactly the sources the oracle DBG renders use), render on the virtual uniform
	 * grid, and save one hyperstack:
	 *   <center>-CUAS-RT-RENDER.tiff  [t: s00..s15, merged][z: scenes]
	 * directly comparable (slice-by-slice, same borrowed poses) to the oracle
	 * -CUAS-INDIVIDUAL-CUAS-DBG / -CUAS-MERGED-CUAS-DBG debugging products.
	 * By Claude on 07/05/2026, from Andrey's plan: test and debug THIS first,
	 * only then return to the pose test.
	 */
	public static void testRenderSequence(
			final CLTParameters clt_parameters,
			final QuadCLT       center_CLT,
			final QuadCLT []    quadCLTs,
			final int           debugLevel) {
		final int margin = clt_parameters.imp.margin;
		final boolean is_aux = center_CLT.isAux();
		final double [][] center_dsi = center_CLT.dsi;
		// RENDERING disparity: the PLAIN slice conditioned to be defined for ALL tiles
		// (sky = infinity/0). Deliberately NOT gated by use_lma_dsi: the LMA slice is NaN
		// wherever the LMA did not fit (the whole sky) and would drop those tiles from the
		// render. By Claude on 07/05/2026, per Andrey.
		final double [] center_disparity = conditionDisparityCuas(
				center_dsi[is_aux ? TwoQuadCLT.DSI_DISPARITY_AUX : TwoQuadCLT.DSI_DISPARITY_MAIN],
				center_CLT.getTilesX(),               // By Claude on 07/07/2026 (small-FG Laplacian fill)
				clt_parameters.curt.disp_fg_tiles,
				clt_parameters.curt.disp_fg_over);
		final ErsCorrection ers_center = center_CLT.getErsCorrection();
		final ImageDtt image_dtt = new ImageDtt(
				center_CLT.getNumSensors(),
				clt_parameters.transform_size,
				clt_parameters.img_dtt,
				center_CLT.isAux(),
				center_CLT.isMonochrome(),
				center_CLT.isLwir(),
				clt_parameters.getScaleStrength(center_CLT.isAux()),
				center_CLT.getGPU());
		final int img_width =  center_CLT.getGPUQuad().getImageWidth();
		final int img_height = center_CLT.getGPUQuad().getImageHeight();
		final int num_sens =   center_CLT.getNumSensors();
		final int num_comps =  num_sens + 1; // 16 sensors + merged
		// Motion-blur compensation, gated by the same imp.mb_en as the oracle DBG renders
		// (mb_en OFF -> compare against the -NOMB oracle products). Vectors from the same
		// getMotionBlur the oracle uses, at the borrowed pose + stored rates. The identity
		// center grid is needed for it. By Claude on 07/05/2026.
		final boolean mb_en =       clt_parameters.imp.mb_en;
		final double  mb_tau =      clt_parameters.imp.mb_tau;
		final double  mb_max_gain = clt_parameters.imp.mb_max_gain_inter;
		final double [][] pXpYD_center = mb_en ? OpticalFlow.transformToScenePxPyD(
				null, center_disparity, ZERO3, ZERO3, center_CLT, center_CLT) : null;
		System.out.println("testRenderSequence(): motion-blur compensation "+
				(mb_en ? ("ON (tau="+mb_tau+" s, max gain "+mb_max_gain+") - compare vs the MB oracle DBG") :
					"OFF - compare vs the -NOMB oracle DBG"));
		// PER-SEQUENCE DC level for the scene conditioning - NOT per-scene: a per-scene sky
		// median oscillates with rotation (terrain sweeps the sky half; 234-count pk-pk at
		// exactly 1 rev period measured in THIS product, POSE-10 forensics) and downstream
		// temporal averages turn it into content ghosts. ONE constant per sequence
		// (CuasRT.appropriateCenter derived it from the borrowed center = the average
		// scene). By Claude on 07/05/2026, Andrey's rule.
		final CuasConditioning.Config cond_cfg = CuasConditioning.configFromCurt(clt_parameters.curt); // rowcol_* now dialog-driven (RT step 2a) // By Claude on 07/07/2026
		cond_cfg.dc_level = CuasRT.sequenceDcLevel(clt_parameters, center_CLT);
		if (Double.isNaN(cond_cfg.dc_level)) {
			System.out.println("testRenderSequence(): WARNING - no per-sequence DC level (RT-STATE missing, curt.dc_center==0):"+
					" falling back to the PER-SCENE sky median (rotation-synced DC flicker!)");
		}
		final ArrayList<String>     ts_names = new ArrayList<String>();
		final ArrayList<float [][]> rendered = renderPass( // By Claude on 07/06/2026 (loop extracted for the LoG A/B second pass)
				clt_parameters, center_CLT, quadCLTs, image_dtt, ers_center,
				center_disparity, pXpYD_center, cond_cfg, margin,
				mb_en, mb_tau, mb_max_gain, ts_names, debugLevel);
		saveRenderHyper(center_CLT, "-CUAS-RT-RENDER", rendered, ts_names,
				img_width, img_height, num_comps);
		// ---- Evaluation product (Andrey 07/07/2026, replaces the manual ImageJ UM
		// workflow): -CUAS-RT-RENDER-EVAL, merged component only, z = average + scenes,
		// t (auto-range order): 1) LoG (no subtraction), 2) LoG of (scene - robust
		// average), 3) difference without LoG; t2/t3 are 0 for the average frame.
		saveRenderEvalLog(clt_parameters, center_CLT, rendered, ts_names,
				img_width, img_height, num_comps); // By Claude on 07/07/2026
		// ---- LoG A/B (curt.log_test / curt.log_ident): the aberration+LoG
		// kernel-fold validation (roadmap RT-seed step 1b). Comparable hyperstacks:
		//   log_test:  -LOG-ORACLE (product above x the EXISTING pre-DNN pixel LoG,
		//              convolve2DLReLU LINEAR alpha=1) + -LOG-FOLDED (same render,
		//              LoG folded into the aberration kernels - no pixel conv).
		//   log_ident: -ID (IDENTITY kernels - aberrations skipped) + -ID-LOG-ORACLE
		//              (pixel LoG of -ID) + -LOG-ONLY (pure-LoG kernels). Isolates
		//              the fold mechanism from the aberration-kernel content
		//              (Andrey 07/06/2026: if -LOG-ONLY vs -ID-LOG-ORACLE residual
		//              persists, the fold/L0 is wrong; if it collapses to the wrap
		//              floor, the aberration-kernel interaction was responsible).
		// Original kernels ALWAYS restored. Offline verdict:
		// attic/CLAUDE/compare_cuda_versions.py --only LOG-,ID. By Claude on 07/06/2026.
		if (clt_parameters.curt.log_test || clt_parameters.curt.log_ident) {
			final GpuQuad gpuQuad = center_CLT.getGPUQuad();
			final double [][][][][][] orig_kernels = gpuQuad.getQuadCLT().getCLTKernels();
			if (orig_kernels == null) {
				System.out.println("testRenderSequence(): LoG A/B ERROR - getCLTKernels()==null (kernels never loaded?), skipping LoG tests");
				return;
			}
			final double [] log_td = CuasLogFold.getLogTd(
					clt_parameters.curt.psf_radius,
					clt_parameters.curt.n_sigma);
			// EXACT pre-DNN LoG: same class, same construction as CuasDetectRT.detectTargets
			final CuasRTUtils cuasRTUtils = new CuasRTUtils (
					clt_parameters.curt.psf_radius,     // double   psf_radius,
					clt_parameters.curt.n_sigma,        // double   n_sigma,
					clt_parameters.curt.temp_weights.clone(), // double[] temporal_weights,
					clt_parameters.curt.pix_decimate,   // int      pixel_decimate,
					clt_parameters.curt.vel_decimate,   // int      velocity_decimate,
					clt_parameters.curt.vel_radius,     // int      velocity_radius,
					clt_parameters.curt.vel_suppr_rad,  // int      vel_suppr_rad,
					img_width,                          // int      width,
					img_height);                        // int      height
			try {
				if (clt_parameters.curt.log_test) {
					System.out.println("testRenderSequence(): LoG A/B - applying pixel-domain LoG (LINEAR) for -LOG-ORACLE");
					applyPixelLog(rendered, cuasRTUtils); // in place: originals already saved
					saveRenderHyper(center_CLT, "-CUAS-RT-RENDER-LOG-ORACLE", rendered, ts_names,
							img_width, img_height, num_comps);
					rendered.clear(); // release ~11GB before the next render pass
					System.out.println("testRenderSequence(): LoG A/B - uploading FOLDED aberration+LoG kernels (CuasLogFold, tap-precompensated)");
					double [][][][][][] folded_kernels=CuasLogFold.foldKernels(orig_kernels, log_td);
					if (clt_parameters.curt.show_kernels) {
						GpuQuad.showCLtKernel(
								folded_kernels, // double [][][][][][] clt_kernels,
								"folded-kernel", // String prefix,
								0, // int    nsens,
								0, // int    ncol,
								clt_parameters.tileX/2, // int    ktilex, // half tileX
								clt_parameters.tileY/2); // int    ktiley) { // half tileY
					}  
//					gpuQuad.setConvolutionKernels(CuasLogFold.foldKernels(orig_kernels, log_td), true);
					gpuQuad.setConvolutionKernels(folded_kernels, true);
					final ArrayList<String> ts_names2 = new ArrayList<String>();
					final ArrayList<float [][]> rendered2 = renderPass(
							clt_parameters, center_CLT, quadCLTs, image_dtt, ers_center,
							center_disparity, pXpYD_center, cond_cfg, margin,
							mb_en, mb_tau, mb_max_gain, ts_names2, debugLevel);
					saveRenderHyper(center_CLT, "-CUAS-RT-RENDER-LOG-FOLDED", rendered2, ts_names2,
							img_width, img_height, num_comps);
				}
				System.out.println("testRenderSequence(): LoG A/B - uploading IDENTITY kernels (aberrations skipped)");
//				gpuQuad.setConvolutionKernels(CuasLogFold.uniformKernels(orig_kernels, null), true);
				double [][][][][][] uniform_kernels=CuasLogFold.uniformKernels(orig_kernels, null);
				if (clt_parameters.curt.show_kernels) {
					GpuQuad.showCLtKernel(
							uniform_kernels, // double [][][][][][] clt_kernels,
							"uniform-kernel", // String prefix,
							0, // int    nsens,
							0, // int    ncol,
							clt_parameters.tileX/2, // int    ktilex, // half tileX
							clt_parameters.tileY/2); // int    ktiley) { // half tileY
				}  
				gpuQuad.setConvolutionKernels(uniform_kernels, true);
				
				final ArrayList<String> ts_id = new ArrayList<String>();
				final ArrayList<float [][]> rend_id = renderPass(
						clt_parameters, center_CLT, quadCLTs, image_dtt, ers_center,
						center_disparity, pXpYD_center, cond_cfg, margin,
						mb_en, mb_tau, mb_max_gain, ts_id, debugLevel);
				saveRenderHyper(center_CLT, "-CUAS-RT-RENDER-ID", rend_id, ts_id,
						img_width, img_height, num_comps);
				System.out.println("testRenderSequence(): LoG A/B - applying pixel-domain LoG (LINEAR) for -ID-LOG-ORACLE");
				applyPixelLog(rend_id, cuasRTUtils);
				saveRenderHyper(center_CLT, "-CUAS-RT-RENDER-ID-LOG-ORACLE", rend_id, ts_id,
						img_width, img_height, num_comps);
				rend_id.clear();
				System.out.println("testRenderSequence(): LoG A/B - uploading PURE-LoG kernels (LoG-only, aberrations skipped)");
				uniform_kernels=CuasLogFold.uniformKernels(orig_kernels, log_td);
				if (clt_parameters.curt.show_kernels) {
					GpuQuad.showCLtKernel(
							uniform_kernels, // double [][][][][][] clt_kernels,
							"uniform-LoG-kernel", // String prefix,
							0, // int    nsens,
							0, // int    ncol,
							clt_parameters.tileX/2, // int    ktilex, // half tileX
							clt_parameters.tileY/2); // int    ktiley) { // half tileY
				}  
//				gpuQuad.setConvolutionKernels(CuasLogFold.uniformKernels(orig_kernels, log_td), true);
				gpuQuad.setConvolutionKernels(uniform_kernels, true);
				final ArrayList<String> ts_lo = new ArrayList<String>();
				final ArrayList<float [][]> rend_lo = renderPass(
						clt_parameters, center_CLT, quadCLTs, image_dtt, ers_center,
						center_disparity, pXpYD_center, cond_cfg, margin,
						mb_en, mb_tau, mb_max_gain, ts_lo, debugLevel);
				saveRenderHyper(center_CLT, "-CUAS-RT-RENDER-LOG-ONLY", rend_lo, ts_lo,
						img_width, img_height, num_comps);
			} finally { // ALWAYS restore the original kernels for anything running after
				gpuQuad.setConvolutionKernels(orig_kernels, true);
				System.out.println("testRenderSequence(): LoG A/B - original kernels restored");
			}
			System.out.println("testRenderSequence(): LoG A/B done - compare the LOG products offline"+
					" (attic/CLAUDE/compare_cuda_versions.py)");
		}
	}

	/** Convolve every slice with the EXISTING pre-DNN pixel-domain LoG (LINEAR
	 *  alpha=1, NaN-aware), replacing slices in place. By Claude on 07/06/2026. */
	private static void applyPixelLog(
			final ArrayList<float [][]> rendered,
			final CuasRTUtils           cuasRTUtils) {
		final double [] kernel2d = cuasRTUtils.getKernel2d();
		for (int ns = 0; ns < rendered.size(); ns++) {
			final float [][] r = rendered.get(ns);
			if (r == null) continue;
			for (int nc = 0; nc < r.length; nc++) {
				if (r[nc] == null) continue;
				final double [] dslice = new double [r[nc].length];
				for (int i = 0; i < dslice.length; i++) dslice[i] = r[nc][i];
				final double [] dlog = cuasRTUtils.convolve2DLReLU(
						dslice,    // final double [] data,
						kernel2d,  // final double [] kernel, // square
						null,      // double []       result_in,
						1.0);      // final double    alpha) - LINEAR (the standing ruling)
				final float [] fslice = new float [dlog.length];
				for (int i = 0; i < fslice.length; i++) fslice[i] = (float) dlog[i];
				r[nc] = fslice;
			}
		}
	}

	/**
	 * One full render pass over the sequence (the loop extracted verbatim from
	 * testRenderSequence so the LoG A/B can run it twice - with the original and
	 * with the folded kernels). By Claude on 07/06/2026.
	 */
	private static ArrayList<float [][]> renderPass(
			final CLTParameters clt_parameters,
			final QuadCLT       center_CLT,
			final QuadCLT []    quadCLTs,
			final ImageDtt      image_dtt,
			final ErsCorrection ers_center,
			final double []     center_disparity,
			final double [][]   pXpYD_center,
			final CuasConditioning.Config cond_cfg,
			final int           margin,
			final boolean       mb_en,
			final double        mb_tau,
			final double        mb_max_gain,
			final ArrayList<String> ts_names, // OUT: aligned scene names
			final int           debugLevel) {
		final ArrayList<float [][]> rendered = new ArrayList<float [][]>();
		int num_ok = 0, num_fail = 0;
		for (int nscene = 0; nscene < quadCLTs.length; nscene++) {
			if (quadCLTs[nscene] == null) continue;
			final String ts_name = quadCLTs[nscene].getImageName();
			// BORROWED pose + BORROWED stored ERS rates - the same sources the oracle
			// DBG renders use (renderSceneSequence: getSceneXYZ/ATR + getSceneErsXYZ_dt/ATR_dt)
			final double [][] pose = ers_center.getSceneXYZATR(ts_name);
			final double [] ers_xyz_dt = ers_center.getSceneErsXYZ_dt(ts_name);
			final double [] ers_atr_dt = ers_center.getSceneErsATR_dt(ts_name);
			float [][] rslt = null;
			if (pose == null) {
				System.out.println("testRenderSequence(): scene "+nscene+" ("+ts_name+") has no stored pose - skipping");
			} else {
				if ((ers_xyz_dt == null) || (ers_atr_dt == null)) {
					System.out.println("testRenderSequence(): scene "+nscene+" ("+ts_name+") has no stored ERS rates - zeros");
				}
				// THE STEP UNDER TEST: raw /jp4/ -> conditioning -> GPU, motion (omegas +
				// velocities provision) entering through the ONE wrapper. Here the rates are
				// borrowed (stored); in production they come analytically from pose + rpm.
				final boolean raw_ok = CuasConditioning.conditionSceneToGpuCuas(
						quadCLTs[nscene],
						cond_cfg,             // Config: defaults + the PER-SEQUENCE dc_level // By Claude on 07/05/2026
						(ers_atr_dt != null) ? ers_atr_dt : ZERO3.clone(), // omegas
						(ers_xyz_dt != null) ? ers_xyz_dt : ZERO3.clone(), // velocities (provision)
						ImageDtt.THREADS_MAX,
						debugLevel - 2);
				if (raw_ok) {
					double [][] mb_vectors = null;
					if (mb_en) { // same oracle machinery + same (borrowed) rates as the DBG renders
						mb_vectors = OpticalFlow.getMotionBlur(
								center_CLT,          // QuadCLT        ref_scene,
								quadCLTs[nscene],    // QuadCLT        scene,
								pXpYD_center,        // double [][]    ref_pXpYD,
								pose[0],             // double []      camera_xyz,
								pose[1],             // double []      camera_atr,
								(ers_xyz_dt != null) ? ers_xyz_dt : ZERO3.clone(), // camera_xyz_dt,
								(ers_atr_dt != null) ? ers_atr_dt : ZERO3.clone(), // camera_atr_dt,
								0,                   // int            shrink_gaps,
								debugLevel - 2);     // int            debug_level
					}
					rslt = renderSceneVirtual(
							clt_parameters, image_dtt, center_CLT, quadCLTs[nscene],
							center_disparity, pose[0], pose[1], margin,
							mb_vectors, mb_tau, mb_max_gain, debugLevel - 2);
				} else {
					System.out.println("testRenderSequence(): scene "+nscene+" ("+ts_name+") raw-jp4 ingest FAILED");
				}
			}
			if (rslt != null) num_ok++; else num_fail++;
			rendered.add(rslt); // null OK - NaN slices keep z aligned
			ts_names.add(ts_name);
			if ((debugLevel > -4) && (((num_ok + num_fail) % 50) == 0)) {
				System.out.println("testRenderSequence(): rendered "+(num_ok+num_fail)+" scenes ("+num_fail+" failed)");
			}
		}
		System.out.println("testRenderSequence(): render pass done, scenes OK="+num_ok+", failed="+num_fail);
		return rendered;
	}

	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):
	 *    t (channels): merged FIRST, then s00..s15;
	 *    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 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,
			final String                suffix,
			final ArrayList<float [][]> rendered,
			final ArrayList<String>     ts_names,
			final int                   img_width,
			final int                   img_height,
			final int                   num_comps) {
		final int npix = img_width * img_height;
		final float [] nan_slice = new float [npix];
		java.util.Arrays.fill(nan_slice, Float.NaN);
		final int last_scene = rendered.size() - 1;
		final int [][] avg_ranges = {{0, last_scene}, {last_scene / 4, (3 * last_scene) / 4}};
		final String [] avg_names = {"average", "average-"+avg_ranges[1][0]+":"+avg_ranges[1][1]};
		final ImageStack stack = new ImageStack(img_width, img_height);
		for (int nc = 0; nc < num_comps; nc++) {
			// 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);
			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], npix));
			}
			for (int ns = 0; ns < rendered.size(); ns++) {
				final float [][] r = rendered.get(ns);
				stack.addSlice(comp+":"+ts_names.get(ns), ((r != null) && (r[nr] != null)) ? r[nr] : nan_slice);
			}
		}
		final int num_frames = rendered.size() + avg_ranges.length;
		final ImagePlus imp = new ImagePlus(center_CLT.getImageName()+suffix, stack);
		// (average slices above: eroded + outlier-rejected; see class constants)
		imp.setDimensions(1, num_frames, num_comps);
		imp.setOpenAsHyperStack(true);
		imp.getProcessor().resetMinAndMax();
		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: all contributors, median/IQR "+AVG_OUTLIER_SIGMA+"-sigma outlier-rejected)");
	}

	/**
	 * Evaluation hyperstack (Andrey's spec 07/07/2026): merged component only,
	 * z (frames) = "average" (the robust all-scenes average) + each scene,
	 * t (channels, ordered so ImageJ auto-ranges on the LoG view):
	 *   1) "log"      - pixel LoG (pre-DNN kernel, LINEAR) of the image itself;
	 *   2) "log-diff" - LoG of (scene - average); 0 for the average frame;
	 *   3) "diff"     - (scene - average) without LoG;   0 for the average frame.
	 * Replaces the manual run-UM-in-ImageJ evaluation workflow. By Claude on 07/07/2026.
	 */
	private static void saveRenderEvalLog(
			final CLTParameters         clt_parameters,
			final QuadCLT               center_CLT,
			final ArrayList<float [][]> rendered,
			final ArrayList<String>     ts_names,
			final int                   img_width,
			final int                   img_height,
			final int                   num_comps) {
		final int npix = img_width * img_height;
		final int nr = num_comps - 1; // merged
		final CuasRTUtils cuasRTUtils = new CuasRTUtils ( // same construction as CuasDetectRT.detectTargets
				clt_parameters.curt.psf_radius,     // double   psf_radius,
				clt_parameters.curt.n_sigma,        // double   n_sigma,
				clt_parameters.curt.temp_weights.clone(), // double[] temporal_weights,
				clt_parameters.curt.pix_decimate,   // int      pixel_decimate,
				clt_parameters.curt.vel_decimate,   // int      velocity_decimate,
				clt_parameters.curt.vel_radius,     // int      velocity_radius,
				clt_parameters.curt.vel_suppr_rad,  // int      vel_suppr_rad,
				img_width,                          // int      width,
				img_height);                        // int      height
		final double [] kernel2d = cuasRTUtils.getKernel2d();
		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][];
		comps[0][0] = logSlice(avg, kernel2d, cuasRTUtils);
		comps[1][0] = new float [npix]; // zeros for the average frame
		comps[2][0] = new float [npix];
		for (int ns = 0; ns < rendered.size(); ns++) {
			final float [][] r = rendered.get(ns);
			final float [] px = ((r != null) && (r[nr] != null)) ? r[nr] : null;
			final float [] diff = new float [npix];
			if (px != null) {
				for (int i = 0; i < npix; i++) {
					diff[i] = px[i] - avg[i]; // NaN where either is NaN
				}
				comps[0][ns+1] = logSlice(px,   kernel2d, cuasRTUtils);
				comps[1][ns+1] = logSlice(diff, kernel2d, cuasRTUtils);
				comps[2][ns+1] = diff;
			} else {
				java.util.Arrays.fill(diff, Float.NaN);
				comps[0][ns+1] = diff;
				comps[1][ns+1] = diff;
				comps[2][ns+1] = diff;
			}
		}
		final String [] comp_names = {"log", "log-diff", "diff"};
		final ImageStack stack = new ImageStack(img_width, img_height);
		for (int nc = 0; nc < comps.length; nc++) {
			stack.addSlice(comp_names[nc]+":average", comps[nc][0]);
			for (int ns = 0; ns < rendered.size(); ns++) {
				stack.addSlice(comp_names[nc]+":"+ts_names.get(ns), comps[nc][ns+1]);
			}
		}
		final ImagePlus imp = new ImagePlus(center_CLT.getImageName()+"-CUAS-RT-RENDER-EVAL", stack);
		imp.setDimensions(1, num_z, comps.length);
		imp.setOpenAsHyperStack(true);
		imp.getProcessor().resetMinAndMax();
		center_CLT.saveImagePlusInModelDirectory(null, imp); // title as filename
		System.out.println("testRenderSequence(): saved -CUAS-RT-RENDER-EVAL (merged only; t = log, log-diff, diff; "+
				"z = robust average + "+rendered.size()+" scenes)");
	}

	/** Pixel LoG (pre-DNN kernel, LINEAR alpha=1) of one float slice. By Claude on 07/07/2026. */
	private static float [] logSlice(final float [] px, final double [] kernel2d, final CuasRTUtils cuasRTUtils) {
		final double [] dslice = new double [px.length];
		for (int i = 0; i < dslice.length; i++) dslice[i] = px[i];
		final double [] dlog = cuasRTUtils.convolve2DLReLU(
				dslice,    // final double [] data,
				kernel2d,  // final double [] kernel, // square
				null,      // double []       result_in,
				1.0);      // final double    alpha) - LINEAR (the standing ruling)
		final float [] out = new float [px.length];
		for (int i = 0; i < out.length; i++) out[i] = (float) dlog[i];
		return out;
	}

	/** 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 float [] robustAverage(
			final ArrayList<float [][]> rendered,
			final int                   nr,
			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() {
					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 avg;
	}
}
