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

/**
 * Lean, self-contained per-sensor conditioning for the CUAS real-time / tile-processor path.
 * By Claude on 06/30/2026.
 *
 * Deliberately bypasses the heavy QuadCLT "conditioning over conditioning" instance-array path:
 * it operates directly on raw jp4-domain per-sensor pixel arrays and the (borrowed) calibration
 * values pulled from a QuadCLT (lwir_scales / lwir_offsets / lwir_scales2 and the per-pixel FPN).
 *
 * Pipeline (per sensor), applied in place:
 *   1) Row/Col denoise (optional): estimate a per-row and a per-column average, optionally
 *      high-pass the 1-D profile (preserve the low frequencies for the downstream LoG), subtract.
 *      The HPF is applied to the length-H / length-W average vectors, NOT the full image
 *      (mathematically ~equivalent, far fewer operations).
 *   2) Photometric + fixed-pattern (2-map model, per-pixel-ready):
 *         corrected = scales2*(raw-offset)^2 + scale*(raw-offset)    - FPN(x,y)
 *      matching the existing model (QuadCLT.calibratePhotometric2 forward form) so that with
 *      scales2==0 (the abandoned 2nd-order term) it is bit-identical to today's additive path.
 *      scale = global per-sensor (uniform) for now; the per-pixel *scale* map (multiplicative
 *      dust) is the future improvement, hence the 2-image (scale, offset) framing.
 */
public class CuasConditioning {

	public static class Config {
		public boolean rowcol_enable        = true;   // subtract per-row and per-col noise
		public boolean rowcol_hpf           = true;   // high-pass the 1-D avg profile before subtracting
		public double  rowcol_hpf_sigma     = 8.0;    // px; ~ the LoG scale (KR3-ish) - noise band the LoG would pass
		public double  rowcol_max_abs       = Double.POSITIVE_INFINITY; // clip |pixel| for the average (outlier reject)
		public double  rowcol_weight_outlier= 0.01;   // small weight for out-of-range pixels
		public boolean photometric          = true;   // apply scale*(raw-offset) [+ scales2 quadratic]
		public boolean apply_fpn            = true;    // subtract per-pixel FPN (the offset map's per-pixel part)
		public Config() {}
	}

	/**
	 * Condition all sensors in place.
	 * @param image_data [sensor][0][pixel] mono LWIR, row-major (idx = y*width + x). Modified in place.
	 * @param width      image width in pixels (height = pixels/width).
	 * @param scales     per-sensor global scale  (lwir_scales;  null -> 1.0).
	 * @param offsets    per-sensor global offset (lwir_offsets; null -> 0.0), subtracted before scaling.
	 * @param scales2    per-sensor 2nd-order term (lwir_scales2; null -> 0.0; normally unused).
	 * @param fpn        [sensor][0][pixel] per-pixel additive FPN, subtracted after photometric (null -> none).
	 * @param cfg        configuration (null -> defaults).
	 */
	public static void condition(
			double [][][] image_data,
			int           width,
			double []     scales,
			double []     offsets,
			double []     scales2,
			double [][][] fpn,
			Config        cfg) {
		if (cfg == null) cfg = new Config();
		final int num_sens = image_data.length;
		for (int s = 0; s < num_sens; s++) if (image_data[s] != null && image_data[s][0] != null) {
			final double [] px = image_data[s][0];
			if (cfg.rowcol_enable) {
				subtractRowColSensor(px, width, cfg);
			}
			if (cfg.photometric || cfg.apply_fpn) {
				final double B0 = (scales  != null) ? scales [s] : 1.0;
				final double C0 = (offsets != null) ? offsets[s] : 0.0;
				final double A0 = (scales2 != null) ? scales2[s] : 0.0;
				final double [] fp = (cfg.apply_fpn && fpn != null && fpn[s] != null) ? fpn[s][0] : null;
				for (int i = 0; i < px.length; i++) {
					double p = px[i];
					if (cfg.photometric) {
						final double d = p - C0;
						p = A0 * d * d + B0 * d;   // == B0*(raw-C0) when A0==0 (linear, matches current additive path)
					}
					if (fp != null) p -= fp[i];
					px[i] = p;
				}
			}
		}
	}

	/** Row/Col denoise for one sensor's flat pixel array (in place): subtract per-row, then per-col. */
	public static void subtractRowColSensor(final double [] px, final int width, final Config cfg) {
		final int height = px.length / width;
		// --- per-row profile (average along x for each row y) ---
		double [] rowProf = lineAverage(px, width, height, true, cfg.rowcol_max_abs, cfg.rowcol_weight_outlier);
		if (cfg.rowcol_hpf) rowProf = highpass(rowProf, cfg.rowcol_hpf_sigma);
		for (int y = 0; y < height; y++) {
			final double sub = rowProf[y];
			if (sub != 0.0) {
				final int row0 = y * width;
				for (int x = 0; x < width; x++) px[row0 + x] -= sub;
			}
		}
		// --- per-col profile (average along y for each col x), recomputed on row-corrected data ---
		double [] colProf = lineAverage(px, width, height, false, cfg.rowcol_max_abs, cfg.rowcol_weight_outlier);
		if (cfg.rowcol_hpf) colProf = highpass(colProf, cfg.rowcol_hpf_sigma);
		for (int y = 0; y < height; y++) {
			final int row0 = y * width;
			for (int x = 0; x < width; x++) px[row0 + x] -= colProf[x];
		}
	}

	/**
	 * Outlier-weighted average producing a 1-D profile.
	 * @param perRow true  -> length-height profile, averaging along x (per-row offset);
	 *               false -> length-width  profile, averaging along y (per-col offset).
	 */
	private static double [] lineAverage(
			final double [] px, final int width, final int height,
			final boolean perRow, final double max_abs, final double weight_outlier) {
		final int n = perRow ? height : width;
		final double [] sum = new double [n];
		final double [] wgt = new double [n];
		int idx = 0;
		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {
				final double d = px[idx++];
				if (Double.isNaN(d)) continue;
				final double w = (Math.abs(d) <= max_abs) ? 1.0 : weight_outlier;
				final int k = perRow ? y : x;
				sum[k] += w * d;
				wgt[k] += w;
			}
		}
		for (int k = 0; k < n; k++) sum[k] = (wgt[k] > 0.0) ? (sum[k] / wgt[k]) : 0.0;
		return sum;
	}

	/** 1-D Gaussian high-pass of a profile: x - smooth(x, sigma), reflect-padded. Returns a new array. */
	public static double [] highpass(final double [] x, final double sigma) {
		if (sigma <= 0.0) return x.clone();
		final int r = Math.max(1, (int) Math.ceil(3.0 * sigma));
		final double [] k = new double [2 * r + 1];
		double ksum = 0.0;
		for (int i = -r; i <= r; i++) { final double v = Math.exp(-0.5 * (i * i) / (sigma * sigma)); k[i + r] = v; ksum += v; }
		for (int i = 0; i < k.length; i++) k[i] /= ksum;
		final int n = x.length;
		final double [] hp = new double [n];
		for (int i = 0; i < n; i++) {
			double sm = 0.0;
			for (int j = -r; j <= r; j++) {
				int ii = i + j;
				if (ii < 0)      ii = -ii - 1;          // reflect
				else if (ii >= n) ii = 2 * n - ii - 1;  // reflect
				if (ii < 0) ii = 0; else if (ii >= n) ii = n - 1;
				sm += k[j + r] * x[ii];
			}
			hp[i] = x[i] - sm;
		}
		return hp;
	}

	/** Build the per-sensor per-pixel OFFSET map (post-scale, subtractive) = scale*offset + FPN,
	 *  so that (scale*raw - offset_map) == scale*(raw-offset) - FPN. For the future 2-map API /
	 *  visualization; not required by {@link #condition}. */
	public static double [][] buildOffsetMaps(
			final double [][][] fpn, final double [] scales, final double [] offsets, final int width) {
		final int num_sens = fpn.length;
		final double [][] out = new double [num_sens][];
		for (int s = 0; s < num_sens; s++) if (fpn[s] != null && fpn[s][0] != null) {
			final double B0 = (scales  != null) ? scales [s] : 1.0;
			final double C0 = (offsets != null) ? offsets[s] : 0.0;
			final double [] fp = fpn[s][0];
			final double [] m = new double [fp.length];
			final double base = B0 * C0;
			for (int i = 0; i < fp.length; i++) m[i] = base + fp[i];
			out[s] = m;
		}
		return out;
	}
}
