Skip to content

I
imagej-elphel
Commit 95e72b22 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

CODEX: Connect CUDA pose Jacobian oracle 

Co-authored-by: 's avatarCodex <codex@elphel.com>
parent e7771dbc
Hide whitespace changes
Inline Side-by-side
Showing with 171 additions and 11 deletions
src/main/java/com/elphel/imagej/cuas/rt/CuasPoseRT.java
View file @ 95e72b22
......@@ -1095,14 +1095,17 @@ public class CuasPoseRT {
}
final ErsCorrection reference_ers = center_CLT.getErsCorrection();
final ErsCorrection scene_ers = scene.getErsCorrection();
final IntersceneLmaFloat.Camera reference_camera = IntersceneLmaFloat.Camera.capture(reference_ers);
final IntersceneLmaFloat.Camera scene_camera = IntersceneLmaFloat.Camera.capture(scene_ers);
final float [] reference_xyz = IntersceneLmaFloat.toFloat3(ZERO3);
final float [] reference_atr = IntersceneLmaFloat.toFloat3(ZERO3);
final float [] scene_xyz = IntersceneLmaFloat.toFloat3(scene_xyzatr0[0]);
final float [] scene_atr = IntersceneLmaFloat.toFloat3(scene_xyzatr0[1]);
final float [] flat_centers = IntersceneLmaFloat.flattenCenters(cm[0]);
final IntersceneLmaFloat.Result float_result = IntersceneLmaFloat.getFxDerivs(
IntersceneLmaFloat.Camera.capture(reference_ers),
IntersceneLmaFloat.Camera.capture(scene_ers),
IntersceneLmaFloat.toFloat3(ZERO3),
IntersceneLmaFloat.toFloat3(ZERO3),
IntersceneLmaFloat.toFloat3(scene_xyzatr0[0]),
IntersceneLmaFloat.toFloat3(scene_xyzatr0[1]),
IntersceneLmaFloat.flattenCenters(cm[0]),
reference_camera, scene_camera,
reference_xyz, reference_atr, scene_xyz, scene_atr,
flat_centers,
float_selection,
null); // validate raw pixel/J first; eigen transform is a separate layer
final IntersceneLmaFloat.Comparison float_comparison = IntersceneLmaFloat.compareToDouble(
......@@ -1110,6 +1113,16 @@ public class CuasPoseRT {
ZERO3, ZERO3, scene_xyzatr0[0], scene_xyzatr0[1], cm[0],
float_selection, null, -1);
System.out.println("CuasPoseRT "+float_comparison);
final IntersceneLmaFloat.Result cuda_result = lmaGpu.execPoseFxJacobian(
reference_camera, scene_camera,
reference_xyz, reference_atr, scene_xyz, scene_atr,
flat_centers, float_selection, true);
if (cuda_result != null) {
final IntersceneLmaFloat.ResultComparison cuda_comparison =
IntersceneLmaFloat.compareResults(cuda_result, float_result);
System.out.println("CuasPoseRT CUDA vs float getFxDerivs: "+
cuda_comparison.format("CUDA", "float"));
}
}
cycle_rms_meas[nlma] = intersceneLma.getLastRms()[1]; // re-measured pure RMS at this cycle's pose // By Claude on 07/12/2026
// pose_lma_debug >= 2: per-inner-step solver lines; header attributes them to
......
src/main/java/com/elphel/imagej/gpu/GpuQuad.java
View file @ 95e72b22
......@@ -52,6 +52,7 @@ import com.elphel.imagej.tileprocessor.Correlation2d;
import com.elphel.imagej.tileprocessor.DttRad2;
import com.elphel.imagej.tileprocessor.GeometryCorrection;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.IntersceneLmaFloat;
import com.elphel.imagej.tileprocessor.QuadCLT;
import Jama.Matrix;
......@@ -2904,6 +2905,24 @@ public class GpuQuad{ // quad camera description
return null;
}
/**
* Primitive-float pose projection/Jacobian gate. The JNA backend launches one
* CUDA thread per tile and returns raw pixel-space fx/J; base/JCuda does not
* carry {@code tp_lma.cu} and returns null.
*/
public IntersceneLmaFloat.Result execPoseFxJacobian(
IntersceneLmaFloat.Camera reference,
IntersceneLmaFloat.Camera scene,
float [] reference_xyz,
float [] reference_atr,
float [] scene_xyz,
float [] scene_atr,
float [] centers,
boolean [] selection,
boolean calculate_jacobian) {
return null;
}
public void execRBGA(
double [] color_weights,
boolean is_lwir,
......
src/main/java/com/elphel/imagej/gpu/jna/GpuQuadJna.java
View file @ 95e72b22
......@@ -6,6 +6,7 @@ import com.elphel.imagej.gpu.TpTask;
import com.elphel.imagej.tileprocessor.Correlation2d;
import com.elphel.imagej.tileprocessor.CorrVector;
import com.elphel.imagej.tileprocessor.GeometryCorrection;
import com.elphel.imagej.tileprocessor.IntersceneLmaFloat;
import com.elphel.imagej.tileprocessor.QuadCLT;
import com.sun.jna.Native;
import com.sun.jna.Pointer;
......@@ -471,6 +472,66 @@ public class GpuQuadJna extends GpuQuad {
return lib.tp_proc_set_peak_debias(proc, debias, (debias == null) ? 0 : debias.length) == 0;
}
// ---- primitive-float pose fx/J (tp_lma.cu, roadmap rung 3-A4b) ---- // By Codex on 07/15/2026
private static float[] poseCameraMeta(final IntersceneLmaFloat.Camera camera) {
return new float[] {
camera.width, camera.height, camera.focalLength, camera.pixelSize,
camera.distortionRadius, camera.disparityRadius, camera.lineTime, camera.stepR};
}
private static void copyPoseVector3(final float[] vector, final float[] packed, final int offset) {
if ((vector == null) || (vector.length < 3)) {
throw new IllegalArgumentException("GpuQuadJna.execPoseFxJacobian: expected 3-vector");
}
System.arraycopy(vector, 0, packed, offset, 3);
}
@Override public IntersceneLmaFloat.Result execPoseFxJacobian(
final IntersceneLmaFloat.Camera reference,
final IntersceneLmaFloat.Camera scene,
final float[] reference_xyz,
final float[] reference_atr,
final float[] scene_xyz,
final float[] scene_atr,
final float[] centers,
final boolean[] selection,
final boolean calculate_jacobian) {
if ((reference == null) || (scene == null) || (centers == null) || ((centers.length % 3) != 0)) {
throw new IllegalArgumentException("GpuQuadJna.execPoseFxJacobian: inconsistent inputs");
}
final int numTiles = centers.length / 3;
if ((numTiles == 0) || ((selection != null) && (selection.length != numTiles))) {
throw new IllegalArgumentException("GpuQuadJna.execPoseFxJacobian: selection length");
}
final float[] poseVectors = new float[12];
copyPoseVector3(reference_xyz, poseVectors, 0);
copyPoseVector3(reference_atr, poseVectors, 3);
copyPoseVector3(scene_xyz, poseVectors, 6);
copyPoseVector3(scene_atr, poseVectors, 9);
final byte[] packedSelection = new byte[numTiles];
for (int tile = 0; tile < numTiles; tile++) {
packedSelection[tile] = (byte) (((selection == null) || selection[tile]) ? 1 : 0);
}
final float[] fx = new float[IntersceneLmaFloat.NUM_COMPONENTS * numTiles];
final float[] nativeJt = calculate_jacobian ?
new float[IntersceneLmaFloat.NUM_PARAMS * IntersceneLmaFloat.NUM_COMPONENTS * numTiles] :
new float[1]; // keep the JNA array argument non-null; native ignores it in fx-only mode
final byte[] valid = new byte[numTiles];
final int rc = lib.tp_proc_exec_pose_fx_jacobian(
proc,
poseCameraMeta(reference), reference.radial,
reference.rByRDist, reference.rByRDist.length, reference.ers, reference.ers.length,
poseCameraMeta(scene), scene.radial,
scene.rByRDist, scene.rByRDist.length, scene.ers, scene.ers.length,
poseVectors, centers, packedSelection, numTiles, calculate_jacobian ? 1 : 0,
fx, nativeJt, valid);
if (rc != 0) {
throw new IllegalStateException("GpuQuadJna.execPoseFxJacobian rc="+rc+": "+lib.tp_last_error());
}
int numValid = 0;
for (byte value : valid) if (value != 0) numValid++;
return new IntersceneLmaFloat.Result(
numTiles, fx, calculate_jacobian ? nativeJt : null, valid, numValid);
}
// ---- LMA normal-equation products (tp_lma.cu, roadmap rung 3) ---- // By Codex on 07/14/2026
@Override public double[] execLmaNormalProducts(
double[] weights, double[][] jt, double[] ymfxWeighted) {
......
src/main/java/com/elphel/imagej/gpu/jna/Stage0.java
View file @ 95e72b22
......@@ -10,7 +10,7 @@ import com.sun.jna.Pointer;
* com.elphel.imagej.gpu.jna.Stage0 [kernel_src_dir] [libcudadevrt.a]
*/
public class Stage0 {
private static final int EXPECTED_KERNELS = 26; // 19 base + 4 consolidate + 2 peak + 1 LMA products
private static final int EXPECTED_KERNELS = 29; // 19 base + 2 conditioning + 4 consolidate + 2 peak + pose fx/J + products
public static void main(String[] args) {
String srcdir = (args.length > 0) ? args[0] : "/home/elphel/git/tile_processor_gpu/src";
String devrt = (args.length > 1) ? args[1] : "/usr/local/cuda/targets/x86_64-linux/lib/libcudadevrt.a";
......
src/main/java/com/elphel/imagej/gpu/jna/TpJna.java
View file @ 95e72b22
......@@ -12,8 +12,9 @@ public interface TpJna extends Library {
/** NVRTC-compile the kernels in srcdir (+ getTpDefines), cuLink(libcudadevrt), load the module.
* Returns an opaque module handle, or null on failure (see tp_last_error()). */
Pointer tp_create_module(String srcdir, String libcudadevrt);
/** Number of kernel functions resolved in the module (28 expected through LWIR conditioning:
* 19 base + 2 conditioning + 4 consolidation + 2 peak + 1 LMA-products), or -1 if handle is null. */
/** Number of kernel functions resolved in the module (29 expected:
* 19 base + 2 conditioning + 4 consolidation + 2 peak + pose-fx/J + LMA-products),
* or -1 if handle is null. */
int tp_module_num_functions(Pointer module);
/** Last error message (NVRTC/CUDA log), empty if none. */
String tp_last_error();
......@@ -136,6 +137,17 @@ public interface TpJna extends Library {
int useFloat);
/** D2H the peak rows (count x 8 floats). Returns the count. */
int tp_proc_get_peaks(Pointer proc, float[] out);
/** Primitive-float raw pose fx/J, before eigen transforms. Camera metadata is
* {width,height,focal,pixel,distortionRadius,disparityRadius,lineTime,stepR}. */
int tp_proc_exec_pose_fx_jacobian(
Pointer proc,
float[] refMeta, float[] refRadial,
float[] refRbr, int refRbrLen, float[] refErs, int refErsLen,
float[] sceneMeta, float[] sceneRadial,
float[] sceneRbr, int sceneRbrLen, float[] sceneErs, int sceneErsLen,
float[] poseVectors, float[] centers, byte[] selection,
int numTiles, int calculateJacobian,
float[] fx, float[] jt, byte[] valid);
/** Deterministic double normal-equation products for the lean pose LMA.
* jt is parameter-major [numParams][numValues]; out is row-major H followed by b.
* Damping and the small solve remain in Java. 0 on success. */
......
src/main/java/com/elphel/imagej/tileprocessor/IntersceneLmaFloat.java
View file @ 95e72b22
......@@ -110,7 +110,7 @@ public final class IntersceneLmaFloat {
public final byte [] valid;
public final int numValid;
private Result(
public Result(
final int numTiles,
final float [] fx,
final float [] jt,
......@@ -124,6 +124,22 @@ public final class IntersceneLmaFloat {
}
}
/** Statistics for comparing two flattened float implementations. */
public static final class ResultComparison {
public int compared;
public int actualOnlyValid;
public int expectedOnlyValid;
public float maxAbsFx;
public float maxAbsJt;
public String format(final String actualLabel, final String expectedLabel) {
return String.format(
"compared=%d, %s-only=%d, %s-only=%d, max|fx|=%g, max|J|=%g",
compared, actualLabel, actualOnlyValid, expectedLabel, expectedOnlyValid,
maxAbsFx, maxAbsJt);
}
}
/** Statistics for the optional float-vs-current-double comparison gate. */
public static final class Comparison {
public int compared;
......@@ -334,6 +350,45 @@ public final class IntersceneLmaFloat {
return comparison;
}
/** Compare raw fx/J and validity from two float implementations. */
public static ResultComparison compareResults(final Result actual, final Result expected) {
if ((actual == null) || (expected == null) || (actual.numTiles != expected.numTiles) ||
(actual.fx.length != expected.fx.length) || (actual.valid.length != expected.valid.length) ||
((actual.jt == null) != (expected.jt == null)) ||
((actual.jt != null) && (actual.jt.length != expected.jt.length))) {
throw new IllegalArgumentException("incompatible float getFxDerivs results");
}
final ResultComparison comparison = new ResultComparison();
for (int tile = 0; tile < actual.numTiles; tile++) {
final boolean actualValid = actual.valid[tile] != 0;
final boolean expectedValid = expected.valid[tile] != 0;
if (actualValid && !expectedValid) {
comparison.actualOnlyValid++;
continue;
}
if (!actualValid && expectedValid) {
comparison.expectedOnlyValid++;
continue;
}
if (!actualValid) continue;
comparison.compared++;
final int vi = NUM_COMPONENTS * tile;
for (int component = 0; component < NUM_COMPONENTS; component++) {
comparison.maxAbsFx = Math.max(comparison.maxAbsFx,
Math.abs(actual.fx[vi + component] - expected.fx[vi + component]));
}
if (actual.jt == null) continue;
for (int par = 0; par < NUM_PARAMS; par++) {
final int ji = par * NUM_COMPONENTS * actual.numTiles + vi;
for (int component = 0; component < NUM_COMPONENTS; component++) {
comparison.maxAbsJt = Math.max(comparison.maxAbsJt,
Math.abs(actual.jt[ji + component] - expected.jt[ji + component]));
}
}
}
return comparison;
}
private static final class WorldDerivs {
final float [] world = new float [3];
final float [] dPixel = new float [9]; // columns: d/d px, py, disparity
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment