package com.elphel.imagej.cuas.rt;

import java.awt.Rectangle;
import java.util.Arrays;
import java.util.concurrent.atomic.AtomicInteger;

import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.cuas.rt.TemporalKernelGenerator.ConsolidationKernel;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.TileNeibs;

import ij.ImagePlus;

public class CuasRTUtils {
	public static double [] IMP_MINMAX= {-5,5}; //default range for
	final double         psf_radius;
	final double         n_sigma; // for LoG kernel size
	final double []      kernel2d;
	final private int    kernel2d_rad;
	// contains pixel indices, approximately in reverse bit order for multithreading efficiency
	// for the center area only, does not need to care about borders 
	final private int [] indx_center_2d;
	// remaining near-border pixel indices
	final private int [] indx_margins_2d;
	// all pixels indices
	final private int [] indx_all_2d;
	// same for 3x3 velocities
	final private int    kernel3d3_rad;
	final private int [] indx_all_3d3;
	final private int [] indx_center_3d3;
	final private int [] indx_margins_3d3;

	final private int    pixel_decimate;
	final private int    velocity_decimate;
	final private int    velocity_radius;
	
	final private double [] temporal_weights;
	final private double [][][][][] kernel5d;
	final private int    kernel5d_rad;
	final private int [] indx_all_5d;
	final private int [] indx_center_5d;
	final private int [] indx_margins_5d;
	
	final private int    vel_suppr_rad;
	final private int [] indx_all_vsupp;
	final private int [] indx_center_vsupp;
	final private int [] indx_margins_vsupp;
	
	final ConsolidationKernel consolidationKernel;

	// Direct pixel-domain fine-velocity kernel (bypasses the 3x3 coarse stage) // By Claude on 06/11/2026
	final ConsolidationKernel consolidationKernelDirect;
	final private int    kernel5dp_rad;
	final private int [] indx_all_5dp;
	final private int [] indx_center_5dp;
	final private int [] indx_margins_5dp;

	// balance kernel - center = +1, sum of other weights = -1.
	final private int    width;
	final private int    height;
//	final private int    tl_offset;
	public CuasRTUtils (
			double   psf_radius,
			double   n_sigma,
			int      kernel3d3_rad,
			double[] temporal_weights,
            int      pixel_decimate,
            int      velocity_decimate,
            int      velocity_radius,
            int      vel_suppr_rad,
			int      width,
			int      height) {
		this(psf_radius, n_sigma, kernel3d3_rad, temporal_weights, pixel_decimate, // By Claude on 06/12/2026
				velocity_decimate, velocity_radius, vel_suppr_rad, width, height, true); // default: matched direct kernel
	}

	public CuasRTUtils ( // By Claude on 06/12/2026
			double   psf_radius,
			double   n_sigma,
			int      kernel3d3_rad,
			double[] temporal_weights,
            int      pixel_decimate,
            int      velocity_decimate,
            int      velocity_radius,
            int      vel_suppr_rad,
			int      width,
			int      height,
			boolean  c5_matched) { // direct/C5P kernel: true = matched-filter half-cosine bump, false = trapezoid optical-analogy
		this.psf_radius = psf_radius;
		this.width =  width;
		this.height = height;
		this.n_sigma = n_sigma; 
		this.vel_suppr_rad=vel_suppr_rad;
		// Setup 2D pre-filter kernel, based on LoG
		kernel2d_rad = getRadiusLoGKernel(psf_radius,n_sigma);
		kernel2d = getLoGKenel(psf_radius, kernel2d_rad);
		int [][] indices = getIndices(
				kernel2d_rad,
				width,
				height);
		indx_all_2d =      indices[0];
		indx_center_2d =   indices[1];
		indx_margins_2d =  indices[2];
		
		// Setup 3D kernel to generate coarse velocities and input to the next pyramid level (center output)
		this.kernel3d3_rad = kernel3d3_rad;
		indices = getIndices(
				kernel3d3_rad,
				width,
				height);
		indx_all_3d3 =     indices[0];
		indx_center_3d3 =  indices[1];
		indx_margins_3d3 = indices[2];
		
		// Setup 5D kernel for generating fine velocities at increased spatial resolution
		this.temporal_weights = temporal_weights.clone();
		normalize1d (this.temporal_weights);
		this.pixel_decimate =    pixel_decimate;
		this.velocity_decimate = velocity_decimate;
		this.velocity_radius =   velocity_radius;

		consolidationKernel = TemporalKernelGenerator.generateKernel(
				this.temporal_weights,  // double[] temporal_weights,
				pixel_decimate,    // int pixel_decimate,
				velocity_decimate, // int velocity_decimate,
				velocity_radius);  // int velocity_radius)
		kernel5d =     consolidationKernel.getKernel();
		kernel5d_rad = consolidationKernel.getSpatialRadius();
		indices = getIndices(
				kernel5d_rad,
				width,
				height);
		indx_all_5d =     indices[0];
		indx_center_5d =  indices[1];
		indx_margins_5d = indices[2];
		// Direct pixel-domain kernel: same trapezoid projection over scalar frames,  // By Claude on 06/11/2026
		// one frame longer to span the same temporal extent as num_hist 3d3 pair-steps
		double [] tw_direct = new double [this.temporal_weights.length + 1];
		for (int h = 0; h < tw_direct.length; h++) {
			tw_direct[h] = this.temporal_weights[Math.min(h, this.temporal_weights.length - 1)];
		}
		normalize1d (tw_direct);
		consolidationKernelDirect = c5_matched ? // By Claude on 06/12/2026
			TemporalKernelGenerator.generateKernelMatched(
				tw_direct,         // double[] temporal_weights,
				pixel_decimate,    // int pixel_decimate,
				velocity_decimate, // int velocity_decimate,
				velocity_radius)   // int velocity_radius)
			:
			TemporalKernelGenerator.generateKernelDirect(
				tw_direct,         // double[] temporal_weights,
				pixel_decimate,    // int pixel_decimate,
				velocity_decimate, // int velocity_decimate,
				velocity_radius);  // int velocity_radius)
		kernel5dp_rad = consolidationKernelDirect.getSpatialRadius();
		indices = getIndices(
				kernel5dp_rad,
				width,
				height);
		indx_all_5dp =     indices[0];
		indx_center_5dp =  indices[1];
		indx_margins_5dp = indices[2];
		indices = getIndices(
				this.vel_suppr_rad,
				width,
				height);
		indx_all_vsupp =     indices[0];
		indx_center_vsupp =  indices[1];
		indx_margins_vsupp = indices[2];
	}
	
	public double getPsfRadius() {
		return psf_radius;
	}
	
	public int getKernel2dRadius() {
		return kernel2d_rad;
	}
	
	public int get3d3Radius() {
		return kernel3d3_rad;
	}

	public int getNumHist() {
		return kernel5d.length;
	}

	/** Temporal taps (frames) of the direct pixel-domain fine-velocity kernel. */
	public int getNumHistDirect() { // By Claude on 06/11/2026
		return consolidationKernelDirect.getKernel().length;
	}

	public double [] getKernel2d() {
		return kernel2d;
	}
	
	private void normalize1d(double [] data) {
		double s = 0;
		for (double d : data) s+= d;
		for (int i = 0; i < data.length; i++) {
			data[i] /= s;
		}
	}
			
	
	private static int [][] getIndices(
			int kernel_rad,
			int width,
			int height) {
		int kernel_size = 2 * kernel_rad+1;
		int [] i_all =     new int [width*height];
		int cent_width = width-kernel_size+1;
		int cent_height = height-kernel_size+1;
		int cent_size = cent_width*cent_height;
		int [] i_center =  new int [cent_size];
		int [] i_margins = new int [i_all.length - i_center.length];
		int n;
		for (n = 0; (1 << n) < i_all.length; n++);
		int shft = 32-n;
		int n_all = 1 << n;
//		int mask = n_all-1; 
		int iall = 0, icent=0, imarg=0;
		for (int i = 0; i < n_all; i++) {
			int ri = Integer.reverse(i)>>>shft;
			if (ri < i_all.length) {
				int y = ri/width;
				int x = ri%width;
				i_all[iall++] = ri;
				if ((x < kernel_rad) || (y < kernel_rad) || (x >= (width-kernel_rad)) ||  (y >= (height-kernel_rad))) {
					i_margins[imarg++] = ri;
				} else {
					i_center[icent++] = ri;
				}
			}
		}
		int [][] indices = {i_all, i_center, i_margins};
		return indices;
	}
	
	
	public void checkData(double [] data) {
		if (data.length != (width * height)) {
			throw new IllegalArgumentException("Data size = "+data.length+" does not match dimensions "+width+"*"+height+"="+(width*height));
		}
	}

	public void checkData(double [][] data) {
		if (data.length != (width * height)) {
			throw new IllegalArgumentException("Data size = "+data.length+" does not match dimensions "+width+"*"+height+"="+(width*height));
		}
	}
	
	public void checkKernel(double [] kernel, int rad) {
		int kernel_size = 2 * rad + 1;
		if (kernel.length != (kernel_size * kernel_size)) {
			throw new IllegalArgumentException("Kernel length = "+kernel.length+" does not match kernel size (2*"+kernel2d_rad+"+1) * (2*"+kernel2d_rad+"+1) ="+(kernel_size * kernel_size));
		}
	}
	
	public double [][][] convolve3d(
			final double [][][] data,         // outer index - historic layers, [0] - latest, second index - pixel, third - direction ([9])
			double [][][]     result_in,
			final double      min_frac_last){
		final int vel_size = 2*velocity_radius + 1;
		final int nlayers = data.length;
		final double [][][] result = (result_in != null) ? result_in : new double [width*height][pixel_decimate*pixel_decimate][vel_size * vel_size];
		final double [][][] last_layer = new double [width*height][pixel_decimate*pixel_decimate][vel_size * vel_size];
		final double [][][][][] kernel = consolidationKernel.getKernel();
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final int [] coarse_vel_indx = consolidationKernel.getCoarseVIdx();
		final int decimate_size = pixel_decimate*pixel_decimate;
		System.out.print("convolve3d() nhist=");
		for (int nhist = 0; nhist < nlayers; nhist++) {
			System.out.print(nhist+" ");
			final int fhist = nhist;
			final int spat_rad = consolidationKernel.getSpatialRadius(); // may be later dependent on nhist for efficiency 
			final double [][][][] kernel_hist = kernel[fhist];
			final double [][]     data_hist = data[fhist];
			ai.set(0);
			// processing center area, no need to care about margins
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_center_5d.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_center_5d[nPix];
							int target_x = ipix_dst % width; 
							int target_y = ipix_dst / width;
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
									int src_coarse_v = coarse_vel_indx[v_out_idx];
									if (src_coarse_v >=0) {
										double accumulated_flux = 0.0;
										// Outer spatial neighborhood limits defined by the pre-computed kernel
										for (int dy = -spat_rad; dy <= spat_rad; dy++) {
											int src_y = target_y + dy;
											for (int dx = -spat_rad; dx <= spat_rad; dx++) {
												int src_x = target_x + dx;
												// 1. Boundary check for source image coordinates - not here
												// 2. Locate the source macro-pixel index in your history array
												int src_pixel_idx = src_y * width + src_x;
												// 3. Pull historical energy from ONLY the mapped coarse velocity channel
												double historical_value = data_hist[src_pixel_idx][src_coarse_v];
												// 4. Multiply by the exact pre-calculated geometric weight for this sub-pixel layout
												int k_dy = dy + spat_rad;
												int k_dx = dx + spat_rad;
												double weight = kernel_hist[v_out_idx][sub_idx][k_dy][k_dx];
												accumulated_flux += historical_value * weight;
											}
										}
										result[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
										if (fhist == 0) {
											last_layer[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
										}
									}
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
			ai.set(0);
			// processing center area, no need to care about margins
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_margins_5d.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_margins_5d[nPix];
							int target_x = ipix_dst % width; 
							int target_y = ipix_dst / width;
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
									int src_coarse_v = coarse_vel_indx[v_out_idx];
									if (src_coarse_v >=0) {
										double accumulated_flux = 0.0;
										// Outer spatial neighborhood limits defined by the pre-computed kernel
										for (int dy = -spat_rad; dy <= spat_rad; dy++) {
											int src_y = target_y + dy;
											for (int dx = -spat_rad; dx <= spat_rad; dx++) {
												int src_x = target_x + dx;
												// 1. Boundary check for source image coordinates
												if (src_x >= 0 && src_x < width && src_y >= 0 && src_y < height) {
													// 2. Locate the source macro-pixel index in your history array
													int src_pixel_idx = src_y * width + src_x;
													// 3. Pull historical energy from ONLY the mapped coarse velocity channel
													double historical_value = data_hist[src_pixel_idx][src_coarse_v];
													// 4. Multiply by the exact pre-calculated geometric weight for this sub-pixel layout
													int k_dy = dy + spat_rad;
													int k_dx = dx + spat_rad;
													double weight = kernel_hist[v_out_idx][sub_idx][k_dy][k_dx];
													accumulated_flux += historical_value * weight;
												}
											}
										}
										result[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
										if (fhist == 0) {
											last_layer[ipix_dst][sub_idx][v_out_idx] += accumulated_flux;
										}
									}
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		// Remove if last_layer < fraction of the result
		if (min_frac_last > 0) {
			final double [][][][] kernel_hist = kernel[0];
			ai.set(0);
			// processing center area, no need to care about margins
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_all_5d.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_all_5d[nPix];
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
									if (last_layer[ipix_dst][sub_idx][v_out_idx] < min_frac_last * result[ipix_dst][sub_idx][v_out_idx]) {
										result[ipix_dst][sub_idx][v_out_idx] = 0;
									}
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		
		System.out.println();
		return result;
	}

	/**
	 * ROI variant: compute fine-velocity convolution for a pixel rectangle only.
	 * Result is indexed ROI-relatively: result[roi_pix][sub][vel] where
	 * roi_pix = (target_y - roi.y) * roi.width + (target_x - roi.x).
	 * Pixels outside roi are skipped; full-image source data is still used for
	 * the spatial neighbourhood reads.
	 */
	public double [][][] convolve3d(
			final double [][][] data,
			final Rectangle     roi,
			double [][][]       result_in,
			final double        min_frac_last){
		final int vel_size = 2*velocity_radius + 1;
		final int nlayers = data.length;
		final int roi_npix = roi.width * roi.height;
		final double [][][] result = (result_in != null) ? result_in :
				new double [roi_npix][pixel_decimate*pixel_decimate][vel_size * vel_size];
		final double [][][] last_layer = new double [roi_npix][pixel_decimate*pixel_decimate][vel_size * vel_size];
		final double [][][][][] kernel = consolidationKernel.getKernel();
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final int [] coarse_vel_indx = consolidationKernel.getCoarseVIdx();
		final int decimate_size = pixel_decimate*pixel_decimate;
		final int roi_x0 = roi.x, roi_y0 = roi.y, roi_x1 = roi.x + roi.width, roi_y1 = roi.y + roi.height;
		for (int nhist = 0; nhist < nlayers; nhist++) {
			final int fhist = nhist;
			final int spat_rad = consolidationKernel.getSpatialRadius();
			final double [][][][] kernel_hist = kernel[fhist];
			final double [][]     data_hist = data[fhist];
			ai.set(0);
			// center pixels — no source boundary checks needed
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_center_5d.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_center_5d[nPix];
							int target_x = ipix_dst % width;
							int target_y = ipix_dst / width;
							if (target_x < roi_x0 || target_x >= roi_x1 ||
								target_y < roi_y0 || target_y >= roi_y1) continue;
							int roi_pix = (target_y - roi_y0) * roi.width + (target_x - roi_x0);
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
									int src_coarse_v = coarse_vel_indx[v_out_idx];
									if (src_coarse_v >= 0) {
										double accumulated_flux = 0.0;
										for (int dy = -spat_rad; dy <= spat_rad; dy++) {
											int src_y = target_y + dy;
											for (int dx = -spat_rad; dx <= spat_rad; dx++) {
												int src_x = target_x + dx;
												int src_pixel_idx = src_y * width + src_x;
												double historical_value = data_hist[src_pixel_idx][src_coarse_v];
												int k_dy = dy + spat_rad;
												int k_dx = dx + spat_rad;
												double weight = kernel_hist[v_out_idx][sub_idx][k_dy][k_dx];
												accumulated_flux += historical_value * weight;
											}
										}
										result[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
										if (fhist == 0) {
											last_layer[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
										}
										
									}
								}
							}
						}
					}
				};
			}
			ImageDtt.startAndJoin(threads);
			ai.set(0);
			// margin pixels — source boundary checks required
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_margins_5d.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_margins_5d[nPix];
							int target_x = ipix_dst % width;
							int target_y = ipix_dst / width;
							if (target_x < roi_x0 || target_x >= roi_x1 ||
								target_y < roi_y0 || target_y >= roi_y1) continue;
							int roi_pix = (target_y - roi_y0) * roi.width + (target_x - roi_x0);
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
									int src_coarse_v = coarse_vel_indx[v_out_idx];
									if (src_coarse_v >= 0) {
										double accumulated_flux = 0.0;
										for (int dy = -spat_rad; dy <= spat_rad; dy++) {
											int src_y = target_y + dy;
											for (int dx = -spat_rad; dx <= spat_rad; dx++) {
												int src_x = target_x + dx;
												if (src_x >= 0 && src_x < width && src_y >= 0 && src_y < height) {
													int src_pixel_idx = src_y * width + src_x;
													double historical_value = data_hist[src_pixel_idx][src_coarse_v];
													int k_dy = dy + spat_rad;
													int k_dx = dx + spat_rad;
													double weight = kernel_hist[v_out_idx][sub_idx][k_dy][k_dx];
													accumulated_flux += historical_value * weight;
												}
											}
										}
										result[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
										if (fhist == 0) {
											last_layer[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
										}
									}
								}
							}
						}
					}
				};
			}
			ImageDtt.startAndJoin(threads);
		}
		// Remove if last_layer < fraction of the result
		if (min_frac_last > 0) {
			final double [][][][] kernel_hist = kernel[0];
			ai.set(0);
			// processing center area, no need to care about margins
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_all_5d.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_all_5d[nPix];
							int target_x = ipix_dst % width;
							int target_y = ipix_dst / width;
							if (target_x < roi_x0 || target_x >= roi_x1 ||
								target_y < roi_y0 || target_y >= roi_y1) continue;
							int roi_pix = (target_y - roi_y0) * roi.width + (target_x - roi_x0);
							
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < kernel_hist.length; v_out_idx++) {
									if (last_layer[roi_pix][sub_idx][v_out_idx] < min_frac_last * result[roi_pix][sub_idx][v_out_idx]) {
										result[roi_pix][sub_idx][v_out_idx] = 0;
									}
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		return result;
	}

	/**
	 * Direct pixel-domain fine-velocity convolution (ROI only) - bypasses the     // By Claude on 06/11/2026
	 * 3x3 coarse-velocity stage: a single-stage velocity-matched filter applied
	 * straight to the LoG+LReLU frames. Optimal in the low-SNR linear regime,
	 * where convolve3D3LReLU adds ~2x background noise at equal contrast.
	 * Output shape and indexing identical to convolve3d(roi):
	 * result[roi_pix][sub][vel], roi_pix = (y-roi.y)*roi.width + (x-roi.x).
	 * @param data          [h][npix] scalar frames, [0] - NEWEST, getNumHistDirect() entries
	 * @param roi           target pixel rectangle (full-image coordinates)
	 * @param result_in     reuse array or null
	 * @param min_frac_last zero outputs where the newest frame contributes less
	 *                      than this fraction (dark-frame feedforward gate)
	 */
	public double [][][] convolve3dDirect( // By Claude on 06/11/2026
			final double [][]   data,
			final Rectangle     roi,
			double [][][]       result_in,
			final double        min_frac_last){
		final int vel_size = 2*velocity_radius + 1;
		final int nlayers = data.length;
		final int roi_npix = roi.width * roi.height;
		final double [][][] result = (result_in != null) ? result_in :
				new double [roi_npix][pixel_decimate*pixel_decimate][vel_size * vel_size];
		final double [][][] last_layer = new double [roi_npix][pixel_decimate*pixel_decimate][vel_size * vel_size];
		final double [][][][][] kernel = consolidationKernelDirect.getKernel();
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final int decimate_size = pixel_decimate*pixel_decimate;
		final int roi_x0 = roi.x, roi_y0 = roi.y, roi_x1 = roi.x + roi.width, roi_y1 = roi.y + roi.height;
		// Sparse tap lists: the (matched) kernel is mostly zeros - per (t, v, sub) the   // By Claude on 06/12/2026
		// bump trajectory covers only a few cells of the (2*rad+1)^2 square. Iterate
		// only nonzero taps (~20-30x fewer multiply-accumulates than the dense loops).
		final int k_rad = consolidationKernelDirect.getSpatialRadius();
		final int [][][][]    tap_dx = new int    [kernel.length][][][];
		final int [][][][]    tap_dy = new int    [kernel.length][][][];
		final double [][][][] tap_wt = new double [kernel.length][][][];
		for (int t = 0; t < kernel.length; t++) {
			int nvout_t = kernel[t].length;
			tap_dx[t] = new int [nvout_t][][];
			tap_dy[t] = new int [nvout_t][][];
			tap_wt[t] = new double [nvout_t][][];
			for (int v = 0; v < nvout_t; v++) {
				tap_dx[t][v] = new int [decimate_size][];
				tap_dy[t][v] = new int [decimate_size][];
				tap_wt[t][v] = new double [decimate_size][];
				for (int sub = 0; sub < decimate_size; sub++) {
					double [][] kw = kernel[t][v][sub];
					int n = 0;
					for (int ky = 0; ky < kw.length; ky++) for (int kx = 0; kx < kw[ky].length; kx++) {
						if (kw[ky][kx] != 0) n++;
					}
					int [] tdx = new int [n];
					int [] tdy = new int [n];
					double [] tw = new double [n];
					n = 0;
					for (int ky = 0; ky < kw.length; ky++) for (int kx = 0; kx < kw[ky].length; kx++) {
						if (kw[ky][kx] != 0) {
							tdy[n] = ky - k_rad;
							tdx[n] = kx - k_rad;
							tw[n++] = kw[ky][kx];
						}
					}
					tap_dx[t][v][sub] = tdx;
					tap_dy[t][v][sub] = tdy;
					tap_wt[t][v][sub] = tw;
				}
			}
		}
		for (int nhist = 0; nhist < nlayers; nhist++) {
			final int fhist = nhist;
			final int [][][]    tap_dx_h = tap_dx[fhist];  // By Claude on 06/12/2026
			final int [][][]    tap_dy_h = tap_dy[fhist];
			final double [][][] tap_wt_h = tap_wt[fhist];
			final double []       data_hist = data[fhist];
			ai.set(0);
			// center pixels - no source boundary checks needed
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_center_5dp.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_center_5dp[nPix];
							int target_x = ipix_dst % width;
							int target_y = ipix_dst / width;
							if (target_x < roi_x0 || target_x >= roi_x1 ||
								target_y < roi_y0 || target_y >= roi_y1) continue;
							int roi_pix = (target_y - roi_y0) * roi.width + (target_x - roi_x0);
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < tap_dx_h.length; v_out_idx++) {
									final int []    tdx = tap_dx_h[v_out_idx][sub_idx]; // sparse taps // By Claude on 06/12/2026
									final int []    tdy = tap_dy_h[v_out_idx][sub_idx];
									final double [] tw =  tap_wt_h[v_out_idx][sub_idx];
									double accumulated_flux = 0.0;
									for (int k = 0; k < tdx.length; k++) {
										accumulated_flux += data_hist[(target_y + tdy[k]) * width + target_x + tdx[k]] * tw[k];
									}
									result[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
									if (fhist == 0) {
										last_layer[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
									}
								}
							}
						}
					}
				};
			}
			ImageDtt.startAndJoin(threads);
			ai.set(0);
			// margin pixels - source boundary checks required
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_margins_5dp.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_margins_5dp[nPix];
							int target_x = ipix_dst % width;
							int target_y = ipix_dst / width;
							if (target_x < roi_x0 || target_x >= roi_x1 ||
								target_y < roi_y0 || target_y >= roi_y1) continue;
							int roi_pix = (target_y - roi_y0) * roi.width + (target_x - roi_x0);
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < tap_dx_h.length; v_out_idx++) {
									final int []    tdx = tap_dx_h[v_out_idx][sub_idx]; // sparse taps // By Claude on 06/12/2026
									final int []    tdy = tap_dy_h[v_out_idx][sub_idx];
									final double [] tw =  tap_wt_h[v_out_idx][sub_idx];
									double accumulated_flux = 0.0;
									for (int k = 0; k < tdx.length; k++) {
										int src_x = target_x + tdx[k];
										int src_y = target_y + tdy[k];
										if (src_x >= 0 && src_x < width && src_y >= 0 && src_y < height) {
											accumulated_flux += data_hist[src_y * width + src_x] * tw[k];
										}
									}
									result[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
									if (fhist == 0) {
										last_layer[roi_pix][sub_idx][v_out_idx] += accumulated_flux;
									}
								}
							}
						}
					}
				};
			}
			ImageDtt.startAndJoin(threads);
		}
		// Remove if last_layer < fraction of the result
		if (min_frac_last > 0) {
			final int num_vout = kernel[0].length;
			ai.set(0);
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPix = ai.getAndIncrement(); nPix < indx_all_5dp.length; nPix = ai.getAndIncrement()) {
							int ipix_dst = indx_all_5dp[nPix];
							int target_x = ipix_dst % width;
							int target_y = ipix_dst / width;
							if (target_x < roi_x0 || target_x >= roi_x1 ||
								target_y < roi_y0 || target_y >= roi_y1) continue;
							int roi_pix = (target_y - roi_y0) * roi.width + (target_x - roi_x0);
							for (int sub_idx = 0; sub_idx < decimate_size; sub_idx++) {
								for (int v_out_idx = 0; v_out_idx < num_vout; v_out_idx++) {
									if (last_layer[roi_pix][sub_idx][v_out_idx] < min_frac_last * result[roi_pix][sub_idx][v_out_idx]) {
										result[roi_pix][sub_idx][v_out_idx] = 0;
									}
								}
							}
						}
					}
				};
			}
			ImageDtt.startAndJoin(threads);
		}
		return result;
	}

	public double[][] convolve3D3LReLU(
			final double [] data,
			final double [] data_prev,
			double [][]     result_in,
			final int       mode_3d3,
			final double    curt_3d3_frac,  //  = 0.5; // for stronger signals give this fraction of the total (of 9) weight to the 1-of-9 winner
			final double    curt_3d3_thrsh, //  = 1.0; // scale down curt_3d3_frac for weaker signals (in ratio of min(data[now]/curt_3d3_thrsh,1.0)
			final double    curt_3d3_power, //  = 2.0; // in mode_3d3=2, redistribute coarse velocities weights as the specified power of their positive values (ignore negative ones)
			final double    w_now, // normally 0.5, so center will propagate to the next pyramid levels
			final double    w_prev,// normally 0.5, so center will propagate to the next pyramid levels
			final double    alpha,
			final int       debugLevel) {
		switch (mode_3d3) {
		case CuasDetectRT.CONV3D3_MODE_PLAIN:
			return convolve3D3LReLU(
				data,         // final double [] data,
				data_prev,    // final double [] data_prev,
				result_in,   // double [][]     result_in,
				w_now,       // final double    w_now, // normally 0.5, so center will propagate to the next pyramid levels
				w_prev,      // final double    w_prev,// normally 0.5, so center will propagate to the next pyramid levels
				alpha,       // final double    alpha,
				debugLevel); // final int       debugLevel) 
		case CuasDetectRT.CONV3D3_MODE_WINNER:
			return convolve3D3LReLU(
					data,           // final double [] data,
					data_prev,      // final double [] data_prev,
					result_in,      // double [][]     result_in,
					curt_3d3_frac,  // final double    curt_3d3_frac, //  =        0.5;   // for stronger signals give this fraction of the total (of 9) weight to the 1-of-9 winner
					curt_3d3_thrsh, //final double    curt_3d3_thrsh, //  =       1.0;   // scale down curt_3d3_frac for weaker signals (in ratio of min(data[now]/curt_3d3_thrsh,1.0)
					w_now,          // final double    w_now, // normally 0.5, so center will propagate to the next pyramid levels
					w_prev,         // final double    w_prev,// normally 0.5, so center will propagate to the next pyramid levels
					alpha,          // final double    alpha,
					debugLevel);    // final int       debugLevel) 
		case CuasDetectRT.CONV3D3_MODE_POWER:
			return convolve3D3LReLU(
					data,           // final double [] data,
					data_prev,      // final double [] data_prev,
					result_in,      // double [][]     result_in,
					curt_3d3_frac,  // final double    curt_3d3_frac, //  =        0.5;   // for stronger signals give this fraction of the total (of 9) weight to the 1-of-9 winner
					curt_3d3_thrsh, // final double    curt_3d3_thrsh, //  =       1.0;   // scale down curt_3d3_frac for weaker signals (in ratio of min(data[now]/curt_3d3_thrsh,1.0)
					curt_3d3_power, // final double    curt_3d3_power, //  = 2.0; // in mode_3d3=2, redistribute coarse velocities weights as the specified power of their positive values (ignore negative ones)
					w_now,          // final double    w_now, // normally 0.5, so center will propagate to the next pyramid levels
					w_prev,         // final double    w_prev,// normally 0.5, so center will propagate to the next pyramid levels
					alpha,          // final double    alpha,
					debugLevel);    // final int       debugLevel) 
		default:
			System.out.println ("convolve3D3LReLU(): mode_3d3="+mode_3d3+" is undefined");
			return null;
			
		}
	}	
	
	
	/**
	 * 
	 * @param data
	 * @param data_prev
	 * @param result_in
	 * @param curt_3d3_frac
	 * @param curt_3d3_thrsh
	 * @param w_now
	 * @param w_prev
	 * @param alpha
	 * @param debugLevel
	 * @return
	 */
	public double[][] convolve3D3LReLU(
			final double [] data,
			final double [] data_prev,
			double [][]     result_in,
			final double   curt_3d3_frac, //  =        0.5;   // for stronger signals give this fraction of the total (of 9) weight to the 1-of-9 winner
			final double   curt_3d3_thrsh, //  =       1.0;   // scale down curt_3d3_frac for weaker signals (in ratio of min(data[now]/curt_3d3_thrsh,1.0)
			final double    w_now, // normally 0.5, so center will propagate to the next pyramid levels
			final double    w_prev,// normally 0.5, so center will propagate to the next pyramid levels
			final double    alpha,
			final int       debugLevel) {
		final int kern_size_3d3 = 2 * kernel3d3_rad + 1;
		final double [][] result = (result_in != null) ? result_in : new double [data.length][kern_size_3d3*kern_size_3d3];
		int dbg_thresh = 10; // -4;
		checkData(data);
		checkData(data_prev);
		checkData(result);
		final int tl_offset =  kernel3d3_rad * (width+1);
		final int center_index = 2 * kernel3d3_rad * (kernel3d3_rad + 1);
		final int row_offset = width - (2*kernel3d3_rad +1);
		final int dbg_pix = (debugLevel > dbg_thresh) ? (300+ 173*width):-1;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		// processing center area, no need to care about margins
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					double [] coarse_vel = new double [kern_size_3d3*kern_size_3d3];
					for (int nPix = ai.getAndIncrement(); nPix < indx_center_3d3.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_center_3d3[nPix];
						if (ipix_dst == dbg_pix) {
							System.out.println("convolve3D3LReLU(): ipix_dst="+ipix_dst);	
						}
						int ipix_prev = ipix_dst- tl_offset;
						int ikern = kern_size_3d3*kern_size_3d3;
						Arrays.fill(coarse_vel, 0);
						double v0 = w_now * data[ipix_dst]; // will be added to all velocity outputs
						for (int iy = 0; iy < kern_size_3d3; iy++) {
							for (int ix = 0; ix < kern_size_3d3; ix++) {
								double v = data_prev[ipix_prev++] * w_prev + v0;
								coarse_vel[--ikern] = v;
							}
							ipix_prev += row_offset;
						}
						// velocity competition
						double frac = Math.min(curt_3d3_frac, data[ipix_dst]/curt_3d3_thrsh);
						if (frac > 0) {
							double sumv = 0;
							int imx = 0;
							double bestd=coarse_vel[imx];
							for (int i = 0; i < coarse_vel.length;i++) {
								double d = coarse_vel[i]; 
								sumv += d;
								coarse_vel[i] = d * (1.0 -frac);
								if (d > bestd) {
									imx = i;
									bestd= d;
								}
							}
							coarse_vel[imx] += sumv * frac;
						}
						// Apply LReLU
						for (int i = 0; i < coarse_vel.length;i++) {
							double v = coarse_vel[i]; 
							result[ipix_dst][i] = Math.max(v, alpha * v);		
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		// marginal area, check indices
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					TileNeibs tn = new TileNeibs(width, height);
					double [] coarse_vel = new double [kern_size_3d3*kern_size_3d3];
					for (int nPix = ai.getAndIncrement(); nPix < indx_margins_3d3.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_margins_3d3[nPix];
						int ikern = kern_size_3d3*kern_size_3d3;
						Arrays.fill(coarse_vel, 0);
						double v0 = w_now * data[ipix_dst]; // will be added to all velocity outputs
						for (int iy = 0; iy < kern_size_3d3; iy++) {
							int dy = iy-kernel3d3_rad;
							for (int ix = 0; ix < kern_size_3d3; ix++) {
								int dx = ix-kernel3d3_rad;
								int ipix_prev = tn.getNeibIndex(ipix_dst, dx, dy);
								double v = v0;
								if (ipix_prev >= 0) {
									v += data_prev[ipix_prev] * w_prev;
								}
								coarse_vel[--ikern] = v;
							}
						}
						// velocity competition
						double frac = Math.min(curt_3d3_frac, data[ipix_dst]/curt_3d3_thrsh);
						if (frac > 0) {
							double sumv = 0;
							int imx = 0;
							double bestd=coarse_vel[imx];
							for (int i = 0; i < coarse_vel.length;i++) {
								double d = coarse_vel[i]; 
								sumv += d;
								coarse_vel[i] = d * (1.0 -frac);
								if (d > bestd) {
									imx = i;
									bestd= d;
								}
							}
							coarse_vel[imx] += sumv * frac;
						}
						// Apply LReLU
						for (int i = 0; i < coarse_vel.length;i++) {
							double v = coarse_vel[i]; 
							result[ipix_dst][i] = Math.max(v, alpha * v);		
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return result;
	}
	
	
	
	public double[][] convolve3D3LReLU(
			final double [] data,
			final double [] data_prev,
			double [][]     result_in,
			final double    curt_3d3_frac, //  =        0.5;   // for stronger signals give this fraction of the total (of 9) weight to the 1-of-9 winner
			final double    curt_3d3_thrsh, //  =       1.0;   // scale down curt_3d3_frac for weaker signals (in ratio of min(data[now]/curt_3d3_thrsh,1.0)
			final double    curt_3d3_power, //  = 2.0; // in mode_3d3=2, redistribute coarse velocities weights as the specified power of their positive values (ignore negative ones)
			final double    w_now, // normally 0.5, so center will propagate to the next pyramid levels
			final double    w_prev,// normally 0.5, so center will propagate to the next pyramid levels
			final double    alpha,
			final int       debugLevel) {
		final int kern_size_3d3 = 2 * kernel3d3_rad + 1;
		final double [][] result = (result_in != null) ? result_in : new double [data.length][kern_size_3d3*kern_size_3d3];
		int dbg_thresh = 10; // -4;
		checkData(data);
		checkData(data_prev);
		checkData(result);
		final int tl_offset =  kernel3d3_rad * (width+1);
//		final int center_index = 2 * kernel3d3_rad * (kernel3d3_rad + 1);
		final int row_offset = width - (2*kernel3d3_rad +1);
		final int dbg_pix = (debugLevel > dbg_thresh) ? (300+ 173*width):-1;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		// processing center area, no need to care about margins
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					double [] coarse_vel =     new double [kern_size_3d3*kern_size_3d3];
					double [] coarse_vel_pwr = new double [kern_size_3d3*kern_size_3d3];
					for (int nPix = ai.getAndIncrement(); nPix < indx_center_3d3.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_center_3d3[nPix];
						if (ipix_dst == dbg_pix) {
							System.out.println("convolve3D3LReLU(): ipix_dst="+ipix_dst);	
						}
						int ipix_prev = ipix_dst- tl_offset;
						int ikern = kern_size_3d3*kern_size_3d3;
						Arrays.fill(coarse_vel, 0);
						Arrays.fill(coarse_vel_pwr, 0);
						double v0 = w_now * data[ipix_dst]; // will be added to all velocity outputs
						for (int iy = 0; iy < kern_size_3d3; iy++) {
							for (int ix = 0; ix < kern_size_3d3; ix++) {
								double v = data_prev[ipix_prev++] * w_prev + v0;
								coarse_vel[--ikern] = v;
							}
							ipix_prev += row_offset;
						}
						// velocity competition
						double frac = Math.min(curt_3d3_frac, data[ipix_dst]/curt_3d3_thrsh);
						if (frac > 0) {
							double sumv = 0;
							double sumv_pwr = 0;
							for (int i = 0; i < coarse_vel.length;i++) {
								double d = coarse_vel[i]; // includes v0;
								double d_pwr = (d > 0) ? Math.pow(d, curt_3d3_power) : 0;
								sumv += d;
								sumv_pwr += d_pwr; 
								coarse_vel_pwr[i] = d_pwr; 
							}
							double w = sumv * frac / sumv_pwr;
							if (sumv_pwr > 0) {
								for (int i = 0; i < coarse_vel.length;i++) {
									coarse_vel[i] = (1.0 -frac) * coarse_vel[i]  + w * coarse_vel_pwr[i]; 
								}								
							}
						}
						// Apply LReLU
						for (int i = 0; i < coarse_vel.length;i++) {
							double v = coarse_vel[i]; 
							result[ipix_dst][i] = Math.max(v, alpha * v);		
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		// marginal area, check indices
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					TileNeibs tn = new TileNeibs(width, height);
					double [] coarse_vel = new double [kern_size_3d3*kern_size_3d3];
					double [] coarse_vel_pwr = new double [kern_size_3d3*kern_size_3d3];
					for (int nPix = ai.getAndIncrement(); nPix < indx_margins_3d3.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_margins_3d3[nPix];
						int ikern = kern_size_3d3*kern_size_3d3;
						Arrays.fill(coarse_vel, 0);
						Arrays.fill(coarse_vel_pwr, 0);
						double v0 = w_now * data[ipix_dst]; // will be added to all velocity outputs
						for (int iy = 0; iy < kern_size_3d3; iy++) {
							int dy = iy-kernel3d3_rad;
							for (int ix = 0; ix < kern_size_3d3; ix++) {
								int dx = ix-kernel3d3_rad;
								int ipix_prev = tn.getNeibIndex(ipix_dst, dx, dy);
								double v = v0;
								if (ipix_prev >= 0) {
									v += data_prev[ipix_prev] * w_prev;
								}
								coarse_vel[--ikern] = v;
							}
						}
						// velocity competition
						double frac = Math.min(curt_3d3_frac, data[ipix_dst]/curt_3d3_thrsh);
						if (frac > 0) {
							double sumv = 0;
							double sumv_pwr = 0;
							for (int i = 0; i < coarse_vel.length;i++) {
								double d = coarse_vel[i]; // includes v0;
								double d_pwr = (d > 0) ? Math.pow(d, curt_3d3_power) : 0;
								sumv += d;
								sumv_pwr += d_pwr; 
								coarse_vel_pwr[i] = d_pwr; 
							}
							double w = sumv * frac / sumv_pwr;
							if (sumv_pwr > 0) {
								for (int i = 0; i < coarse_vel.length;i++) {
									coarse_vel[i] = (1.0 -frac) * coarse_vel[i]  + w * coarse_vel_pwr[i]; 
								}								
							}
						}
						// Apply LReLU
						for (int i = 0; i < coarse_vel.length;i++) {
							double v = coarse_vel[i]; 
							result[ipix_dst][i] = Math.max(v, alpha * v);		
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return result;
	}
	
	
	
	public double[][] convolve3D3LReLU(
			final double [] data,
			final double [] data_prev,
			double [][]     result_in,
			final double    w_now, // normally 0.5, so center will propagate to the next pyramid levels
			final double    w_prev,// normally 0.5, so center will propagate to the next pyramid levels
			final double    alpha,
			final int       debugLevel) {
		final int kern_size_3d3 = 2 * kernel3d3_rad + 1;
		final double [][] result = (result_in != null) ? result_in : new double [data.length][kern_size_3d3*kern_size_3d3];
		int dbg_thresh = 10; // -4;
		checkData(data);
		checkData(data_prev);
		checkData(result);
		final int tl_offset =  kernel3d3_rad * (width+1);
		final int row_offset = width - (2*kernel3d3_rad +1);
		final int dbg_pix = (debugLevel > dbg_thresh) ? (300+ 173*width):-1;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		// processing center area, no need to care about margins
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < indx_center_3d3.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_center_3d3[nPix];
						if (ipix_dst == dbg_pix) {
							System.out.println("convolve3D3LReLU(): ipix_dst="+ipix_dst);	
						}
						int ipix_prev = ipix_dst- tl_offset;
						int ikern = kern_size_3d3*kern_size_3d3;
						double v0 = w_now * data[ipix_dst]; // will ve added to all velocity outputs
						for (int iy = 0; iy < kern_size_3d3; iy++) {
							for (int ix = 0; ix < kern_size_3d3; ix++) {
								double v = data_prev[ipix_prev++] * w_prev + v0;
								result[ipix_dst][--ikern] = Math.max(v, alpha * v);								
							}
							ipix_prev += row_offset;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		// marginal area, check indices
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					TileNeibs tn = new TileNeibs(width, height);
					for (int nPix = ai.getAndIncrement(); nPix < indx_margins_3d3.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_margins_3d3[nPix];
						int ikern = kern_size_3d3*kern_size_3d3;
						double v0 = w_now * data[ipix_dst]; // will be added to all velocity outputs
						for (int iy = 0; iy < kern_size_3d3; iy++) {
							int dy = iy-kernel3d3_rad;
							for (int ix = 0; ix < kern_size_3d3; ix++) {
								int dx = ix-kernel3d3_rad;
								int ipix_prev = tn.getNeibIndex(ipix_dst, dx, dy);
								double v = v0;
								if (ipix_prev >= 0) {
									v += data_prev[ipix_prev] * w_prev;
								}
								result[ipix_dst][--ikern] = Math.max(v, alpha * v);								
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return result;
	}
	
	/**
	 * Combine (pairwise average)2 pixel arrays 		
	 * @param data_now  pixel data [width*height]
	 * @param data_prev pixel data [width*height]
	 * @param result_in allocated pixel array or null (to be allocated)
	 * @param alpha for LReLU (0 - classic ReLU, 1.0 - linear transfer, no ReLU)
	 * @return result pixel array, same as result_in if not null
	 */
	public double[] temporalAverageLReLU(
			final double [] data_now,
			final double [] data_prev,
			double []       result_in,
			final double    alpha) {
		// assuming symmetrical kernel, not inverting indices
		final double [] result = (result_in != null) ? result_in : new double [data_now.length];
		final int kernel_size = 2 * kernel2d_rad + 1;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		// processing center area, no need to care about margins
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < indx_all_2d.length; nPix = ai.getAndIncrement()) {
						int ipix = indx_all_2d[nPix];
						double v = 0.5*(data_now[ipix]+data_prev[ipix]);
						result[ipix] = Math.max(v, alpha * v);
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return result;
	}
	
	
	
	public double[] convolve2DLReLU(
			final double [] data,
			final double [] kernel, // square
			double []       result_in,
			final double    alpha) {
		// assuming symmetrical kernel, not inverting indices
		final double [] result = (result_in != null) ? result_in : new double [data.length];
		// verify dimensions;
		checkData(data);
		checkData(result);
		checkKernel(kernel, kernel2d_rad);
		final int tl_offset = kernel2d_rad * (width+1);
		final int kernel_size = 2 * kernel2d_rad + 1;
		final int row_inc = width- kernel_size;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		// processing center area, no need to care about margins
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPix = ai.getAndIncrement(); nPix < indx_center_2d.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_center_2d[nPix];
						int ipix_src = ipix_dst- tl_offset;
						int ikern = 0;
						double v = 0;
						for (int iy = 0; iy < kernel_size; iy++) {
							for (int ix = 0; ix < kernel_size; ix++) {
								v+=data[ipix_src++]*kernel[ikern++];
							}
							ipix_src+=row_inc;
						}
						result[ipix_dst] = Math.max(v, alpha * v);
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		// marginal area, check indices
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					TileNeibs tn = new TileNeibs(width, height);
					for (int nPix = ai.getAndIncrement(); nPix < indx_margins_2d.length; nPix = ai.getAndIncrement()) {
						int ipix_dst = indx_margins_2d[nPix];
						double v = 0;
						for (int iy = 0; iy < kernel_size; iy++) {
							int dy = iy-kernel2d_rad;
							for (int ix = 0; ix < kernel_size; ix++) {
								int dx = ix-kernel2d_rad;
								int ipix_src = tn.getNeibIndex(ipix_dst, dx, dy);
								if (ipix_src >= 0) {
									v+=data[ipix_src]*kernel[ix + iy*kernel_size];	
								}
							}
						}
						result[ipix_dst] = Math.max(v, alpha * v);
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return result;
	}
	
    /**
     * Calculates the integer radius for the LoG kernel based on the PSF radius
     * and the desired truncation boundary (n_sigma).
     */
    public static int getRadiusLoGKernel(double psf_radius, double n_sigma) {
        double sigma = psf_radius / Math.sqrt(2.0);
        return (int) Math.ceil(sigma * n_sigma);
    }

    /**
     * Generates a 1D flattened flux-normalized 2D Laplacian of Gaussian kernel.
     * The sum of positive elements is exactly 1.0, and the sum of negative elements is exactly -1.0.
     */
    public static double[] getLoGKenel(double psf_radius, int kernel_radius) {
        double sigma = psf_radius / Math.sqrt(2.0);
        double sigma2 = sigma * sigma;
       
        int width = 2 * kernel_radius + 1;
        double[] kernel = new double[width * width];
       
        double sumPos = 0.0;
        double sumNeg = 0.0;
       
        // 1. Generate raw inverted LoG shape and accumulate sums
        for (int y = -kernel_radius; y <= kernel_radius; y++) {
            for (int x = -kernel_radius; x <= kernel_radius; x++) {
                double r2 = x * x + y * y;
               
                // Positive center, negative surround for bright blob detection
                double raw = (1.0 - (r2 / (2.0 * sigma2))) * Math.exp(-r2 / (2.0 * sigma2));
               
                int idx = (y + kernel_radius) * width + (x + kernel_radius);
                kernel[idx] = raw;
               
                if (raw > 0) {
                    sumPos += raw;
                } else {
                    sumNeg += raw;
                }
            }
        }
       
        // 2. Apply strict flux-normalization
        double absSumNeg = Math.abs(sumNeg);
        for (int i = 0; i < kernel.length; i++) {
            if (kernel[i] > 0) {
                kernel[i] /= sumPos;
            } else if (kernel[i] < 0) {
                kernel[i] /= absSumNeg;
            }
        }
       
        return kernel;
    }

// debug methods
    public ImagePlus showKernel2d(String title) {
    	if (title==null) {
    		title = "kernel2d";
    	}
    	ImagePlus imp_kernel2d = ShowDoubleFloatArrays.makeArrays(
    			kernel2d,
				2 * kernel2d_rad + 1,
				2 * kernel2d_rad + 1,
				title);
    	imp_kernel2d.setDisplayRange(IMP_MINMAX[0],IMP_MINMAX[1]);
    	return imp_kernel2d;
    }
    
    public double [] render3d3Flat(double [][] coarse_vel) {
    	int width41 =  4 * width+1;
    	int height41 = 4 * height+1;
    	double [] render_flat = new double [width41 * height41];
    	Arrays.fill(render_flat, Double.NaN);
    	for (int npix=0; npix < coarse_vel.length; npix++) if (coarse_vel[npix] !=null){
    		int py = npix/width;
    		int px = npix%width;
    		int pix4 = (px * 4 + 1) + (py * 4 + 1) * width41;
    		int cvi = 0;
    		for (int dy = 0; dy < 3; dy++) {
        		for (int dx = 0; dx < 3; dx++) {
        			render_flat[pix4++] = coarse_vel[npix][cvi++];
        		}
        		pix4 += width41 - 3; // next row and 3 pixels back
    		}
    	}
    	return render_flat;
    }

    public double [][] render3d3Sliced(double [][] coarse_vel) {
    	double [][] render_sliced = new double [9][width * height];
    	for (double [] slice:render_sliced) {
    		Arrays.fill(slice, Double.NaN);
    	}
    	for (int npix=0; npix < coarse_vel.length; npix++) if (coarse_vel[npix] !=null){
    		for (int s = 0; s < coarse_vel[npix].length; s++) {
    			render_sliced[s][npix] = coarse_vel[npix][s];
    		}
    	}
    	return render_sliced;
    }
    
    public ImagePlus render3d3Hyperstack(
    		double [][][] coarse_vel_seq, // [nser][npix][nvel]
    		String []     scene_names,
    		String        title) {
    	String [] arrows9 = get9Arrows();
    	double [][][] render_sliced = new double [arrows9.length][coarse_vel_seq.length][width * height];
    	for (double [][] slices:render_sliced) {
    		for (double [] slice:slices) {
    			Arrays.fill(slice, Double.NaN);
    		}
    	}
    	for (int nseq = 0; nseq < coarse_vel_seq.length; nseq++) {
    		for (int npix=0; npix < coarse_vel_seq[nseq].length; npix++) if (coarse_vel_seq[nseq][npix] !=null){
    			for (int s = 0; s < coarse_vel_seq[nseq][npix].length; s++) {
    				render_sliced[s][nseq][npix] = coarse_vel_seq[nseq][npix][s];
    			}
    		}
    	}
		ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
				render_sliced, // double[][][] pixels, 
				width,         // int          width, 
				title,         // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				scene_names,   // all slices*frames titles or just slice titles or null
				arrows9,       // String []    frame_titles, // frame titles or null
				false);        // boolean      show)
		imp.setDisplayRange(IMP_MINMAX[0],IMP_MINMAX[1]);
    	return imp;
    }
    
    
    public String [] get9Arrows() {
    	String [] arrows= {
    			"\u2191\u2190",  // Up-Left   
    			"\u2191",        // Up        
    			"\u2192\u2191",  // Right-Up  
    			"\u2190",        // Left      
    			"\u25E6",        // Bullet    
    			"\u2192",        // Right     
    			"\u2193\u2190",  // Down-Left 
    			"\u2193",        // Down
    			"\u2192\u2193"}; // Right-Down
    	return arrows;
    }
    public String [] get9Arrows_diags() {
    	String [] arrows= {
    			"\u2196",  // Left-Up   
    			"\u2191",  // Up        
    			"\u2197",  // Right-Up  
    			"\u2190",  // Left      
    			"\u25E6",  // Bullet    
    			"\u2192",  // Right     
    			"\u2199",  // Left-Down
    			"\u2193",  // Down
    			"\u2198"};  // Right-Down
    	return arrows;
    }
    
    public ImagePlus showKernel5d(String title) {
    	// consolidationKernel 
    	// [num_subs][num_temps][pixels]
    	if (title == null) {
    		title = "kernel5d";
    	}
    	int spat_dim = consolidationKernel.getSpatialDim();
    	int coarse_dim = 3;
    	int coarse_num = coarse_dim*coarse_dim;
    	int render_dim = spat_dim * (coarse_dim + 1);
    	int num_hist = kernel5d.length;
    	int num_vout = kernel5d[0].length;
    	int vout_dim = (int) Math.sqrt(num_vout);
    	int vout_rad = (vout_dim-1)/2;
    	
    	int num_sub =  kernel5d[0][0].length;
    	int sub_dim = (int) Math.sqrt(num_sub);
    	
    	double [][][] rendered = new double [num_sub*num_hist][num_vout][(render_dim + 1) * (render_dim + 1)];
    	String [] top_titles =   new String[num_sub*num_hist];
    	String [] slice_titles = new String[num_vout];
    	for (double [][] slice_hist:rendered) {
    		for (double [] slice:slice_hist) {
    			Arrays.fill(slice, Double.NaN);
    		}
    	}
    	for (int nvout = 0; nvout < num_vout; nvout++) {
    		int vout_v = nvout / vout_dim - vout_rad ; 
    		int vout_h = nvout % vout_dim - vout_rad;
    		slice_titles[nvout] = nvout+ " ("+((vout_v >0) ?"+":"")+vout_v+":"+((vout_h >0) ?"+":"")+vout_h+")";
    	}
    	for (int nsub = 0; nsub < num_sub; nsub++) {
    		int sub_v = nsub / sub_dim; 
    		int sub_h = nsub % sub_dim;
        	for (int nhist = 0; nhist < num_hist; nhist++) {
        		int num_tslice = nsub * num_hist + nhist;
        		top_titles[num_tslice] = sub_v+":"+sub_h+" - " + nhist;
            	for (int nvout = 0; nvout < num_vout; nvout++) {
            		int coarse_index = consolidationKernel.getCoarseVIdx()[nvout];
            		for (int py = 0; py < spat_dim; py++) {
            			int py0_rend = py * (coarse_dim + 1) +1; // ) * render_dim;
            			for (int px = 0; px < spat_dim; px++) {
            				int px0_rend = px * (coarse_dim + 1) +1; //	
            				for (int coarse_i = 0; coarse_i < coarse_num; coarse_i++) {
            					double w = 0; // maybe opposite - use NaN here?
            					int coarse_v = coarse_i / coarse_dim; 
            					int coarse_h = coarse_i % coarse_dim;
            					if (coarse_i == coarse_index) {
            						w = kernel5d[nhist][nvout][nsub][py][px]; // unused will be 0, not NaN;
            					}
            					int pix_rend = (py0_rend + coarse_v) * (render_dim+1) + (px0_rend + coarse_h);
            					rendered[num_tslice][nvout][pix_rend] = w;
            				}
            			}
            		}
            	}        		
        	}    		
    	}
		ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
				rendered,      // double[][][] pixels, 
				(render_dim+1),// int          width, 
				title,         // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				slice_titles,  // String []    titles, // all slices*frames titles or just slice titles or null
				top_titles,    // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
				false);        // boolean      show)
		imp.setDisplayRange(IMP_MINMAX[0],IMP_MINMAX[1]);
    	return imp;
    }
   /*
    public ImagePlus showKernel5d_1(String title) {
    	// consolidationKernel 
    	// [num_subs][num_temps][pixels]
    	if (title == null) {
    		title = "kernel5d";
    	}
    	int spat_dim = consolidationKernel.getSpatialDim();
    	int coarse_dim = 3;
    	int coarse_num = coarse_dim*coarse_dim;
    	int render_dim = spat_dim * (coarse_dim+1);
    	int num_hist = kernel5d.length;
    	int num_vout = kernel5d[0].length;
    	int vout_dim = (int) Math.sqrt(num_vout);
    	int vout_rad = (vout_dim-1)/2;
    	
    	int num_sub =  kernel5d[0][0].length;
    	int sub_dim = (int) Math.sqrt(num_sub);
    	
//    	double [][][] rendered = new double [num_sub*num_hist][num_vout][(render_dim + 1) * (render_dim + 1)];
    	double [][][] rendered = new double [num_hist][num_vout][(render_dim + 1) * (render_dim + 1)];
    	String [] top_titles =   new String[num_sub*num_hist];
    	String [] slice_titles = new String[num_vout];
    	for (double [][] slice_hist:rendered) {
    		for (double [] slice:slice_hist) {
    			Arrays.fill(slice, Double.NaN);
    		}
    	}
    	for (int nvout = 0; nvout < num_vout; nvout++) {
    		int vout_v = nvout / vout_dim - vout_rad ; 
    		int vout_h = nvout % vout_dim - vout_rad;
    		slice_titles[nvout] = nvout+ " ("+((vout_v >0) ?"+":"")+vout_v+":"+((vout_h >0) ?"+":"")+vout_h+")";
    	}
//    	for (int nsub = 0; nsub < num_sub; nsub++) {
//    		int sub_v = nsub / sub_dim; 
//   		int sub_h = nsub % sub_dim;
        	for (int nhist = 0; nhist < num_hist; nhist++) {
        		int num_tslice = nsub * num_hist + nhist;
        		top_titles[num_tslice] = sub_v+":"+sub_h+" - " + nhist;
            	for (int nvout = 0; nvout < num_vout; nvout++) {
            		int coarse_index = consolidationKernel.getCoarseVIdx()[nvout];
            		for (int py = 0; py < spat_dim; py++) {
            			int py0_rend = py * (coarse_dim + 1) +1; // ) * render_dim;
            			for (int px = 0; px < spat_dim; px++) {
            				int px0_rend = px * (coarse_dim + 1) +1; //	
            				for (int coarse_i = 0; coarse_i < coarse_num; coarse_i++) {
            					double w = 0; // maybe opposite - use NaN here?
            					int coarse_v = coarse_i / coarse_dim; 
            					int coarse_h = coarse_i % coarse_dim;
            					if (coarse_i == coarse_index) {
            						w = kernel5d[nhist][nvout][nsub][py][px]; // unused will be 0, not NaN;
            					}
            					int pix_rend = (py0_rend + coarse_v) * render_dim + (px0_rend + coarse_h);
            					rendered[num_tslice][nvout][pix_rend] = w;
            				}
            			}
            		}
            	}        		
        	}    		
    	
		ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
				rendered,      // double[][][] pixels, 
				width,         // int          width, 
				title,         // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				slice_titles,  // String []    titles, // all slices*frames titles or just slice titles or null
				top_titles,    // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
				false);        // boolean      show)
		imp.setDisplayRange(IMP_MINMAX[0],IMP_MINMAX[1]);
    	return imp;
    }
*/    
    
    
    // Display convolution result [vout][2*width*2*height]. Or hyperstack with x and y sliders?
    /**
     * Display a sequence of velocities per subpixel (or pixel if no decimation) layers as
     * a hyperstack. The bottom slider (outer index) selects the output velocity, the top slider
     * (and a scroll wheel) - scenes (or averaged scenes named by the last integrated) 
     * @param data [nscene][npix][nsubpix][nvelocity]
     * @param scene_titles scene names or null
     * @param title  image title or null
     * @return
     */
    public ImagePlus showConvKernel5d(
    		double [][][][] data,
    		String []       scene_titles,
    		String          title) {
    	int num_scenes = data.length;
    	int num_vout = kernel5d[0].length;
    	int vout_dim = (int) Math.sqrt(num_vout);
    	int vout_rad = (vout_dim-1)/2;
    	int num_sub =  kernel5d[0][0].length;
    	int sub_dim = (int) Math.sqrt(num_sub);
    	int full_with =   width *  sub_dim;
    	int full_height = height * sub_dim;
    	if (title==null) {
    		title = "OutputConvKernel5D";
    	}
    	if (scene_titles == null) {
    		scene_titles = new String[num_scenes];
    		for (int i = 0; i < scene_titles.length; i++) {
    			scene_titles[i]= String.format("scene-%03d", i);
    		}
    	}
		String [] v_titles = new String[num_vout];

    	for (int nvout = 0; nvout < num_vout; nvout++) {
    		int vout_v = nvout / vout_dim - vout_rad ; 
    		int vout_h = nvout % vout_dim - vout_rad;
    		v_titles[nvout] = nvout+ " ("+((vout_v >0) ?"+":"")+vout_v+":"+((vout_h >0) ?"+":"")+vout_h+")";
    	}
    	
    	double [][][] rendered = new double [num_vout][num_scenes][full_with*full_height];
    	for (double [][] d2 : rendered) for (double [] d1 : d2) Arrays.fill(d1, Double.NaN);
	
		ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
				rendered,      // double[][][] pixels, 
				width,         // int          width, 
				title,         // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
				scene_titles,  // slice_titles,  // String []    titles, // all slices*frames titles or just slice titles or null
				v_titles,      // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
				false);        // boolean      show)
		imp.setDisplayRange(IMP_MINMAX[0],IMP_MINMAX[1]);
    	return imp;
    }

    /**
     * Render fine-velocity output for a pixel ROI as a flat image stack (one slice per scene).
     * For each pixel in the ROI a (vout_dim × vout_dim) velocity-map square is drawn.
     * Squares are tiled in the output image with a 1-pixel NaN separator column/row
     * (renders as visible grid in ImageJ). Sub-pixels (pixel_decimate > 1) tile inside each square.
     *
     * Output dimensions:
     *   width  = roi.width  * (vout_dim + 1) * sub_dim
     *   height = roi.height * (vout_dim + 1) * sub_dim
     *
     * @param data         [nscene][width*height][nsub][nvel] — output of convolve3d
     * @param roi          pixel rectangle to render (full-image coordinates)
     * @param scene_titles per-slice labels, or null
     * @param title        image title, or null
     */
    public ImagePlus showConvKernel5d(
    		double [][][][] data,
    		Rectangle       roi,
    		String []       scene_titles,
    		String          title) {
    	int num_scenes = data.length;
    	int num_vout   = kernel5d[0].length;
    	int vout_dim   = (int) Math.sqrt(num_vout);
    	int num_sub    = kernel5d[0][0].length;
    	int sub_dim    = (int) Math.sqrt(num_sub);

    	int roi_width  = roi.width  * (vout_dim + 1) * sub_dim;
    	int roi_height = roi.height * (vout_dim + 1) * sub_dim;
    	int block_w    = (vout_dim + 1) * sub_dim;
    	int block_h    = (vout_dim + 1) * sub_dim;

    	if (title == null) title = "CONV5D-ROI";
    	if (scene_titles == null) {
    		scene_titles = new String[num_scenes];
    		for (int i = 0; i < scene_titles.length; i++) {
    			scene_titles[i] = String.format("scene-%03d", i);
    		}
    	}

    	double [][] rendered = new double[num_scenes][roi_width * roi_height];
    	for (double [] slice : rendered) Arrays.fill(slice, Double.NaN);

    	for (int nscene = 0; nscene < num_scenes; nscene++) {
    		for (int py_roi = 0; py_roi < roi.height; py_roi++) {
    			int py_img = roi.y + py_roi;
    			if (py_img < 0 || py_img >= height) continue;
    			for (int px_roi = 0; px_roi < roi.width; px_roi++) {
    				int px_img = roi.x + px_roi;
    				if (px_img < 0 || px_img >= width) continue;
    				int pixel_idx = py_roi * roi.width + px_roi;
    				int block_x   = px_roi * block_w;
    				int block_y   = py_roi * block_h;
    				double [][]   pix_data = data[nscene][pixel_idx]; // [nsub][nvel]
    				for (int sy = 0; sy < sub_dim; sy++) {
    					for (int sx = 0; sx < sub_dim; sx++) {
    						int sub_idx = sy * sub_dim + sx;
    						int slot_x  = block_x + sx * (vout_dim + 1);
    						int slot_y  = block_y + sy * (vout_dim + 1);
    						for (int vy = 0; vy < vout_dim; vy++) {
    							for (int vx = 0; vx < vout_dim; vx++) {
    								int v_out_idx = vy * vout_dim + vx;
    								int rx = slot_x + vx;
    								int ry = slot_y + vy;
    								rendered[nscene][ry * roi_width + rx] = pix_data[sub_idx][v_out_idx];
    							}
    						}
    					}
    				}
    			}
    		}
    	}

    	ImagePlus imp = ShowDoubleFloatArrays.makeArrays(
    			rendered,
    			roi_width,
    			roi_height,
    			title,
    			scene_titles);
    	imp.setDisplayRange(IMP_MINMAX[0], IMP_MINMAX[1]);
    	return imp;
    }

    /**
     * Render fine-velocity ROI output as an ImageJ hyperstack.
     * Bottom slider (T) iterates over all num_vout velocity outputs, labelled "+vy:+vx".
     * Second slider (Z) iterates over scenes.
     * Image size: roi.width * pixel_decimate × roi.height * pixel_decimate — no gap grid.
     * Sub-pixels are spatially embedded: pixel (px,py), sub-pixel (sx,sy) → image pixel
     * (px*sub_dim+sx, py*sub_dim+sy).
     *
     * @param data         [nscene][roi_npix][nsub][nvel] — ROI-indexed output of convolve3d(roi)
     * @param roi          pixel rectangle (full-image coordinates)
     * @param scene_titles per-scene labels, or null
     * @param title        image title, or null
     */
    public ImagePlus showConvKernel5dHyperRect(
    		double [][][][] data,
    		Rectangle       roi,
    		String []       scene_titles,
    		String          title) {
    	int num_scenes = data.length;
    	int num_vout   = kernel5d[0].length;
    	int vout_dim   = (int) Math.sqrt(num_vout);
    	int vout_rad   = (vout_dim - 1) / 2;
    	int num_sub    = kernel5d[0][0].length;
    	int sub_dim    = (int) Math.sqrt(num_sub);

    	int img_width  = roi.width  * sub_dim;
    	int img_height = roi.height * sub_dim;

    	if (title == null) title = "CONV5D-HYPER-RECT";
    	if (scene_titles == null) {
    		scene_titles = new String[num_scenes];
    		for (int i = 0; i < scene_titles.length; i++) {
    			scene_titles[i] = String.format("scene-%03d", i);
    		}
    	}
    	// Extra leading velocity slice (index 0) = max over all velocities, a quick locator    // By Claude on 06/13/2026
    	// of non-empty pixels. It is the first slice ImageJ shows, and the display range below
    	// is scaled to it so the small (~0.1) maxima are immediately visible on open.
    	String[] v_titles = new String[num_vout + 1];
    	v_titles[0] = "MAX (all v)";
    	for (int nvout = 0; nvout < num_vout; nvout++) {
    		int vout_v = nvout / vout_dim - vout_rad;
    		int vout_h = nvout % vout_dim - vout_rad;
    		v_titles[nvout + 1] = nvout + " (" + ((vout_v > 0) ? "+" : "") + vout_v +
    				":" + ((vout_h > 0) ? "+" : "") + vout_h + ")";
    	}

    	// outer dim = velocities (T = bottom slider; index 0 = MAX-over-v), inner dim = scenes (Z = second slider)
    	double [][][] pixels_hyper = new double [num_vout + 1][num_scenes][img_width * img_height];
    	for (double [][] d2 : pixels_hyper) for (double [] d1 : d2) Arrays.fill(d1, Double.NaN);

    	for (int nscene = 0; nscene < num_scenes; nscene++) {
    		for (int py_roi = 0; py_roi < roi.height; py_roi++) {
    			int py_img = roi.y + py_roi;
    			if (py_img < 0 || py_img >= height) continue;
    			for (int px_roi = 0; px_roi < roi.width; px_roi++) {
    				int px_img = roi.x + px_roi;
    				if (px_img < 0 || px_img >= width) continue;
    				int roi_pix = py_roi * roi.width + px_roi;
    				double [][] pix_data = data[nscene][roi_pix]; // [nsub][nvel]
    				for (int sy = 0; sy < sub_dim; sy++) {
    					for (int sx = 0; sx < sub_dim; sx++) {
    						int sub_idx = sy * sub_dim + sx;
    						int ix      = px_roi * sub_dim + sx;
    						int iy      = py_roi * sub_dim + sy;
    						int img_pix = iy * img_width + ix;
    						double vmax = Double.NEGATIVE_INFINITY;                              // By Claude on 06/13/2026
    						for (int v_out_idx = 0; v_out_idx < num_vout; v_out_idx++) {
    							double v = pix_data[sub_idx][v_out_idx];
    							pixels_hyper[v_out_idx + 1][nscene][img_pix] = v;                 // velocities shifted by 1
    							if (v > vmax) vmax = v;
    						}
    						pixels_hyper[0][nscene][img_pix] = vmax;                              // leading MAX-over-v slice
    					}
    				}
    			}
    		}
    	}

    	// Display range from the MAX-over-v slice (the one that opens first), so non-empty      // By Claude on 06/13/2026
    	// pixels are auto-visible; velocity slices share it (their values are <= the max).
    	double disp_min = Double.POSITIVE_INFINITY, disp_max = Double.NEGATIVE_INFINITY;
    	for (double [] sc : pixels_hyper[0]) for (double v : sc) if (!Double.isNaN(v)) {
    		if (v < disp_min) disp_min = v;
    		if (v > disp_max) disp_max = v;
    	}
    	if (!(disp_max > disp_min)) { disp_min = IMP_MINMAX[0]; disp_max = IMP_MINMAX[1]; }       // fallback (empty/flat)

    	ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
    			pixels_hyper, // double[num_vout+1][num_scenes][img_w*img_h]
    			img_width,
    			title,
    			scene_titles, // Z slider (second)
    			v_titles,     // T slider (bottom)
    			false);
    	imp.setDisplayRange(disp_min, disp_max);
    	return imp;
    }

    /**
     * 4D non-maximum suppression on a ROI fine-velocity result.
     * For each (pixel, velocity) point, find the maximum M over the ±vel_suppr_rad neighbourhood
     * in both pixel space and velocity space (vx/vy), then:
     *   thresh = max(athresh, M*rthresh); suppress if (M - val) >= thresh
     *   otherwise keep val unchanged (including the weak-signal passthrough when M < athresh).
     *
     * Uses indx_center_vsupp / indx_margins_vsupp for bit-reversed multithreaded iteration
     * (kernel radius = vel_suppr_rad matches both pixel and velocity neighbourhood).
     *
     * @param data     [roi_npix][nsub][nvel] — ROI-indexed output of convolve3d(roi)
     * @param roi      pixel rectangle (full-image coordinates)
     * @param athresh  absolute threshold: M must reach this for suppression to activate
     * @param rthresh  relative threshold: val/M below this → zeroed (e.g. 0.5)
     * @param dbg_min  debug feature - zero if weaker to test recurrent with strong targets
     * @return new filtered array, same shape as data
     */
    public double [][][] filterConv5dROI(
    		final double [][][] data,
    		final Rectangle     roi,
    		final double        athresh,
    		final double        rthresh,
    		final double        dbg_min) { // zero if below, for testing with strong signals
    	final int num_vout = kernel5d[0].length;
    	final int vout_dim = (int) Math.sqrt(num_vout);
    	final int num_sub  = kernel5d[0][0].length;
    	final int roi_x0 = roi.x, roi_y0 = roi.y;
    	final int roi_x1 = roi.x + roi.width, roi_y1 = roi.y + roi.height;
    	final double [][][] result = new double [data.length][num_sub][num_vout];
    	final Thread[] threads = ImageDtt.newThreadArray();
    	final AtomicInteger ai = new AtomicInteger(0);

    	// Center pixels: all 8 pixel neighbours are within image bounds; only ROI check needed.
    	for (int ithread = 0; ithread < threads.length; ithread++) {
    		threads[ithread] = new Thread() { public void run() {
    			for (int nPix = ai.getAndIncrement(); nPix < indx_center_vsupp.length; nPix = ai.getAndIncrement()) {
    				int ipix = indx_center_vsupp[nPix];
    				int px = ipix % width, py = ipix / width;
    				if (px < roi_x0 || px >= roi_x1 || py < roi_y0 || py >= roi_y1) continue;
    				int roi_pix = (py - roi_y0) * roi.width + (px - roi_x0);
    				for (int sub_idx = 0; sub_idx < num_sub; sub_idx++) {
    					for (int v_out_idx = 0; v_out_idx < num_vout; v_out_idx++) {
    						double val = data[roi_pix][sub_idx][v_out_idx];
    						int vy = v_out_idx / vout_dim, vx = v_out_idx % vout_dim;
    						double M = 0;
    						for (int dpy = -vel_suppr_rad; dpy <= vel_suppr_rad; dpy++) {
    							int ny = py + dpy;
    							for (int dpx = -vel_suppr_rad; dpx <= vel_suppr_rad; dpx++) {
    								int nx = px + dpx;
    								if (nx < roi_x0 || nx >= roi_x1 || ny < roi_y0 || ny >= roi_y1) continue;
    								int nroi = (ny - roi_y0) * roi.width + (nx - roi_x0);
    								for (int dvy = -vel_suppr_rad; dvy <= vel_suppr_rad; dvy++) {
    									int nvy = vy + dvy;
    									if (nvy < 0 || nvy >= vout_dim) continue;
    									for (int dvx = -vel_suppr_rad; dvx <= vel_suppr_rad; dvx++) {
    										int nvx = vx + dvx;
    										if (nvx < 0 || nvx >= vout_dim) continue;
    										double d = data[nroi][sub_idx][nvy * vout_dim + nvx];
    										if (d > M) M = d;
    									}
    								}
    							}
    						}
    						double thresh = Math.max(athresh, M * rthresh);
    						result[roi_pix][sub_idx][v_out_idx] =((val < dbg_min) || ((M - val) >= thresh)) ? 0 : val;
    					}
    				}
    			}
    		}};
    	}
    	ImageDtt.startAndJoin(threads);

    	// Margin pixels: additionally check image bounds before ROI bounds.
    	ai.set(0);
    	for (int ithread = 0; ithread < threads.length; ithread++) {
    		threads[ithread] = new Thread() { public void run() {
    			for (int nPix = ai.getAndIncrement(); nPix < indx_margins_vsupp.length; nPix = ai.getAndIncrement()) {
    				int ipix = indx_margins_vsupp[nPix];
    				int px = ipix % width, py = ipix / width;
    				if (px < roi_x0 || px >= roi_x1 || py < roi_y0 || py >= roi_y1) continue;
    				int roi_pix = (py - roi_y0) * roi.width + (px - roi_x0);
    				for (int sub_idx = 0; sub_idx < num_sub; sub_idx++) {
    					for (int v_out_idx = 0; v_out_idx < num_vout; v_out_idx++) {
    						double val = data[roi_pix][sub_idx][v_out_idx];
    						int vy = v_out_idx / vout_dim, vx = v_out_idx % vout_dim;
    						double M = 0;
    						for (int dpy = -vel_suppr_rad; dpy <= vel_suppr_rad; dpy++) {
    							int ny = py + dpy;
    							if (ny < 0 || ny >= height) continue;
    							for (int dpx = -vel_suppr_rad; dpx <= vel_suppr_rad; dpx++) {
    								int nx = px + dpx;
    								if (nx < 0 || nx >= width) continue;
    								if (nx < roi_x0 || nx >= roi_x1 || ny < roi_y0 || ny >= roi_y1) continue;
    								int nroi = (ny - roi_y0) * roi.width + (nx - roi_x0);
    								for (int dvy = -vel_suppr_rad; dvy <= vel_suppr_rad; dvy++) {
    									int nvy = vy + dvy;
    									if (nvy < 0 || nvy >= vout_dim) continue;
    									for (int dvx = -vel_suppr_rad; dvx <= vel_suppr_rad; dvx++) {
    										int nvx = vx + dvx;
    										if (nvx < 0 || nvx >= vout_dim) continue;
    										double d = data[nroi][sub_idx][nvy * vout_dim + nvx];
    										if (d > M) M = d;
    									}
    								}
    							}
    						}
    						double thresh = Math.max(athresh, M * rthresh);
    						result[roi_pix][sub_idx][v_out_idx] =((val < dbg_min) || ((M - val) >= thresh)) ? 0 : val;
    					}
    				}
    			}
    		}};
    	}
    	ImageDtt.startAndJoin(threads);
    	return result;
    }

    // === C5P Bayesian posterior (softmax over the matched-filter bank) ===              // By Claude on 06/12/2026
    // Exact single-segment posterior: P(theta|D) = exp(z^2/2 + ln pi) / (Lambda0 + sum), //
    // where the matched filter output is the sufficient statistic, z^2/2 the GLRT        //
    // exponent, pi(v) the velocity prior, and the shared denominator (competition over   //
    // the local 4D hypothesis window + the no-target hypothesis) is what enforces the    //
    // single-target prior. Hard NMS (filterConv5dROI) is the lambda->0 limit.            //

    private double [][] kernelDirectEnergy = null; // lazy [nsub][nvout] sum of squared taps // By Claude on 06/12/2026

    /**
     * Per-channel template energy of the direct (C5P) kernel: E[sub][v] = sum_{t,dy,dx} w^2. // By Claude on 06/12/2026
     * For white input noise sigma the convolve3dDirect output variance is sigma^2*E, so
     * the normalized SNR of a channel is z = out / (sigma*sqrt(E)).
     */
    public double [][] getKernelDirectEnergy() { // By Claude on 06/12/2026
    	if (kernelDirectEnergy == null) {
    		double [][][][][] k = consolidationKernelDirect.getKernel();
    		int num_vout = k[0].length;
    		int num_sub =  k[0][0].length;
    		double [][] e = new double [num_sub][num_vout];
    		for (int t = 0; t < k.length; t++) {
    			for (int v = 0; v < num_vout; v++) {
    				for (int sub = 0; sub < num_sub; sub++) {
    					double s = 0;
    					for (double [] row : k[t][v][sub]) for (double w : row) s += w*w;
    					e[sub][v] += s;
    				}
    			}
    		}
    		kernelDirectEnergy = e;
    	}
    	return kernelDirectEnergy;
    }

    /**
     * Radially-symmetric log-prior ln(pi) over the fine-velocity grid: 1.0 inside r0     // By Claude on 06/12/2026
     * (physical max speed in velocity cells), cos^2 taper r0..r1 (cosine bump wing in
     * velocity space), ~0 beyond (physically impossible). Radial - the speed bound is
     * isotropic, so the square-grid corners (e.g. (4,4)=1.41 px/frame) are excluded.
     * Falling prior matches the temporal-pyramid handoff: the same target appears first
     * at lower levels where its pixel velocity is lower.
     */
    public double [] getVelLogPrior(double r0, double r1) { // By Claude on 06/12/2026
    	int v_dim = 2 * velocity_radius + 1;
    	double [] log_prior = new double [v_dim * v_dim];
    	for (int vy = -velocity_radius; vy <= velocity_radius; vy++) {
    		for (int vx = -velocity_radius; vx <= velocity_radius; vx++) {
    			double r = Math.sqrt(vx*vx + vy*vy);
    			double pi;
    			if (r <= r0)      pi = 1.0;
    			else if (r >= r1) pi = 0.0;
    			else {
    				double c = Math.cos(0.5 * Math.PI * (r - r0) / (r1 - r0));
    				pi = c * c;
    			}
    			log_prior[(vy + velocity_radius) * v_dim + (vx + velocity_radius)] =
    					Math.log(Math.max(pi, 1e-12));
    		}
    	}
    	return log_prior;
    }

    /**
     * Per-cell log-likelihood exponent ("score") of the single-segment Bayes model:      // By Claude on 06/12/2026
     *   score = z^2/2 + ln pi(v),   z = max(out,0) / (lambda*sqrt(E[sub][v]))
     * out is the (amplitude-calibrated) convolve3dDirect output, lambda the assumed
     * white-noise sigma of the input frames - the "fat zero" sweep parameter: small
     * lambda -> sharp confident posterior (and amplified noise), large lambda -> blur
     * toward the prior. One-sided z (amplitude prior s>=0) is the principled LReLU.
     * @param data      [roi_npix][nsub][nvel] - raw convolve3dDirect output (pre-NMS)
     * @param lambda    assumed input noise sigma
     * @param log_prior [nvel] per-velocity ln pi (getVelLogPrior()), or null for flat
     * @return same shape as data
     */
    /** log of the standard normal CDF, numerically stable for z << 0. */ // By Claude on 06/13/2026
    private static double lnPhi(double z) {
    	if (z > -10.0) {
    		return Math.log(0.5 * erfc(-z / Math.sqrt(2.0)));
    	}
    	// asymptotic: Phi(z) ~ phi(z)/(-z), ln Phi ~ -z^2/2 - ln(-z) - 0.5 ln(2pi)
    	return -0.5 * z * z - Math.log(-z) - 0.5 * Math.log(2.0 * Math.PI);
    }

    /** complementary error function (Numerical Recipes erfcc, |err| < 1.2e-7). */ // By Claude on 06/13/2026
    private static double erfc(double x) {
    	double t = 1.0 / (1.0 + 0.5 * Math.abs(x));
    	double ans = t * Math.exp(-x * x - 1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 +
    			t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 +
    			t * (-0.82215223 + t * 0.17087277)))))))));
    	return (x >= 0.0) ? ans : 2.0 - ans;
    }

    /**
     * Per-cell log-likelihood "score" of the single-segment Bayes model.                 // By Claude on 06/12/2026
     * z = out / (lambda*sqrt(E)). exact=false: score = z^2/2 (z one-sided via max(out,0))
     * + ln pi - the quick GLRT profile. exact=true: the full strength(s>=0)-marginalized   // By Claude on 06/13/2026
     * log-posterior, score = z^2/2 + ln Phi(z) - 0.5 ln E + ln pi:
     *   - ln Phi(z): exact s>=0 one-sidedness (smooth suppression of negative/near-zero
     *     correlation, replacing the hard max);
     *   - -0.5 ln E: per-channel template-energy normalization (fair cross-velocity
     *     comparison; E varies with sub-pixel phase).
     * Note: -0.5 ln E is lambda-independent while z^2/2 ~ 1/lambda^2, so at large lambda the
     * energy term carries relatively more weight (correct: weak signal -> normalization
     * dominates). The constant 0.5 ln(2 pi sigma^2) cancels in the softmax.
     * @param exact full marginalization (true) vs the quick z^2/2 profile (false)
     */
    public double [][][] scoreC5PROI( // By Claude on 06/12/2026
    		final double [][][] data,
    		final double        lambda,
    		final double []     log_prior,
    		final boolean       exact) { // By Claude on 06/13/2026
    	final double [][] energy = getKernelDirectEnergy();
    	final int num_sub  = energy.length;
    	final int num_vout = energy[0].length;
    	final double [][][] result = new double [data.length][num_sub][num_vout];
    	final Thread[] threads = ImageDtt.newThreadArray();
    	final AtomicInteger ai = new AtomicInteger(0);
    	for (int ithread = 0; ithread < threads.length; ithread++) {
    		threads[ithread] = new Thread() { public void run() {
    			for (int nPix = ai.getAndIncrement(); nPix < data.length; nPix = ai.getAndIncrement()) {
    				for (int sub_idx = 0; sub_idx < num_sub; sub_idx++) {
    					for (int v = 0; v < num_vout; v++) {
    						double lp = (log_prior != null) ? log_prior[v] : 0.0;
    						double val = data[nPix][sub_idx][v];
    						double E = energy[sub_idx][v];
    						double sc;
    						if (E > 0) {
    							double z = val / (lambda * Math.sqrt(E));
    							if (exact) {
    								sc = 0.5 * z * z + lnPhi(z) - 0.5 * Math.log(E); // full s>=0 marginal // By Claude on 06/13/2026
    							} else {
    								sc = (val > 0) ? (0.5 * z * z) : 0.0; // quick GLRT profile (hard one-sided)
    							}
    						} else {
    							sc = exact ? Double.NEGATIVE_INFINITY : 0.0; // no template energy - impossible // By Claude on 06/13/2026
    						}
    						result[nPix][sub_idx][v] = sc + lp;
    					}
    				}
    			}
    		}};
    	}
    	ImageDtt.startAndJoin(threads);
    	return result;
    }

    /**
     * Competition (softmax) stage of the Bayes posterior: every (pixel,velocity)         // By Claude on 06/12/2026
     * hypothesis in the +/-vel_suppr_rad 4D window competes with its neighbors and the
     * no-target hypothesis for the same probability mass:
     *   P_c = exp(s_c - m) / (exp(s0 - m) + sum_window exp(s - m)),  m = window max
     * Soft version of filterConv5dROI 4D NMS (its lambda->0 limit is argmax). Sub-pixel
     * phases compete only within their own sub-grid (consistent with filterConv5dROI).
     * @param score      output of scoreC5PROI [roi_npix][nsub][nvel]
     * @param roi        pixel rectangle (full-image coordinates)
     * @param null_score s0 = null_z^2/2: a detection must beat a null_z-sigma fluctuation
     * @return posterior in [0,1], same shape as score
     */
    public double [][][] softmaxC5PROI( // By Claude on 06/12/2026; separable rewrite (same math, ~1000x fewer exp) // By Claude on 06/13/2026
    		final double [][][] score,
    		final Rectangle     roi,
    		final double        null_score) {
    	// Exact decomposition of the per-cell 4D logsumexp:
    	//   exp(s[np][nv] - K) = exp(s - M[np]) * exp(M[np] - K),  M[np] = per-pixel max over all velocities
    	// so the velocity part collapses to per-pixel velocity-box sums B (stage A, one exp
    	// per cell), and each output cell needs only 2*vsr+1 squared weighted adds (stage B)
    	// instead of (2*vsr+1)^4 exp calls. K = max(window M, null) keeps every exp <= 1.
    	final int num_vout = kernel5d[0].length;
    	final int vout_dim = (int) Math.sqrt(num_vout);
    	final int num_sub  = kernel5d[0][0].length;
    	final int rw = roi.width, rh = roi.height;
    	final int npix = score.length; // == rw*rh (ROI-indexed, as produced by convolve3dDirect/scoreC5PROI)
    	final int vsr = vel_suppr_rad;
    	final int win = 2 * vsr + 1;
    	final double [][][] result = new double [npix][num_sub][num_vout];
    	final double [][]   pix_max = new double [num_sub][npix];          // M: per-pixel max over velocities
    	final double [][][] vel_box = new double [npix][num_sub][num_vout]; // B: per-pixel velocity-box sums of exp(s-M)
    	final double [][][] vel_bmax = new double [npix][num_sub][num_vout]; // M2: per-pixel velocity-box MAX of s (for the underflow path) // By Claude on 06/13/2026
    	final Thread[] threads = ImageDtt.newThreadArray();
    	final AtomicInteger ai = new AtomicInteger(0);
    	// Stage A (independent per pixel): M, e = exp(s-M), B = clipped +-vsr velocity box sum of e,
    	// M2 = clipped +-vsr velocity box max of s (no exp - feeds the per-cell underflow path)
    	for (int ithread = 0; ithread < threads.length; ithread++) {
    		threads[ithread] = new Thread() { public void run() {
    			double [] e   = new double [num_vout];
    			double [] tmp = new double [num_vout];
    			for (int nPix = ai.getAndIncrement(); nPix < npix; nPix = ai.getAndIncrement()) {
    				for (int sub_idx = 0; sub_idx < num_sub; sub_idx++) {
    					double [] s = score[nPix][sub_idx];
    					double m = Double.NEGATIVE_INFINITY;
    					for (int v = 0; v < num_vout; v++) if (s[v] > m) m = s[v];
    					pix_max[sub_idx][nPix] = m;
    					for (int v = 0; v < num_vout; v++) e[v] = Math.exp(s[v] - m);
    					for (int vy = 0; vy < vout_dim; vy++) {      // 1D box along vx
    						int row = vy * vout_dim;
    						for (int vx = 0; vx < vout_dim; vx++) {
    							double bsum = 0;
    							int x0 = Math.max(0, vx - vsr), x1 = Math.min(vout_dim - 1, vx + vsr);
    							for (int nvx = x0; nvx <= x1; nvx++) bsum += e[row + nvx];
    							tmp[row + vx] = bsum;
    						}
    					}
    					double [] bx = vel_box[nPix][sub_idx];
    					for (int vx = 0; vx < vout_dim; vx++) {      // 1D box along vy
    						for (int vy = 0; vy < vout_dim; vy++) {
    							double bsum = 0;
    							int y0 = Math.max(0, vy - vsr), y1 = Math.min(vout_dim - 1, vy + vsr);
    							for (int nvy = y0; nvy <= y1; nvy++) bsum += tmp[nvy * vout_dim + vx];
    							bx[vy * vout_dim + vx] = bsum;
    						}
    					}
    					for (int vy = 0; vy < vout_dim; vy++) {      // 1D max along vx // By Claude on 06/13/2026
    						int row = vy * vout_dim;
    						for (int vx = 0; vx < vout_dim; vx++) {
    							double bmax = Double.NEGATIVE_INFINITY;
    							int x0 = Math.max(0, vx - vsr), x1 = Math.min(vout_dim - 1, vx + vsr);
    							for (int nvx = x0; nvx <= x1; nvx++) if (s[row + nvx] > bmax) bmax = s[row + nvx];
    							tmp[row + vx] = bmax;
    						}
    					}
    					double [] bm = vel_bmax[nPix][sub_idx];      // By Claude on 06/13/2026
    					for (int vx = 0; vx < vout_dim; vx++) {      // 1D max along vy
    						for (int vy = 0; vy < vout_dim; vy++) {
    							double bmax = Double.NEGATIVE_INFINITY;
    							int y0 = Math.max(0, vy - vsr), y1 = Math.min(vout_dim - 1, vy + vsr);
    							for (int nvy = y0; nvy <= y1; nvy++) if (tmp[nvy * vout_dim + vx] > bmax) bmax = tmp[nvy * vout_dim + vx];
    							bm[vy * vout_dim + vx] = bmax;
    						}
    					}
    				}
    			}
    		}};
    	}
    	ImageDtt.startAndJoin(threads);
    	// Stage B: per pixel - window max K, neighbor weights exp(M-K), then per velocity
    	// D = exp(s0-K) + sum_np w*B and P = exp(s-K)/D
    	ai.set(0);
    	for (int ithread = 0; ithread < threads.length; ithread++) {
    		threads[ithread] = new Thread() { public void run() {
    			int    [] np_idx = new int    [win * win];
    			double [] np_w   = new double [win * win];
    			double [] acc    = new double [num_vout];
    			for (int nPix = ai.getAndIncrement(); nPix < npix; nPix = ai.getAndIncrement()) {
    				int px = nPix % rw, py = nPix / rw;
    				for (int sub_idx = 0; sub_idx < num_sub; sub_idx++) {
    					int n_np = 0;
    					double k_ref = null_score;
    					for (int dpy = -vsr; dpy <= vsr; dpy++) {
    						int ny = py + dpy;
    						if (ny < 0 || ny >= rh) continue;
    						for (int dpx = -vsr; dpx <= vsr; dpx++) {
    							int nx = px + dpx;
    							if (nx < 0 || nx >= rw) continue;
    							int nroi = ny * rw + nx;
    							double mn = pix_max[sub_idx][nroi];
    							if (mn > k_ref) k_ref = mn;
    							np_idx[n_np] = nroi;
    							np_w  [n_np] = mn;
    							n_np++;
    						}
    					}
    					double d0 = Math.exp(null_score - k_ref);
    					Arrays.fill(acc, 0.0);
    					for (int i = 0; i < n_np; i++) {
    						double w = Math.exp(np_w[i] - k_ref);
    						if (w == 0) continue;
    						double [] bn = vel_box[np_idx[i]][sub_idx];
    						for (int v = 0; v < num_vout; v++) acc[v] += w * bn[v];
    					}
    					double [] s   = score [nPix][sub_idx];
    					double [] res = result[nPix][sub_idx];
    					for (int v = 0; v < num_vout; v++) {
    						double den = d0 + acc[v];
    						if (den > 0) {
    							res[v] = Math.min(1.0, Math.exp(s[v] - k_ref) / den);
    						} else {
    							// k_ref (per-pixel-window, all velocities) towers >745 over this cell's
    							// own 4D window near a strong off-velocity score (with noisy background
    							// at small lambda this is a large fraction of cells, not a rare case):
    							// numerator and denominator both underflow to 0 -> 0/0. Redo with the
    							// cell's own window max as reference, reusing precomputed box sums.
    							res[v] = cellPosteriorOwnRef(score, pix_max, vel_box, vel_bmax, // By Claude on 06/13/2026
    									sub_idx, px, py, v, vout_dim, rw, rh, null_score);
    						}
    					}
    				}
    			}
    		}};
    	}
    	ImageDtt.startAndJoin(threads);
    	return result;
    }

    /**
     * Per-cell posterior with the cell's OWN window max as reference, for cells where     // By Claude on 06/13/2026
     * the shared k_ref underflows everything. Reuses the stage-A precomputes:
     *   m_c = max(null, max_np M2[np][v])                  (M2 = velocity-box max)
     *   den = exp(null-m_c) + sum_np exp(M[np]-m_c)*B[np][v]
     * Every product is finite: exp(M[np]-m_c)*B <= 49 mathematically (B is relative to
     * M[np] and the box scores are <= m_c); when M[np]-m_c > 600 the pair (huge exp x
     * underflowed B) is recomputed directly over the 2*vsr+1 squared velocity box.
     * ~50 ops per cell instead of the (2*vsr+1)^4 brute force.
     */
    private double cellPosteriorOwnRef( // By Claude on 06/13/2026
    		final double [][][] score,
    		final double [][]   pix_max,
    		final double [][][] vel_box,
    		final double [][][] vel_bmax,
    		final int sub_idx,
    		final int px,
    		final int py,
    		final int v_out_idx,
    		final int vout_dim,
    		final int rw,
    		final int rh,
    		final double null_score) {
    	final int vsr = vel_suppr_rad;
    	double m_c = null_score;
    	for (int dpy = -vsr; dpy <= vsr; dpy++) {
    		int ny = py + dpy;
    		if (ny < 0 || ny >= rh) continue;
    		for (int dpx = -vsr; dpx <= vsr; dpx++) {
    			int nx = px + dpx;
    			if (nx < 0 || nx >= rw) continue;
    			double m2 = vel_bmax[ny * rw + nx][sub_idx][v_out_idx];
    			if (m2 > m_c) m_c = m2;
    		}
    	}
    	int vy = v_out_idx / vout_dim, vx = v_out_idx % vout_dim;
    	double den = Math.exp(null_score - m_c);
    	for (int dpy = -vsr; dpy <= vsr; dpy++) {
    		int ny = py + dpy;
    		if (ny < 0 || ny >= rh) continue;
    		for (int dpx = -vsr; dpx <= vsr; dpx++) {
    			int nx = px + dpx;
    			if (nx < 0 || nx >= rw) continue;
    			int nroi = ny * rw + nx;
    			double d = pix_max[sub_idx][nroi] - m_c;
    			if (d <= 600.0) {
    				if (d >= -745.0) den += Math.exp(d) * vel_box[nroi][sub_idx][v_out_idx];
    				// else negligible (< 1e-323 relative to the m_c term)
    			} else {
    				// exp(d) overflows while B underflowed - sum this neighbor's box directly
    				double [] sn = score[nroi][sub_idx];
    				for (int dvy = -vsr; dvy <= vsr; dvy++) {
    					int nvy = vy + dvy;
    					if (nvy < 0 || nvy >= vout_dim) continue;
    					int row = nvy * vout_dim;
    					for (int dvx = -vsr; dvx <= vsr; dvx++) {
    						int nvx = vx + dvx;
    						if (nvx < 0 || nvx >= vout_dim) continue;
    						den += Math.exp(sn[row + nvx] - m_c);
    					}
    				}
    			}
    		}
    	}
    	return Math.min(1.0, Math.exp(score[py * rw + px][sub_idx][v_out_idx] - m_c) / den);
    }

    /**
     * Linear 4D "whitening"/unsharp on the conv output: reduces the moment of inertia   // By Claude on 06/13/2026
     * (the spatial + velocity spread = the template-bank Gram autocorrelation, the
     * ~2x-bump blobs seen in -SCORE/-RECT) at the cost of high-frequency noise gain.
     * Separable 3-tap Laplacian unsharp  s' = (1+2a)s - a(s_-1 + s_+1)  applied in vx,vy
     * (within each pixel's velocity grid) and x,y (across ROI pixels), edge-clamped.
     * Pure linear, applied BEFORE the nonlinearities - the freezable sharpening stage.
     * a=0 on an axis = identity there. The exact-deconvolution upgrade (measure the
     * autocorrelation, Wiener-invert) replaces this with the same plug point.
     * Caveats: (1) does NOT recompute the score-normalization energy E through the
     * filter, so RE-TUNE lambda after changing whitening (optimal lambda rises to absorb
     * the noise gain - exactly the budget the whitening recovers). (2) Sub-phase planes
     * are sharpened independently (exact at pixel_decimate=1). Output keeps amplitude
     * units (signed - applied to the raw signed conv output, before any positivity).
     * @param data [roi_npix][nsub][nvel] raw convolve3dDirect output (roi_pix = (y-y0)*rw+(x-x0))
     */
    public double [][][] whitenC5PROI( // By Claude on 06/13/2026
    		final double [][][] data,
    		final Rectangle     roi,
    		final double        a_sp,
    		final double        a_vel) {
    	final int npix     = data.length;
    	final int num_sub  = data[0].length;
    	final int num_vout = data[0][0].length;
    	final int vdim     = (int) Math.round(Math.sqrt(num_vout));
    	final int rw = roi.width, rh = roi.height;
    	final double cv = 1.0 + 2.0 * a_vel, cs = 1.0 + 2.0 * a_sp;
    	final Thread[] threads = ImageDtt.newThreadArray();
    	final AtomicInteger ai = new AtomicInteger(0);
    	// pass 1: velocity unsharp (vx then vy) per (pixel, sub)
    	final double [][][] a1 = new double [npix][num_sub][num_vout];
    	for (int it = 0; it < threads.length; it++) {
    		threads[it] = new Thread() { public void run() {
    			double [] tmp = new double [num_vout];
    			for (int np = ai.getAndIncrement(); np < npix; np = ai.getAndIncrement()) {
    				for (int sub = 0; sub < num_sub; sub++) {
    					double [] s = data[np][sub];
    					for (int vy = 0; vy < vdim; vy++) { int row = vy * vdim;
    						for (int vx = 0; vx < vdim; vx++) {
    							double c = s[row + vx];
    							double l = (vx > 0)        ? s[row + vx - 1] : c;
    							double r = (vx < vdim - 1) ? s[row + vx + 1] : c;
    							tmp[row + vx] = cv * c - a_vel * (l + r);
    						}
    					}
    					double [] o = a1[np][sub];
    					for (int vx = 0; vx < vdim; vx++) {
    						for (int vy = 0; vy < vdim; vy++) { int idx = vy * vdim + vx;
    							double c = tmp[idx];
    							double u = (vy > 0)        ? tmp[idx - vdim] : c;
    							double d = (vy < vdim - 1) ? tmp[idx + vdim] : c;
    							o[idx] = cv * c - a_vel * (u + d);
    						}
    					}
    				}
    			}
    		}};
    	}
    	ImageDtt.startAndJoin(threads);
    	// pass 2: spatial unsharp (x then y) per (sub, velocity) across the rw x rh ROI grid
    	final double [][][] result = new double [npix][num_sub][num_vout];
    	final int n_sv = num_sub * num_vout;
    	ai.set(0);
    	for (int it = 0; it < threads.length; it++) {
    		threads[it] = new Thread() { public void run() {
    			double [] tmpx = new double [rw * rh];
    			for (int sv = ai.getAndIncrement(); sv < n_sv; sv = ai.getAndIncrement()) {
    				int sub = sv / num_vout, v = sv % num_vout;
    				for (int py = 0; py < rh; py++) { int rowp = py * rw;
    					for (int px = 0; px < rw; px++) {
    						double c = a1[rowp + px][sub][v];
    						double l = (px > 0)      ? a1[rowp + px - 1][sub][v] : c;
    						double r = (px < rw - 1) ? a1[rowp + px + 1][sub][v] : c;
    						tmpx[rowp + px] = cs * c - a_sp * (l + r);
    					}
    				}
    				for (int px = 0; px < rw; px++) {
    					for (int py = 0; py < rh; py++) { int idx = py * rw + px;
    						double c = tmpx[idx];
    						double u = (py > 0)      ? tmpx[idx - rw] : c;
    						double d = (py < rh - 1) ? tmpx[idx + rw] : c;
    						result[idx][sub][v] = cs * c - a_sp * (u + d);
    					}
    				}
    			}
    		}};
    	}
    	ImageDtt.startAndJoin(threads);
    	return result;
    }

    /**
     * Posterior-gated amplitude: the lambda-controlled SOFT analog of filterConv5dROI    // By Claude on 06/13/2026
     * (hard 4D NMS) that feeds the recurrent layer. Computes the single-segment Bayes
     * posterior at the operative lambda, then keeps the RAW matched-filter amplitude
     * out[theta] where P[theta] >= keep and zeroes the rest. Amplitude units (not the
     * [0,1] probability) so the recurrent layer's bump fits / energy / rs_min thresholds
     * keep their meaning - lambda only decides WHICH cells survive (the competition),
     * the surviving values are the unmodified linear out. keep=0.5: kept iff this
     * hypothesis is more probable than all its 4D-window rivals AND the null combined.
     * @param raw       [roi_npix][nsub][nvel] - unfiltered convolve3dDirect output
     * @param lambda    operative input-noise sigma (one value, not the sweep list)
     * @param log_prior getVelLogPrior() output, or null for flat
     * @param null_score 0.5*null_z^2
     * @param keep      posterior threshold to retain a cell
     */
    public double [][][] posteriorGateC5PROI( // By Claude on 06/13/2026
    		final double [][][] raw,
    		final Rectangle     roi,
    		final double        lambda,
    		final double []     log_prior,
    		final double        null_score,
    		final double        keep,
    		final boolean       exact) { // By Claude on 06/13/2026
    	double [][][] score = scoreC5PROI(raw, lambda, log_prior, exact);
    	double [][][] post  = softmaxC5PROI(score, roi, null_score);
    	final int num_sub  = raw[0].length;
    	final int num_vout = raw[0][0].length;
    	double [][][] result = new double [raw.length][num_sub][num_vout];
    	for (int np = 0; np < raw.length; np++) {
    		for (int sub = 0; sub < num_sub; sub++) {
    			for (int v = 0; v < num_vout; v++) {
    				result[np][sub][v] = (post[np][sub][v] >= keep) ? raw[np][sub][v] : 0.0;
    			}
    		}
    	}
    	return result;
    }

    /**
     * Lambda-sweep hyperstack of C5P score/posterior arrays for visual tuning            // By Claude on 06/12/2026
     * (resolution-vs-noise balance, like fat zero in 2D phase correlation).
     * Bottom (T) slider - lambda values, Z slider - scenes.
     * marginal=false: RECT layout (velocity square + 1-px gap per pixel, as
     * showConvKernel5d(roi)); marginal=true: 2D sub-resolution pixel map, each pixel =
     * sum over velocities (posterior probability that a target crosses this pixel).
     * @param data          [nlam][nscene][roi_npix][nsub][nvel]
     * @param display_range {min,max} initial display range, or null
     */
    public ImagePlus showC5PLambdaHyper( // By Claude on 06/12/2026
    		double [][][][][] data,
    		Rectangle         roi,
    		String []         scene_titles,
    		String []         lam_titles,
    		String            title,
    		boolean           marginal,
    		double []         display_range) {
    	int num_lam    = data.length;
    	int num_scenes = data[0].length;
    	int num_vout   = kernel5d[0].length;
    	int vout_dim   = (int) Math.sqrt(num_vout);
    	int num_sub    = kernel5d[0][0].length;
    	int sub_dim    = (int) Math.sqrt(num_sub);
    	int img_width  = marginal ? (roi.width  * sub_dim) : (roi.width  * (vout_dim + 1) * sub_dim);
    	int img_height = marginal ? (roi.height * sub_dim) : (roi.height * (vout_dim + 1) * sub_dim);
    	int block_w    = (vout_dim + 1) * sub_dim;
    	if (scene_titles == null) {
    		scene_titles = new String[num_scenes];
    		for (int i = 0; i < num_scenes; i++) scene_titles[i] = String.format("scene-%03d", i);
    	}
    	double [][][] rendered = new double [num_lam][num_scenes][img_width * img_height];
    	if (!marginal) for (double [][] d2 : rendered) for (double [] d1 : d2) Arrays.fill(d1, Double.NaN);
    	for (int nlam = 0; nlam < num_lam; nlam++) {
    		for (int nscene = 0; nscene < num_scenes; nscene++) {
    			for (int py_roi = 0; py_roi < roi.height; py_roi++) {
    				if (roi.y + py_roi < 0 || roi.y + py_roi >= height) continue;
    				for (int px_roi = 0; px_roi < roi.width; px_roi++) {
    					if (roi.x + px_roi < 0 || roi.x + px_roi >= width) continue;
    					double [][] pix_data = data[nlam][nscene][py_roi * roi.width + px_roi]; // [nsub][nvel]
    					for (int sy = 0; sy < sub_dim; sy++) {
    						for (int sx = 0; sx < sub_dim; sx++) {
    							int sub_idx = sy * sub_dim + sx;
    							if (marginal) {
    								double s = 0;
    								for (int v = 0; v < num_vout; v++) s += pix_data[sub_idx][v];
    								int ix = px_roi * sub_dim + sx;
    								int iy = py_roi * sub_dim + sy;
    								rendered[nlam][nscene][iy * img_width + ix] = s;
    							} else {
    								int slot_x = px_roi * block_w + sx * (vout_dim + 1);
    								int slot_y = py_roi * block_w + sy * (vout_dim + 1);
    								for (int vy = 0; vy < vout_dim; vy++) {
    									for (int vx = 0; vx < vout_dim; vx++) {
    										rendered[nlam][nscene][(slot_y + vy) * img_width + slot_x + vx] =
    												pix_data[sub_idx][vy * vout_dim + vx];
    									}
    								}
    							}
    						}
    					}
    				}
    			}
    		}
    	}
    	ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
    			rendered,     // [num_lam][num_scenes][pix]
    			img_width,
    			title,
    			scene_titles, // Z slider (second)
    			lam_titles,   // T slider (bottom)
    			false);
    	if (display_range != null) {
    		imp.setDisplayRange(display_range[0], display_range[1]);
    	}
    	return imp;
    }

}
