CLAUDE: self-documenting comments: MCLT overview, quad algebra conj/no-conj,...

CLAUDE: self-documenting comments: MCLT overview, quad algebra conj/no-conj, tp_task bridge, GPU geometry chain

Comment-only. From Andrey's explanation: MCLT perfect-reconstruction overview
(TileProcessor.h, with blog ref); convolveTiles = quad product WITHOUT
conjugation, kernels DTT-II/III family vs type-IV image tiles;
correlateAccumulateTiles = conj(T1) x T2, order-sensitive; shiftTileHor/Vert
= lossless fractional-shift phase rotators (differential rectification);
struct tp_task bridge role + composed-offset semantics; calc_rot_deriv /
calculate_tiles_offsets = GPU side of historic getPortsCoordinates().
Pre-existing dirty jna/tp_jna.cpp intentionally NOT included.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

src/TileProcessor.cu

@@ -4132,6 +4132,16 @@ __global__ void imclt_rbg(
  * Fractional pixel shift (phase rotation), horizontal. In-place. uses 8 threads (.x)
  * Used in convert_direct() -> convert_correct_tiles() -> convertCorrectTile
  *
+ * This (with shiftTileVert) is the key to "differential rectification": the composed
+ * per-sensor shift (factory kernel center offset + post-calibration misalignment + target
+ * disparity + pose relative to the reference pose, see struct tp_task) is split into a
+ * nearest-integer part (tile selection from the source image) and this +/-0.5 pix residual,
+ * applied in the transform domain as a phase rotation of the (cos,sin) quad pairs.
+ * Fractional shifts are LOSSLESS here and through all TD operations — resampling loss
+ * appears only on the final inverse (imclt) back to pixel domain, so no warping/upsampling
+ * is needed to preserve deep-subpixel (0.05-0.1 pix) disparity resolution.
+ * (Comment by Claude on 07/03/2026, from Andrey's explanation.)
+ *
  * @param clt_tile             transform domain representation of a tile: [4][8][8+1], // +1 to alternate column ports
  * @param residual_shift       fractional pixel shift [-0.5, +0.5)
  */
@@ -4208,6 +4218,17 @@ __device__ void shiftTileVert(
  * Convolve image tile with the kernel tile in transform domain
  * Used in convert_direct() -> convert_correct_tiles() -> convertCorrectTile
  *
+ * Element-wise product in the CLT quad algebra: each tile is 4 quadrants (CC,SC,CS,SS) of
+ * horizontal/vertical DCT-IV/DST-IV — a 2D analog of a complex number, one (cos,sin) pair
+ * per axis — so the product below is (ch1+i*sh1)(cv1+j*sv1) x (ch2+i*sh2)(cv2+j*sv2) with
+ * NO conjugation (compare correlateAccumulateTiles(), which conjugates the first factor).
+ * Implements space-variant aberration-correction deconvolution: the per-tile factory kernel
+ * is the inverted PSF. NOTE the asymmetry: image tiles are type-IV (from MCLT window
+ * folding), while kernels are prepared in the type-II/III DTT family ("DTT3 converted" —
+ * DCT-II/DST-II spectra, see ImageDttCPU.clt_dtt3_kernel() in imagej-elphel) so that this
+ * quad product implements spatial convolution across the half-sample-shifted basis pair.
+ * (Comment by Claude on 07/03/2026, from Andrey's explanation.)
+ *
  * @param clt_tile             transform domain representation of a tile [4][8][8+1], // +1 to alternate column ports
  * @param kernel               transform domain representation of a kernel [4][8][8+1], // +1 to alternate column ports
  */
@@ -4268,6 +4289,15 @@ __device__ void convolveTiles(
  * Calculate 2D correlation of a pair from CLT representation and accumulate with a specified color weight
  * Called from correlate2D()->correlate2D_inner()
  *
+ * In the CLT quad algebra (4 quadrants CC,SC,CS,SS = one (cos,sin) pair per axis) this is
+ * the product conj(T1) x T2, conjugating BOTH axes of the first factor (the sign-negated
+ * terms sit on clt_tile1's sine components) — the transform-domain phase-correlation
+ * numerator, analogous to conj(F1)*F2 in FFT correlation. Correlation is therefore
+ * order-sensitive: swapping the tiles mirrors the correlation peak through zero.
+ * Compare convolveTiles() which multiplies WITHOUT conjugation. Amplitude normalization
+ * (with fat zero) is applied separately (normalizeTileAmplitude / corr_normalize).
+ * (Comment by Claude on 07/03/2026, from Andrey's explanation.)
+ *
  * @param scale                weight of the current component for accumulation.
  * @param clt_tile1            transform domain representation of a tile [4][8][8+1], 4 quadrants of the clt data 1,
  *                             rows extended to optimize shared ports

src/TileProcessor.h

@@ -36,6 +36,33 @@
  * \brief header file for  the Tile Processor for frequency domain
 
  */
+
+/*
+ * OVERVIEW — the transform behind the Tile Processor (by Claude on 07/03/2026, from
+ * Andrey's explanation; details + essential figures:
+ * https://blog.elphel.com/2018/01/complex-lapped-transform-bayer/ ):
+ *
+ * Images are decomposed into 16x16 tiles at stride 8 (50% overlap) with the MCLT — the
+ * lapped-transform family from audio coding (mp3/AAC MDCT) with the PERFECT RECONSTRUCTION
+ * property: half-cosine 16-long windows applied twice (analysis and synthesis) overlap-add
+ * as cos^2+sin^2==1, so processed tiles reassemble seamlessly, with none of the block
+ * artifacts of non-overlapping JPEG-style DCT. The 16-long lapped segments fold into 8-long
+ * transforms, so each 16x16 tile is represented as 4x8x8 quadrants (CC,SC,CS,SS) of
+ * separable DCT-IV/DST-IV — a 2D generalization of a complex value, one (cos,sin) pair per
+ * axis. Like the FFT, the CLT has the convolution/multiplication property, used here for
+ * (a) space-variant aberration-correction deconvolution (convolveTiles; kernels prepared in
+ * the type-II/III DTT family) and (b) inter-channel/inter-scene phase correlation
+ * (correlateAccumulateTiles; conjugated first factor), with lossless fractional-pixel
+ * shifts via TD phase rotation (shiftTileHor/Vert). This underpins deep-subpixel
+ * (0.05-0.1 pix) disparity resolution at extreme distance/baseline ratios: many pixels
+ * (a tile, and in CUAS multiple consecutive frames) participate coherently in each
+ * correlation — the reason correlation-based matching is used rather than sparse features
+ * (SIFT etc.). Developed for RGB 5MPix sensors where lateral chromatic aberration in the
+ * corners ruined the RG/GB order (deconvolution essential); with Boson640 LWIR the larger
+ * pixels make lenses relatively better, so aberration correction matters less but the same
+ * machinery provides differential rectification and correlation (see struct tp_task in
+ * geometry_correction.h).
+ */
 #pragma once
 #ifndef TILE_PROCESSOR_H_
 #define TILE_PROCESSOR_H_

src/geometry_correction.cu

@@ -168,6 +168,15 @@ inline __device__ int get_task_txy_gc(int num_tile, float * gpu_ftasks, int num_
 /**
  * Calculate rotation matrices and derivatives by az, tilt, roll, zoom
  * NUM_CAMS blocks of 3,3,3 tiles
+ *
+ * First step of the GPU-side geometry chain (replacing the historic host-side
+ * GeometryCorrection.getPortsCoordinates()): per-sensor rotations are built from the
+ * relative-pose correction vector, then calculate_tiles_offsets() uses them (and the
+ * derivative matrices, e.g. for lazy-eye fitting) to fill per-tile per-sensor xy offsets
+ * and disp_dist Jacobians into the GPU-resident task array (struct tp_task).
+ * Individual-sensor poses are rotation-only (translations are good enough from the
+ * mechanical construction), combined as an approximately ortho-normalized composition.
+ * (Comment by Claude on 07/03/2026, from Andrey's explanation.)
  */
 extern "C" __global__ void calc_rot_deriv(
 		int                  num_cams,
@@ -325,6 +334,17 @@ extern "C" __global__ void calc_rot_deriv(
 }
 
 
+/**
+ * Fill the per-tile, per-sensor geometry into the GPU-resident task array: for each task
+ * computes xy[cam][2] (the composed differential-rectification shift: average-camera
+ * projection + per-sensor difference + target disparity + relative pose, see struct tp_task
+ * in geometry_correction.h) and disp_dist[cam][4] (d(pixel x,y)/d(disparity,non-disparity)
+ * Jacobian). GPU side of the historic Java getPortsCoordinates(); host now supplies only
+ * txy/target_disparity (+ optional centerXY when uniform_grid==0, e.g. interscene
+ * accumulation) and reads results back via GpuQuad.updateTasks() when needed on the host
+ * (disparity LMA, lazy-eye). Consumed by convert_direct via the integer+fractional shift
+ * split. (Comment by Claude on 07/03/2026, from Andrey's explanation.)
+ */
 extern "C" __global__ void calculate_tiles_offsets(
 		int                  uniform_grid, //==0: use provided centers (as for interscene) , !=0 calculate uniform grid
 		int                  num_cams,

src/geometry_correction.h

@@ -60,6 +60,25 @@
 
 #define SCENE_UNITS_SCALE  0.001 // meters from mm
 #define MIN_DISPARITY      0.01  // minimal disparity to try to convert to world coordinates
+/**
+ * Per-tile task descriptor — the Java<->CUDA bridge (host-side mirror: imagej-elphel
+ * com.elphel.imagej.gpu.TpTask, marshalled as a flat float array, ints bit-cast to float;
+ * layout must match get_task_size() below). An array of these IS the GPU work list:
+ * the host decides which 16x16 (stride 8) tiles to process, kernels consume it.
+ *
+ * xy[NUM_CAMS][2] implements "differential rectification": no image warping — each sensor
+ * is corrected only by its difference from a common ideal average camera (pure radial
+ * distortion, average focal length), expressed per tile as a center-symmetrical deconvolution
+ * kernel plus a pure 2D shift. The full shift = factory kernel center offset + post-calibration
+ * misalignment + target disparity + pose relative to the reference pose; convert_direct splits
+ * it into integer (tile selection) + fractional +/-0.5 pix (transform-domain phase rotators
+ * shiftTileHor/Vert — lossless until return to pixel domain).
+ *
+ * Filled either by the host (historic path, GeometryCorrection.getPortsCoordinates()) or by
+ * the GPU itself from pose + per-tile disparity (calc_rot_deriv() + calculate_tiles_offsets()),
+ * then optionally read back D2H (Java GpuQuad.updateTasks()).
+ * (Documentation by Claude on 07/03/2026, from Andrey's explanation.)
+ */
 struct tp_task {
 	int   task;
 	union {
@@ -71,8 +90,10 @@ struct tp_task {
 	                            // if isnan(centerXY[0]), then txy is used to calculate centerXY and all xy
 	// scale == 0 - old way, just set. Scale !=0 - accumulate. Or make > 0 - set too? only negative - subtract?
 	float scale;                // multiply during direct conversion before accumulating in TD - used for motion blur correction
-	float xy       [NUM_CAMS][2];
-	float disp_dist[NUM_CAMS][4]; // calculated with getPortsCoordinates()
+	float xy       [NUM_CAMS][2]; // per-sensor pixel offset (the composed full shift, see struct comment)
+	float disp_dist[NUM_CAMS][4]; // calculated with getPortsCoordinates(); per-sensor 2x2 Jacobian
+	                              // d(pixel x,y)/d(disparity, non-disparity) — used by the disparity
+	                              // LMA (Corr2dLMA) and lazy-eye/ERS correction. By Claude on 07/03/2026
 };
 
 #define get_task_size(x) (sizeof(struct tp_task)/sizeof(float) - 6 * (NUM_CAMS - x))