diff --git a/.ci/Brewfile b/.ci/Brewfile index 8f46e304916c..600a45cfbc76 100644 --- a/.ci/Brewfile +++ b/.ci/Brewfile @@ -48,6 +48,7 @@ brew 'lua@5.4', link: true, overwrite: true brew 'libomp', link: true brew 'ninja' brew 'onnxruntime' +brew 'opencv' brew 'openexr' brew 'openjpeg' brew 'osm-gps-map' diff --git a/.ci/Dockerfile b/.ci/Dockerfile index c27f227b1c1d..d33c31e652ed 100644 --- a/.ci/Dockerfile +++ b/.ci/Dockerfile @@ -89,6 +89,11 @@ RUN rm -rf /var/lib/apt/lists/* && apt-get update && \ liblua5.2-dev \ liblua5.3-dev \ libomp-$LLVM_VER-dev \ + libopencv-calib3d-dev \ + libopencv-core-dev \ + libopencv-features2d-dev \ + libopencv-flann-dev \ + libopencv-imgproc-dev \ libopenexr-dev \ libopenjp2-7-dev \ libosmgpsmap-1.0-dev \ diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e9cf441c1ae7..bed00c6b4f72 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -143,6 +143,11 @@ jobs: liblensfun-dev \ liblua5.4-dev \ libonnxruntime-dev \ + libopencv-calib3d-dev \ + libopencv-core-dev \ + libopencv-features2d-dev \ + libopencv-flann-dev \ + libopencv-imgproc-dev \ libopenexr-dev \ libopenjp2-7-dev \ libosmgpsmap-1.0-dev \ @@ -269,6 +274,7 @@ jobs: libxml2:p lua54:p omp:p + opencv:p openexr:p openjpeg2:p osm-gps-map:p diff --git a/.github/workflows/nightly.yml b/.github/workflows/nightly.yml index 012a33be27d0..841a0386c23e 100644 --- a/.github/workflows/nightly.yml +++ b/.github/workflows/nightly.yml @@ -72,6 +72,11 @@ jobs: liblensfun-bin \ liblensfun-data-v1 \ liblensfun1 \ + libopencv-calib3d-dev \ + libopencv-core-dev \ + libopencv-features2d-dev \ + libopencv-flann-dev \ + libopencv-imgproc-dev \ libosmgpsmap-1.0-dev \ libpango1.0-dev \ libpng-dev \ @@ -258,6 +263,7 @@ jobs: libxml2:p lua54:p omp:p + opencv:p openexr:p openjpeg2:p osm-gps-map:p diff --git a/DefineOptions.cmake b/DefineOptions.cmake index 15fcad411853..e35970c7eba8 100644 --- a/DefineOptions.cmake +++ b/DefineOptions.cmake @@ -35,6 +35,7 @@ option(BUILD_NOISE_TOOLS "Build tools for generating noise profiles" OFF) option(BUILD_CURVE_TOOLS "Build tools for generating base and tone curves" OFF) option(USE_GMIC "Use G'MIC image processing framework." ON) option(USE_ICU "Use ICU - International Components for Unicode." ON) +option(USE_OPENCV "Use OpenCV for HDR exposure-bracket auto-alignment." ON) option(FORCE_COLORED_OUTPUT "Always produce ANSI-colored output (GNU/Clang only)." OFF) option(USE_SDL2 "Enable SDL2 support" ON) diff --git a/README.md b/README.md index ff91b26b1312..d83eb2509a65 100644 --- a/README.md +++ b/README.md @@ -260,6 +260,7 @@ Optional dependencies (minimum version): * libheif 1.13.0 *(for HEIF import & export; also for AVIF import if no libavif)* * libjxl 0.7.0 *(for JPEG XL import & export)* * WebP 0.3.0 *(for WebP import & export)* +* OpenCV 4 — only the core, imgproc, features2d, flann and calib3d modules *(for HDR exposure-bracket auto-alignment)* Optional dependencies (no version requirement): * colord, Xatom *(for fetching the system display color profile)* diff --git a/build.sh b/build.sh index eb3f5fc28cda..30d982ce0d0c 100755 --- a/build.sh +++ b/build.sh @@ -35,7 +35,7 @@ CMAKE_OPTIONS_FROM_CMDLINE="" PRINT_HELP=0 -FEATURES="AI CAMERA COLORD GMIC GRAPHICSMAGICK IMAGEMAGICK JXL KWALLET LIBSECRET LUA MAC_INTEGRATION MAP OPENCL OPENEXR OPENMP UNITY WEBP" +FEATURES="AI CAMERA COLORD GMIC GRAPHICSMAGICK IMAGEMAGICK JXL KWALLET LIBSECRET LUA MAC_INTEGRATION MAP OPENCL OPENCV OPENEXR OPENMP UNITY WEBP" # prepare a lowercase version with a space before and after # it's very important for parse_feature, has no impact in for loop expansions diff --git a/data/darktableconfig.xml.in b/data/darktableconfig.xml.in index 11b99ae117d4..cd2f62e7fc08 100644 --- a/data/darktableconfig.xml.in +++ b/data/darktableconfig.xml.in @@ -14,6 +14,7 @@
+
@@ -384,6 +385,55 @@ do not set the 'uncategorized' entry for tags do not set the 'uncategorized' entry for tags which do not have children + + plugins/lighttable/hdr_merge_auto_align + bool + true + auto-align frames when merging HDR brackets + register and warp exposure brackets onto a reference frame before merging, removing ghosting caused by a shaky tripod or handheld shots.\nonly effective when darktable is built with OpenCV. + + + plugins/lighttable/hdr_merge_auto_reference + bool + false + pick the HDR alignment reference automatically + before merging, examine every bracket and align to the one richest in detail instead of the first frame.\nimproves alignment when the first exposure is a poor template, at the cost of one extra decode pass over the frames.\nrequires auto-align and an OpenCV build. + + + plugins/lighttable/hdr_merge_proxy_scale + float + 0.625 + alignment detail (proxy scale) + fraction of the full sensor resolution at which feature detection and alignment are computed.\nhigher values find more, more distinctive features (better on repetitive or low-contrast scenes) but cost roughly (scale / 0.5)^2 more time and memory.\n0.5 restores the legacy speed; 0.625 matches the reference prototype. + + + plugins/lighttable/hdr_merge_feature_gamma + float + 2.2 + shadow lift for feature detection (gamma) + display-gamma applied to the linear raw before feature detection, to lift shadow detail into the detector's range.\n1.0 disables it. higher values help under-exposed frames but can amplify shadow noise; the percentile stretch already does part of this, so the effect is scene-dependent. + + + plugins/lighttable/hdr_merge_clahe_clip + float + 0.0 + local contrast for feature detection (CLAHE clip) + clip limit of the local contrast enhancement (CLAHE) applied before feature detection. 0 disables it.\nenabling it (e.g. 2.0) helps detect features on extremely under-exposed frames, but on repetitive textures (façades, railings, foliage) it can change descriptor signatures between exposures and cause mis-alignment, so it is off by default. + + + plugins/lighttable/hdr_merge_sift_keypoints + int + 5000 + feature budget per frame (keypoints) + maximum number of SIFT keypoints kept per frame after spatial balancing, before matching.\nbalancing both frames to a common budget keeps a feature-dense frame from flooding the matcher with ambiguous candidates. + + + plugins/lighttable/hdr_merge_debug_images + bool + false + save alignment debug images + write per-frame diagnostic images while merging: the feature-detection input, the detected keypoints and all feature matches (green = inlier, red = outlier).\nsaved as Netpbm (.pgm/.ppm) in a darktable_hdr_align_debug folder in the system temp directory (or the folder named by the DT_HDR_DEBUG_IMAGE_DIR environment variable). for diagnosing alignment problems; off by default. + plugins/lighttable/tagging/case_sensitivity diff --git a/dev-doc/HDR_Alignment_Design.md b/dev-doc/HDR_Alignment_Design.md new file mode 100644 index 000000000000..7191232d131d --- /dev/null +++ b/dev-doc/HDR_Alignment_Design.md @@ -0,0 +1,614 @@ +# HDR Merge Auto-Alignment — Design Document + +## 1. Current HDR merge in darktable + +### 1.1 Entry point and data flow + +``` +lighttable selected-images menu + └─ dt_control_merge_hdr() src/control/jobs/control_jobs.c:2047 + └─ _control_merge_hdr_job_run() control_jobs.c:634 + └─ for each selected image: + dt_imageio_export_with_flags(..., "pre:rawprepare", ...) + └─ _control_merge_hdr_process() control_jobs.c:468 (per-image callback) + └─ normalize accumulator by white level control_jobs.c:678 + └─ dt_imageio_dng_write_float() → "-hdr.dng" + (or "-hdr_NN.dng" if it exists) + └─ dt_image_import() + collection reload +``` + +The export pipeline is run with the stop point `"pre:rawprepare"`, so the buffer +handed to `_control_merge_hdr_process()` via `ivoid` is the **raw CFA mosaic**: + +- single channel (`image.buf_dsc.channels == 1`), +- `float` values, black-subtracted and rescaled to 1.0 saturation, +- still in the Bayer / X-Trans mosaic layout (`image.buf_dsc.filters != 0`). + +### 1.2 What the merge does + +`_control_merge_hdr_process()` is a **streaming accumulator**. The first image +seeds the geometry (`d->first_imgid`, `d->first_filter`, `d->first_xtrans`, +dimensions, orientation, white-balance, color matrix) and allocates two +full-resolution buffers, `d->pixels` and `d->weight`. Every image — including the +first — is then folded into the accumulator with an exposure/saturation-aware +weight (`_envelope()`), per CFA photosite: + +```c +d->pixels[x + wd*y] += w * in * cal; +d->weight[x + wd*y] += w; +``` + +After all frames, the accumulator is normalized and written as a floating-point +DNG mosaic. + +### 1.3 The problem + +The accumulation is **purely positional**: photosite `(x, y)` of every frame is +added to accumulator cell `(x, y)`. If the camera moved between exposures, the +scene content at `(x, y)` differs frame-to-frame, so the merged mosaic is a +weighted average of *misaligned* content → ghosting and loss of resolution. + +There is currently **no registration step**. The only guard is that all frames +must share size and orientation. + +**Goal:** before a non-reference frame is accumulated, estimate the geometric +transform that maps it onto the reference frame and resample its mosaic +accordingly. A shaky tripod / handheld shot needs at least an affine, more +generally a **projective (homography)** transform. + + +## 2. The darktable-specific challenge: aligning a RAW CFA mosaic + +### 2.1 You cannot run SIFT directly on a mosaic + +Adjacent samples in a Bayer/X-Trans mosaic belong to different color channels, so +the raw buffer has a strong high-frequency CFA modulation that is *not* scene +structure. A feature detector would lock onto the Bayer grid, not the image. + +**Solution — luma proxy.** Collapse the CFA to a single **CFA-free grayscale** +luma image, then resample it to a configurable working scale `DT_HDR_PROXY_SCALE` +(default 0.625× of full resolution): + +- The luma at any pixel is the average of a **stride-1 2×2 window**. For *any* 2×2 + patch of a Bayer mosaic, whatever its phase, the four photosites are exactly one + R, one B and two G, so `0.25·(sum) = (R + 2G + B)/4` — full-resolution luma with + no interpolation and no colour bias. That full-res luma is then area-averaged + down to the proxy size. +- X-Trans: a 2×2 patch is not aligned to the 6×6 tile, but the average still + strongly attenuates the CFA modulation; structural detail (what alignment needs) + survives. Documented as an approximation; see §8 open questions. + +The default scale was raised from the legacy 0.5× toward the Python prototype's +0.627× because the linear raw proxy is far less feature-rich than the prototype's +tone-mapped JPEG (e.g. 4242 vs 43626 raw SIFT keypoints on the same frame); a +higher-resolution proxy recovers distinctive detail and discriminates periodic +structure better, at ≈ `(scale/0.5)²` cost in SIFT time. All feature work happens +on this proxy; the estimated homography is in *proxy coordinates* and is rescaled +to full-resolution mosaic coordinates before being applied (§5.4). + +### 2.2 You cannot bilinearly warp a mosaic + +Bilinear interpolation of a mosaic mixes neighboring photosites of **different +colors** → false colors and channel cross-talk. The warped buffer must remain a +**valid CFA mosaic in the reference frame's phase** so the existing accumulator +(which assumes photosite `(x,y)` has color `fcol(y,x)`) stays correct. + +**Solution — CFA-aware (same-color) resampling.** Both frames come from the same +camera/crop, so they share the CFA layout and phase. To fill output photosite +`(x, y)` (whose color is `c = fcol(y, x)`) we: + +1. map it through the full-res homography to a source location `(sx, sy)` in the + moving mosaic, +2. interpolate **only from moving-mosaic photosites whose color is also `c`**, + using a width-2 separable tent weight (≈ bilinear on that color's sublattice). + +Because we always read color `c` from the source and write it to a cell that is +declared color `c`, the mosaic phase is preserved exactly. This is the key piece +that has no analog in the Python script (which warps RGB). + +```mermaid +flowchart LR + A[output photosite x,y
color c = fcol y,x] --> B[H_full · x,y → sx,sy] + B --> C[gather moving photosites with fcol==c
in a small window around sx,sy] + C --> D[width-2 separable tent weights] + D --> E[normalize → out x,y] +``` + +For Bayer the width-2 tent reduces to exact bilinear on the color's period-2 +sublattice. For X-Trans (irregular green lattice) it is a smooth same-color +weighted average — documented as an approximation. + +### 2.3 Streaming reference model + +The merge processes frames one at a time and never holds all frames in memory. +The alignment state therefore caches the reference's **detected SIFT features** +(keypoints + descriptors, plus the CLAHE'd 8-bit proxy for the optional debug +visuals) and CFA metadata — not the raw proxy. `dt_hdr_alignment_set_reference()` +runs SIFT on the reference **once**, and every non-reference frame matches its own +descriptors against that cache (so the reference is detected once per merge, not +N−1 times) and is accumulated immediately. Both float and 8-bit proxy buffers are +freed after the reference features are built. + +**Auto-reference adaptation.** The Python preset picks the richest-feature frame as +reference via a pre-pass (`--auto-reference`). darktable's streaming export does +not naturally allow a cheap pre-pass without exporting twice. Three options were +considered: + +- **(A) First frame as reference** — zero extra cost, matches current + `d->first_imgid` seeding. The default. +- **(B) Middle frame** — the Python default when `--auto-reference` is off. +- **(C) True auto-reference** — a SIFT probe pre-pass over all frames + (`AUTO_REFERENCE_PROBE_DIM = 1500`), then the normal accumulation pass. Most + faithful to the preset, one extra decode pass. + +**Implemented (T12):** option (A) is the default; option (C) is available behind +the opt-in preference `plugins/lighttable/hdr_merge_auto_reference`. When on, a +pre-pass runs the merge export in `probe_mode` (the process callback just calls +`dt_hdr_alignment_probe_features()` and tracks the richest frame), then the +job moves that `imgid` to the front of `params->index` so it becomes both the +geometry seed and the alignment reference. Default **off** because it doubles the +raw-decode cost; the Python HDR preset has it on, so enabling it reproduces that +behavior. + +--- + +## 3. The image alignment pipeline (as built) + +This section is the authoritative, end-to-end description of how a *merge HDR* +action registers and accumulates frames in darktable, as implemented across +`control_jobs.c`, `hdr_alignment.c` and `hdr_alignment_cv.cc`. Sections 1–3 +motivate it; sections 5–6 give the code/structure. + +### 3.1 Data representations + +The pipeline moves a frame through four representations. Knowing which one each +stage consumes is the key to reading the rest of this section. + +| # | Representation | Dims | Built by | Consumed by | +|---|----------------|------|----------|-------------| +| R0 | **CFA mosaic** (single-channel float, pre-demosaic, black-subtracted) | full `wd×ht` | `pre:rawprepare` export | accumulation; warp input/output | +| R1 | **luma proxy** (CFA-free, stride-1 2×2 luma, area-downscaled) | `≈ s·wd × s·ht`, `s` = `DT_HDR_PROXY_SCALE` (0.625) | `_build_proxy` (C, OMP) | transient basis for R2 (not cached) | +| R2 | **8-bit SIFT proxy** (percentile-normalized R1, display-gamma, ×255) | same as R1 | `_proxy_to_u8` (C, OMP) | SIFT feature init | +| H | **homography** (row-major `double[9]`) | — | feature init | warp; gate | + +`H` exists in two coordinate systems: **proxy** coords (what the feature stage +produces) and **full-res** coords (`H_full = S·H_proxy·S⁻¹`, `S=diag(sx,sy,1)` +with `sx,sy ≈ 1/s`, §5.4) used to warp R0. + +### 3.2 Top-level flow + +```mermaid +flowchart TD + A[lighttable: select brackets → merge HDR] --> B[_control_merge_hdr_job_run] + B --> C{align_enabled
and auto_reference?} + C -- yes --> D[probe pre-pass:
export each frame in probe_mode
count SIFT features on R1
→ richest imgid to front of list] + C -- no --> E + D --> E[main pass: for each frame
export pre:rawprepare → _control_merge_hdr_process] + E --> F{first frame?} + F -- yes --> G[seed geometry + set_reference
cache R1/R2 of reference] + F -- no --> H[align_frame: warp R0 onto reference] + G --> I[accumulate in_buf into d->pixels/d->weight] + H --> I + I --> J{more frames?} + J -- yes --> E + J -- no --> K[normalize by white level] + K --> L[dt_imageio_dng_write_float → name-hdr.dng] + L --> M[import + refresh collection] +``` + +`in_buf` is the reference/raw mosaic for the reference frame and the +warped scratch buffer for every aligned frame; on any alignment failure it falls +back to the raw mosaic, so the accumulator math is untouched. + +### 3.3 Per-frame alignment (the core) + +`dt_hdr_alignment_align_frame()` registers one moving frame's R0 onto the cached +reference and returns the warped R0. The SIFT feature-init homography is the +final warp (when it is reliable and sane); there is no separate refinement stage. + +```mermaid +flowchart TD + A[moving CFA mosaic R0] --> B[_build_proxy → R1 luma
stride-1 2×2, area-downscale, OMP] + B --> C[_proxy_to_u8 → R2 8-bit
percentile-normalize + display-gamma, OMP] + C --> D[Stage 1 — feature init
dt_hdr_cv_feature_homography] + D -- H_feature, inliers --> F{Stage 2 — gate:
reliable and sane?} + F -- yes --> H[H_final = H_feature
status OK] + F -- no --> J[H_final = identity
status IDENTITY] + H --> K[Stage 3 — rescale proxy→full §5.4] + K --> L{corner motion < 0.5px
or identity?} + L -- yes --> M[passthrough: copy R0 unchanged] + L -- no --> N[_warp_mosaic_cfa:
CFA-aware same-color warp
C, OMP collapse2] + J --> M + N --> O[warped CFA mosaic R0'] + M --> O +``` + +### 3.4 Stage 1 — feature initialization (`hdr_alignment_cv.cc`, OpenCV) + +Mirrors `estimate_initial_warp_feature_ransac()` for the HDR path. The reference +side of the pipeline (CLAHE → detect → describe → scale floor → spatial balance) +runs **once** in `dt_hdr_cv_features_create()` (called from `set_reference`) and is +cached; only the **moving** frame is detected per call in +`dt_hdr_cv_feature_homography()`. So in the diagram below the "SIFT detect … +balance" chain executes per frame only for the moving proxy — the reference +keypoints/descriptors come straight from the cache: + +```mermaid +flowchart TD + A[moving R2 proxy
display-gamma encoded
+ cached reference features] --> A2[CLAHE local-contrast
optional, off by default] + A2 --> B[SIFT detect
contrastThreshold 0.04] + B --> C[scale floor: drop kp.size < 6px] + C --> C2[spatial balance keypoints
to 5000 each] + C2 --> D[SIFT compute descriptors] + D --> E[FLANN kNN k=2, both directions] + E --> F[Lowe ratio 0.75
+ mutual-NN consistency] + F --> G[spatial subsample to 1800
over 6x6 grid] + G --> H[findHomography RANSAC
reproj 2.5px, conf 0.995] + H --> I{inliers ≥ 50
and ratio ≥ 0.40?} + I -- no --> J[estimateAffine2D fallback
→ promote to 3×3] + I -- yes --> L + J --> L{inliers in ≤ 2 cells
and MAD ≤ 5px?} + L -- yes --> M[refit translation-only
from median displacement] + L -- no --> K[H_feature, inliers] + M --> K +``` + +Output: `H_feature` in proxy coords + inlier count. If SIFT/matches are too few, +returns identity with 0 inliers (the frame is then left unaligned). + +**SIFT spatial keypoint balancing (`spatialBalanceKeypoints`, parity with the +Python "SIFT spatial balance" step).** After the scale floor, both keypoint sets +are capped to `kSiftSpatialBalanceTarget` (5000), keeping the strongest (by SIFT +response) in each cell of a frame-spanning grid of ~5000 cells. This was the +decisive fix for a real bracket where the C port locked onto a translation one +structure-period short of the true (large) camera motion: the reference frame had +**12888** surviving keypoints against the moving frame's **5075**, a 2.5× +asymmetry that floods the mutual-NN matcher with near-duplicate candidates and +lets RANSAC consense on an *aliased* small shift (−18 px, with sub-pixel reproj) +instead of the correct large one (−243 px, which Python recovers). Balancing both +sets to a common, spatially-even 5000 (matching Python's `7016→5000 / 6626→5000`) +removes the asymmetry and the local over-density that feed the alias. + +**How the linear raw is made detector-friendly (and why *not* CLAHE).** The Python +tool runs on tone-mapped 8-bit JPEGs; darktable runs on the *linear* raw mosaic. +The right transform to bridge that gap is a **display gamma** on the proxy +(`DT_HDR_PROXY_FEATURE_GAMMA`, §3.1): it reproduces the JPEG's tonal encoding and +makes the keypoint count comparable to the prototype's. CLAHE — which an earlier +revision applied unconditionally — is **off by default**, matching the Python HDR +preset (§2.1, "CLAHE: False"): the prototype explicitly disables it because on +*repetitive* textures it changes descriptor signatures between exposures and +manufactures false matches. That is precisely the period-aliasing failure mode +seen on a real façade pair, where CLAHE inflated detection to **279 581** raw +keypoints (vs the prototype's 53 789 on the same frame) and the matcher consensed +on a period-aliased small shift. It remains a runtime knob (`hdr_merge_clahe_clip`) +for genuinely feature-starved / extreme-DR brackets. + +**Two further robustness additions beyond the Python HDR path:** + +1. **Spatial subsampling before RANSAC** (`spatialSubsample`, 6×6 grid, budget + `kMaxMatchesForRansac` = 1800). Replaces the old distance cap: distributing + the matches across the frame keeps RANSAC from being dominated by one dense + bright region, which otherwise biases the homography. +2. **Cluster degradation** (`degradeClusteredToTranslation`): if the RANSAC + inliers still fall into ≤ `kClusterDegradeMaxCells` (2) grid cells and their + per-axis displacement MAD ≤ `kClusterTranslationMaxMad` (5 px), an 8-DOF fit + overfits scale/shear and extrapolates wildly across the unconstrained rest of + the frame. The model is then refit as a pure translation from the median + inlier displacement — the correct, conservative model for a near-translation + handheld bracket. + +These are gated so they cannot hurt the clean case: on a well-featured synthetic +pair neither the subsample nor the cluster-degrade fired (inliers already spanned +the grid). + +### 3.5 Stage 2 — reliability gate (`hdr_alignment.c`, C) + +A single gate decides whether to apply the feature-init warp `H_feature`: + +1. **Reliable** — at least `DT_HDR_FEATURE_MIN_INLIERS` (50) RANSAC inliers. +2. **Sane** (`_warp_is_sane`) — translation < 0.30×diagonal and the 2×2 column-norm + scale ∈ [0.5, 2.0]. + +If both hold the warp is applied (`status OK`); otherwise the frame is accumulated +unwarped (`status IDENTITY` — never worse than the legacy, alignment-free merge). + +### 3.6 Stage 3 — CFA-aware warp (`hdr_alignment.c`, C, OMP) + +`H_final` is rescaled to full-res (§5.4) and applied by `_warp_mosaic_cfa`: for +each output photosite `(x,y)` of CFA color `c`, map through `H_full` to a source +point in the moving R0 and interpolate **only from same-color photosites** with a +width-2 separable tent (≈ bilinear on color `c`'s sublattice). This preserves the +reference frame's mosaic phase, so the warped R0' drops straight into the existing +accumulator. A near-identity `H_final` (corner motion < 0.5 px) short-circuits to +a plain copy to avoid needlessly resampling a frame that did not move. + +### 3.7 Where the time goes / parallelism + +| Stage | Cost | Threading | +|-------|------|-----------| +| `_build_proxy`, `_proxy_to_u8`, `_percentile_bounds` | full-res read | `DT_OMP_FOR` (C) | +| SIFT detect/describe + RANSAC | seconds (proxy res) | OpenCV internal | +| FLANN kNN (two directions) | proxy res | `omp parallel sections num_threads(2)` | +| `_warp_mosaic_cfa` | full-res, the hotspot | `DT_OMP_FOR collapse(2)` (C) | + +The SIFT detect/describe and RANSAC steps run on OpenCV's own internal threading. +The two FLANN match directions (image→template and template→image, for the +mutual-NN test) are independent, so darktable runs them concurrently in an `omp +parallel sections` block capped to two threads. Their KD-trees are **not** rebuilt +symmetrically per frame: the reference (template) index is trained **once** in +`dt_hdr_cv_features_create()` and cached on `RefFeatures`, so the image→template +direction reuses it and only the template→image direction builds a fresh index +over the moving descriptors. darktable's `DT_OMP_FOR` covers the per-pixel C hot +paths, the heaviest of which is the full-resolution CFA warp. The +`_percentile_bounds` stretch used by `_proxy_to_u8` is two O(n) parallel reductions +(min/max, then a histogram array-reduction) rather than a serial sort. + +A per-frame that does not move (reliable warp, sub-pixel corner motion) is +accumulated straight from the caller's source buffer — `align_frame` returns +FALSE and skips the full-resolution copy into the scratch buffer. The +auto-reference probe builds its proxy directly at the probe resolution +(`AUTO_REFERENCE_PROBE_DIM`) instead of the full `DT_HDR_PROXY_SCALE` proxy, since +SIFT there runs at the smaller size anyway. + +--- + +## 4. Proposed architecture + +### 4.1 Module layout + +``` +src/common/hdr_alignment.h Public C API + OpenCV backend seam (extern "C") +src/common/hdr_alignment.c Pure-C: proxy build, normalization, warp math, + CFA-aware resample, orchestration, gates +src/common/hdr_alignment_cv.cc OpenCV (C++) backend: SIFT, FLANN, findHomography, + low-res feature probe +``` + +### 4.2 The C / C++ split + +OpenCV 4's API is **C++ only** (the legacy C `cv.h` API was removed). +A `.c` translation unit cannot include `` or call `cv::SIFT`. + +The proposed resolution is a **narrow C-ABI seam**: + +- `hdr_alignment.c` holds everything that is naturally C and darktable-shaped: CFA + proxy extraction, percentile normalization, warp matrix algebra, the CFA-aware + warp, the reliability gate and the per-frame orchestration. +- The primitives that genuinely need OpenCV are declared `extern "C"` in + `hdr_alignment.h` and implemented in `hdr_alignment_cv.cc`. The reference is + detected once and cached behind an opaque handle, so the seam is four functions + rather than a single symmetric detect-both call: + + ```c + // Detect + cache the reference frame's SIFT features (once per merge). + void *dt_hdr_cv_features_create(const uint8_t *proxy, int width, int height, + int balance_target, double clahe_clip); + void dt_hdr_cv_features_destroy(void *features); + + // Match a moving proxy against the cached reference features -> homography. + int dt_hdr_cv_feature_homography(const void *ref_features, const uint8_t *img, + int width, int height, int balance_target, + double clahe_clip, const char *debug_dir, + int frame_index, double H[9], + dt_hdr_cv_feature_stats_t *stats); + + // Auto-reference probe: SIFT keypoint count on a low-res copy. + int dt_hdr_cv_count_features(const uint8_t *luma, int width, int height, int probe_dim); + ``` + + All image data crosses the seam as plain pointers + dimensions; the homography + crosses as a `double[9]` (row-major); the cached reference crosses as an opaque + `void *`. No OpenCV type appears in any `.c` or public header. `debug_dir` / + `frame_index` are owned by the per-merge caller, so no global debug state is + kept and concurrent merges stay independent. + +```mermaid +flowchart TD + subgraph C [hdr_alignment.c — pure C] + P1[CFA → reduced-res luma proxy] --> P2[percentile normalize + gamma → u8] + P2 --> P3[orchestration] + P3 --> P8[reliability gate + sanity] + P3 --> P6[scale H proxy→full-res] + P6 --> P7[CFA-aware same-color warp] + end + subgraph CC [hdr_alignment_cv.cc — OpenCV C++] + B0[(cached reference features
detected once in set_reference)] + B1[SIFT detect+compute
moving proxy only] + B2[FLANN kNN + Lowe + mutual] + B3[findHomography RANSAC + affine fallback
+ cluster→translation degrade] + end + P3 -- "reference u8 proxy (once)" --> B0 + P3 -- "moving u8 proxy" --> B1 + B0 --> B2 + B1 --> B2 --> B3 -- "H, inliers" --> P3 +``` + +### 4.3 Public API (`hdr_alignment.h`) + +```c +typedef struct dt_hdr_align_t dt_hdr_align_t; // opaque alignment state + +// Runtime-tunable parameters, seeded from the compile-time defaults and +// overridable per run from the user's preferences (see §5). +typedef struct dt_hdr_align_params_t +{ + double proxy_scale; // feature-proxy size as a fraction of full res + double feature_gamma; // display-gamma applied to the 8-bit SIFT proxy (1.0 = off) + double clahe_clip; // pre-SIFT CLAHE clip limit (0 = off) + int sift_keypoints; // per-frame SIFT budget after spatial balancing + int debug_images; // write per-frame alignment debug visuals (0 = off) +} dt_hdr_align_params_t; + +void dt_hdr_alignment_default_params(dt_hdr_align_params_t *p); + +// Create/destroy the per-merge state. `params` may be NULL (use defaults); the +// values are clamped to sane ranges and copied. Registration always uses a +// projective (homography) model with an affine fallback on weak support — there +// is no motion-model argument. +dt_hdr_align_t *dt_hdr_alignment_new(const dt_hdr_align_params_t *params); +void dt_hdr_alignment_free(dt_hdr_align_t *a); + +// Cache the reference frame: builds the reduced-res 8-bit luma proxy and runs +// SIFT on it once, storing the detected features (see §2.3). +gboolean dt_hdr_alignment_set_reference(dt_hdr_align_t *a, + const float *mosaic, int width, int height, + uint32_t filters, const uint8_t (*xtrans)[6]); + +// Align one non-reference frame onto the reference. On success writes the warped +// mosaic into `out` (caller-allocated, width*height) and returns TRUE. On failure +// (too few features, drift gate, etc.) returns FALSE and copies input → out so the +// caller can fall back to the unaligned frame. +gboolean dt_hdr_alignment_align_frame(dt_hdr_align_t *a, + const float *mosaic, float *out, + int width, int height, + uint32_t filters, const uint8_t (*xtrans)[6], + dt_hdr_align_result_t *info); + +// Auto-reference pre-pass: SIFT keypoint count of a frame's proxy at the probe +// resolution, used to pick the richest-feature frame as reference. +int dt_hdr_alignment_probe_features(const float *mosaic, int width, int height); +``` + +### 4.4 Proxy → full-resolution homography rescale + +The proxy is `DT_HDR_PROXY_SCALE` of full resolution, so a full-res point `p_full` +maps to proxy point `p_proxy = S⁻¹ p_full` with `S = diag(sx, sy, 1)`, where +`sx = wd/pw`, `sy = ht/ph` are the actual full-res / proxy size ratios (≈ `1/s`; +kept independent so the integer rounding of `pw,ph` cannot skew the warp). The +homography estimated in proxy coordinates (`H_proxy`, reference→moving) is +converted to full-res with the conjugation: + +``` +H_full = S · H_proxy · S^-1 , S = diag(sx, sy, 1) +``` + +Concretely, for a row-major 3×3 (`_h_scale_proxy_to_full`): + +``` +H[0][1] *= sx/sy ; H[1][0] *= sy/sx (cross terms, by aspect) +H[0][2] *= sx ; H[1][2] *= sy (translation) +H[2][0] /= sx ; H[2][1] /= sy (perspective) +``` + +the isotropic case (`sx = sy`) reduces to the same translation×k / perspective÷k +scaling between resolution levels. + +### 4.5 Reliability gate (`_warp_is_sane`) + +Applied in `hdr_alignment.c` before committing the feature-init warp: + +- **Reliable**: at least `DT_HDR_FEATURE_MIN_INLIERS` (50) RANSAC inliers. +- **Sanity** (`_warp_is_sane`): translation < 0.30 × image diagonal, and the + upper-left 2×2 column-norm scale within [0.5, 2.0]. +- If either fails the frame falls back to **identity** (i.e. accumulate unaligned — + never worse than today's behavior). + +--- + +## 5. Integration into the merge job + +The patch to `_control_merge_hdr_process()` / `_control_merge_hdr_job_run()` +(`src/control/jobs/control_jobs.c`) is additive and behavior-preserving: + +1. `dt_control_merge_hdr_t` gains `dt_hdr_align_t *align`, `gboolean + align_enabled`, a `float *aligned_buf` scratch, and the auto-reference probe + fields (`probe_mode`, `probe_best_count`, `probe_best_imgid`). (`first_imgid` + was also retyped `uint32_t → dt_imgid_t`, which it always semantically was.) +2. `_control_merge_hdr_job_run()` first calls `_control_merge_hdr_validate()` on + the selection, giving a friendly early failure on a non-raw / already-merged / + monochrome / mismatched-geometry selection before any decode. It then reads the + opt-in preferences and, only under `#ifdef HAVE_OPENCV`, creates the state, + seeding the runtime parameters from preferences: + ```c + #ifdef HAVE_OPENCV + d.align_enabled = dt_conf_get_bool("plugins/lighttable/hdr_merge_auto_align"); + if(d.align_enabled) + { + dt_hdr_align_params_t p; + dt_hdr_alignment_default_params(&p); + p.proxy_scale = dt_conf_get_float("plugins/lighttable/hdr_merge_proxy_scale"); + p.feature_gamma = dt_conf_get_float("plugins/lighttable/hdr_merge_feature_gamma"); + p.clahe_clip = dt_conf_get_float("plugins/lighttable/hdr_merge_clahe_clip"); + p.sift_keypoints = dt_conf_get_int ("plugins/lighttable/hdr_merge_sift_keypoints"); + p.debug_images = dt_conf_get_bool ("plugins/lighttable/hdr_merge_debug_images"); + d.align = dt_hdr_alignment_new(&p); + } + #endif + ``` +3. In `_control_merge_hdr_process()`, an authoritative per-frame guard re-checks + the *decoded* descriptor (`filters != 0`, single channel, `TYPE_UINT16`) — the + flags-only pre-check can lag the decode (§2 of the review). Then, just before + the accumulation loop, a single `const float *in_buf` selects the source: + - first frame (`imgid == d->first_imgid`): `dt_hdr_alignment_set_reference()` + and allocate `aligned_buf`; `in_buf = ivoid`. + - later frames: `dt_hdr_alignment_align_frame(ivoid → aligned_buf)`; on success + `in_buf = aligned_buf`, else `in_buf = ivoid`. + The three existing `((float *)ivoid)[…]` reads in the loop now read `in_buf[…]`. +4. The output filename is de-duplicated: `-hdr.dng`, then `-hdr_01.dng`, + `-hdr_02.dng`, … so a re-merge never overwrites an earlier result. +5. `_control_merge_hdr_job_run()` cleanup frees `aligned_buf` and `align`. + +```mermaid +flowchart TD + A[frame arrives in _control_merge_hdr_process] --> B{first image?} + B -- yes --> C[seed d->first_*
set_reference proxy
alloc aligned_buf] + C --> H[accumulate in_buf = ivoid] + B -- no --> D{align_enabled
and aligned_buf?} + D -- no --> H + D -- yes --> E[align_frame ivoid → aligned_buf] + E --> F{success?} + F -- yes --> G[accumulate in_buf = aligned_buf] + F -- no --> H +``` + +The patch is intentionally additive: without OpenCV the `#ifdef` leaves +`align_enabled == FALSE`, `in_buf == ivoid`, and the path is byte-for-byte the +current behavior. With OpenCV, it is gated by the preference +`plugins/lighttable/hdr_merge_auto_align` (default **on**, registered in +`data/darktableconfig.xml.in`). Both `d->first_filter` and `d->first_xtrans` are +passed to the aligner — for Bayer `first_filter` is the crop-shifted pattern, for +X-Trans it is `9u` and `first_xtrans` is the crop-adjusted 6×6, which is exactly +what the CFA color lookup needs. + +Per-frame diagnostics are emitted via `dt_print(DT_DEBUG_HDR_MERGE, …)` +(inliers / corner drift; enable with `-d hdr_merge`), and a one-shot +`dt_control_log()` informs the user when auto-align is active. + +### 5.1 Preferences + +All keys live under `plugins/lighttable/` and are registered in +`data/darktableconfig.xml.in` (section *processing → HDR alignment*). The four +numeric knobs seed `dt_hdr_align_params_t` and are clamped in +`dt_hdr_alignment_new()` to the `DT_HDR_*_MIN/MAX` ranges, which are kept in sync +with the `min`/`max` advertised in the XML. + +| Key | Type (range) | Default | Effect | +|-----|--------------|---------|--------| +| `hdr_merge_auto_align` | bool | **true** | master switch for the whole feature | +| `hdr_merge_auto_reference` | bool | false | run the auto-reference probe pre-pass (§5.2) | +| `hdr_merge_proxy_scale` | float [0.25, 1.0] | 0.625 | feature-proxy scale (`DT_HDR_PROXY_SCALE`) | +| `hdr_merge_feature_gamma` | float [1.0, 6.0] | 2.2 | proxy display-gamma (`1.0` = off) | +| `hdr_merge_clahe_clip` | float [0.0, 16.0] | 0.0 | pre-SIFT CLAHE clip (`0` = off) | +| `hdr_merge_sift_keypoints` | int [500, 20000] | 5000 | per-frame SIFT budget after balancing | +| `hdr_merge_debug_images` | bool | false | write per-frame debug visuals (§5.3) | + +### 5.2 Auto-reference pre-pass + +When `hdr_merge_auto_reference` is on and there is more than one frame, an extra +decode sweep runs the export in `probe_mode`: the process callback only calls +`dt_hdr_alignment_probe_features()` and tracks the richest-feature `imgid`, which +is then moved to the front of `params->index` so it becomes both the geometry seed +and the alignment reference (mirrors `align_and_blend.py --auto-reference`). The +sweep advances the progress bar and honors cancellation between frames; it is +opt-in because it doubles the raw-decode cost. + +### 5.3 Debug visuals + +With `hdr_merge_debug_images` on (or the `DT_HDR_DEBUG_IMAGE_DIR` environment +override), each aligned moving frame writes a numbered set of Netpbm images +(feature-detection input, detected keypoints, and colour-coded matches — green = +inlier, red = outlier) to a per-merge directory resolved once in +`dt_hdr_alignment_new()`. The directory (env or the default under the system temp +dir) is created if missing; the reference-side visuals are written once. The path +and frame index are threaded through the seam, so no global debug state is kept +and concurrent merges do not interfere. Entirely diagnostic. diff --git a/packaging/windows/README.md b/packaging/windows/README.md index e4c642b717a7..6c7c870380d7 100644 --- a/packaging/windows/README.md +++ b/packaging/windows/README.md @@ -15,7 +15,7 @@ The steps to build the darktable Windows executable and make the installer are a * Install required and recommended dependencies for darktable: ```bash - pacman -S --needed mingw-w64-ucrt-x86_64-{libxslt,python-jsonschema,curl,drmingw,exiv2,gettext,gmic,graphicsmagick,gtk3,icu,imath,iso-codes,lcms2,lensfun,libavif,libarchive,libgphoto2,libheif,libjpeg-turbo,libjxl,libpng,libraw,librsvg,libsecret,libtiff,libwebp,libxml2,lua54,openexr,openjpeg2,osm-gps-map,portmidi,potrace,pugixml,SDL2,sqlite3,webp-pixbuf-loader,zlib} + pacman -S --needed mingw-w64-ucrt-x86_64-{libxslt,python-jsonschema,curl,drmingw,exiv2,gettext,gmic,graphicsmagick,gtk3,icu,imath,iso-codes,lcms2,lensfun,libavif,libarchive,libgphoto2,libheif,libjpeg-turbo,libjxl,libpng,libraw,librsvg,libsecret,libtiff,libwebp,libxml2,lua54,opencv,openexr,openjpeg2,osm-gps-map,portmidi,potrace,pugixml,SDL2,sqlite3,webp-pixbuf-loader,zlib} ``` * Install the optional tool for building an installer image (currently x64 only): diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 40edc6397a3c..dfc0fd45b8c3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -54,6 +54,8 @@ FILE(GLOB SOURCE_FILES "common/gpx.c" "common/grouping.c" "common/guided_filter.c" + "common/hdr_alignment.c" + "common/hdr_alignment_cv.cc" "common/heal.c" "common/histogram.c" "common/history.c" @@ -668,6 +670,53 @@ if(USE_GMIC) endif(GMIC_FOUND) endif(USE_GMIC) +if(USE_OPENCV) + # Only these OpenCV modules are needed (the "bare minimum"): + # core/imgproc - image ops; features2d - SIFT + FLANN matcher; + # flann - KD-tree index params; + # calib3d - findHomography / estimateAffine2D. + set(DT_OPENCV_MODULES core imgproc features2d flann calib3d) + # First try the standard config (custom builds, or distros shipping the + # unified OpenCVConfig.cmake). QUIET: absence is fine, the HDR auto-alignment + # then compiles to a no-op (hdr_alignment_cv.cc is empty without HAVE_OPENCV). + find_package(OpenCV 4 COMPONENTS ${DT_OPENCV_MODULES} QUIET) + if(NOT OpenCV_FOUND) + # Fallback for distros (e.g. Debian/Ubuntu) that fragment OpenCV into + # per-module libopencv--dev packages but keep OpenCVConfig.cmake only in + # the heavyweight libopencv-dev metapackage. Discover the modules directly so + # users can install just the handful above instead of the whole library. + find_path(OpenCV_INCLUDE_DIRS NAMES opencv2/core.hpp PATH_SUFFIXES opencv4) + set(OpenCV_LIBS "") + set(_dt_opencv_ok TRUE) + foreach(_mod ${DT_OPENCV_MODULES}) + find_library(DT_OPENCV_${_mod}_LIBRARY NAMES opencv_${_mod}) + if(DT_OPENCV_${_mod}_LIBRARY) + list(APPEND OpenCV_LIBS ${DT_OPENCV_${_mod}_LIBRARY}) + else() + set(_dt_opencv_ok FALSE) + endif() + endforeach() + if(OpenCV_INCLUDE_DIRS AND _dt_opencv_ok + AND EXISTS "${OpenCV_INCLUDE_DIRS}/opencv2/core/version.hpp") + file(STRINGS "${OpenCV_INCLUDE_DIRS}/opencv2/core/version.hpp" _dt_cv_major + REGEX "CV_VERSION_MAJOR[ \t]+[0-9]+") + string(REGEX MATCH "[0-9]+" _dt_cv_major "${_dt_cv_major}") + if(_dt_cv_major GREATER_EQUAL 4) + set(OpenCV_FOUND TRUE) + message(STATUS "OpenCV ${_dt_cv_major}.x found via per-module libraries") + endif() + endif() + endif() + if(OpenCV_FOUND) + add_definitions(-DHAVE_OPENCV) + include_directories(SYSTEM ${OpenCV_INCLUDE_DIRS}) + list(APPEND LIBS ${OpenCV_LIBS}) + message(STATUS "OpenCV: HDR exposure-bracket auto-alignment enabled") + else() + message(STATUS "OpenCV not found: HDR exposure-bracket auto-alignment disabled") + endif(OpenCV_FOUND) +endif(USE_OPENCV) + if(USE_ICU) find_package(ICU COMPONENTS i18n data uc) if(ICU_FOUND) diff --git a/src/common/darktable.c b/src/common/darktable.c index 32a69a14c552..da903a72147e 100644 --- a/src/common/darktable.c +++ b/src/common/darktable.c @@ -264,9 +264,9 @@ static int usage(const char *argv0) " Enable debug output to the terminal (or to the log file if on Windows).\n" " Valid channels are:\n\n" " act_on, ai, cache, camctl, camsupport, control, dev, expose,\n" - " imageio, input, ioporder, lighttable, lua, masks, memory,\n" - " nan, opencl, params, perf, pipe, print, pwstorage, signal,\n" - " sql, tiling, picker, undo\n" + " hdr_merge, imageio, input, ioporder, lighttable, lua, masks,\n" + " memory, nan, opencl, params, perf, pipe, print, pwstorage,\n" + " signal, sql, tiling, picker, undo\n" "\n" " It is also possible to specify names that activate all channels\n" " or a certain subset, as well as increase verbosity:\n" @@ -1157,6 +1157,7 @@ int dt_init(int argc, !strcmp(darg, "expose") ? DT_DEBUG_EXPOSE : !strcmp(darg, "picker") ? DT_DEBUG_PICKER : !strcmp(darg, "ai") ? DT_DEBUG_AI : // AI related stuff. + !strcmp(darg, "hdr_merge") ? DT_DEBUG_HDR_MERGE : // HDR bracket merge + auto-alignment 0; if(dadd) darktable.unmuted |= dadd; diff --git a/src/common/darktable.h b/src/common/darktable.h index 1979f437b403..f4f972deaddb 100644 --- a/src/common/darktable.h +++ b/src/common/darktable.h @@ -381,6 +381,7 @@ typedef enum dt_debug_thread_t DT_DEBUG_EXPOSE = 1 << 26, DT_DEBUG_PICKER = 1 << 27, DT_DEBUG_AI = 1 << 28, + DT_DEBUG_HDR_MERGE = 1 << 29, // HDR exposure-bracket merge + auto-alignment DT_DEBUG_ALL = 0xffffffff & ~DT_DEBUG_VERBOSE, DT_DEBUG_COMMON = DT_DEBUG_OPENCL | DT_DEBUG_PARAMS | DT_DEBUG_IMAGEIO | DT_DEBUG_PIPE | DT_DEBUG_LUA | DT_DEBUG_AI, DT_DEBUG_RESTRICT = DT_DEBUG_VERBOSE | DT_DEBUG_PERF, diff --git a/src/common/hdr_alignment.c b/src/common/hdr_alignment.c new file mode 100644 index 000000000000..5caa36b14733 --- /dev/null +++ b/src/common/hdr_alignment.c @@ -0,0 +1,996 @@ +/* + This file is part of darktable, + Copyright (C) 2026 darktable developers. + + darktable 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. + + darktable 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 darktable. If not, see . +*/ + +/* C port of the HDR use-case (`--preset hdr`) of the reference + * align_and_blend.py. This translation unit owns everything that is naturally + * C and darktable-shaped: + * + * - CFA mosaic -> reduced-resolution, CFA-free luma proxy + * - percentile normalization (np.percentile(1,99) + clip equivalent) + * - proxy-coordinate -> full-resolution homography rescale + * - CFA-aware (same-color) resampling of the full-resolution mosaic + * - warp sanity / corner-drift reliability gates + * - per-frame orchestration (set_reference / align_frame) + * + * The OpenCV-dependent primitives (SIFT, FLANN, findHomography) are reached + * through the C-ABI seam declared in hdr_alignment.h and implemented in + * hdr_alignment_cv.cc. + * + * See dev-doc/HDR_Alignment_Design.md for the full design and the mapping back + * to align_and_blend.py. + */ + +#include "common/hdr_alignment.h" + +#include "common/darktable.h" + +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +/* CFA color lookup. Reproduced from develop/imageop_math.h (FC / FCNxtrans) so + * this lean common file does not have to pull in the OpenCL / imageop headers + * that imageop_math.h depends on. */ +static inline int _fc(const size_t row, const size_t col, const uint32_t filters) +{ + return filters >> (((row << 1 & 14) + (col & 1)) << 1) & 3; +} + +static inline int _fcol(const int row, const int col, const uint32_t filters, + const uint8_t (*const xtrans)[6]) +{ + if(filters == 9u) + // +600 (a multiple of the 6x6 X-Trans period) keeps the index non-negative. + return xtrans[(row + 600) % 6][(col + 600) % 6]; + return _fc((size_t)row, (size_t)col, filters); +} + +/* ------------------------------------------------------------------------- * + * Constants ported from the HDR path of align_and_blend.py. + * (Only the subset reached by `--preset hdr` is reproduced here.) + * ------------------------------------------------------------------------- */ + +// percentile_normalize() bounds +#define DT_HDR_PERCENTILE_LOW 1.0 +#define DT_HDR_PERCENTILE_HIGH 99.0 +#define DT_HDR_SMALL_EPS 1e-6 + +// Default per-frame SIFT keypoint budget after spatial balancing (the backend's +// kSiftSpatialBalanceTarget). Kept here so it can seed the runtime parameter. +#define DT_HDR_SIFT_KEYPOINTS 5000 + +// Feature-proxy downscale relative to the full-resolution mosaic. The Python +// prototype does its SIFT work at ~0.627x of the (already demosaiced) image; +// darktable historically used a 2x2-block proxy (0.5x). A higher-resolution +// proxy carries more distinctive detail -> more SIFT keypoints survive and the +// matcher discriminates periodic structure better (fewer aliased lock-ons), at a +// roughly (scale/0.5)^2 cost in SIFT time. Built CFA-free (see _build_proxy). +// This is the robust half of the distinctiveness work (the linear raw proxy +// yields far fewer SIFT keypoints than the prototype's tone-mapped JPEG -- e.g. +// 4242 vs 43626 raw on the same frame). Override at build time +// (-DDT_HDR_PROXY_SCALE=0.5) to restore the legacy behaviour / speed. +#ifndef DT_HDR_PROXY_SCALE +#define DT_HDR_PROXY_SCALE 0.625 +#endif + +// Perceptual (display-gamma) encoding applied to the 8-bit SIFT proxy after the +// percentile stretch, before CLAHE, to make the linear raw proxy more +// display-like (closer to the prototype's tone-mapped JPEG input). Note the +// preceding percentile stretch already performs the global [1,99] stretch, so +// this is mostly a *redistribution* on top of it: in testing it lifted keypoints +// on mid-key scenes but slightly reduced them on noisy deep-shadow frames (where +// the stretch already amplifies shadow noise). Hence it is exposed as a +// build-time knob: 1.0 disables it (-DDT_HDR_PROXY_FEATURE_GAMMA=1.0); raise it +// for a stronger shadow lift. +#ifndef DT_HDR_PROXY_FEATURE_GAMMA +#define DT_HDR_PROXY_FEATURE_GAMMA 2.2 +#endif + +// CLAHE clip limit applied to the 8-bit SIFT proxy before detection. 0 disables +// it. Off by default: the display-gamma proxy already mimics the prototype's +// tone-mapped JPEG, and -- as align_and_blend.py warns (SIFT_USE_CLAHE = False) +// -- CLAHE on repetitive textures changes descriptor signatures and manufactures +// false matches (the period-aliasing failure mode). Raise it (e.g. 2.0) for +// genuinely feature-starved / extreme-DR brackets. +#ifndef DT_HDR_CLAHE_CLIP +#define DT_HDR_CLAHE_CLIP 0.0 +#endif + +// Longest-side resolution of the auto-reference SIFT probe (AUTO_REFERENCE_PROBE_DIM). +#define DT_HDR_AUTO_REFERENCE_PROBE_DIM 1500 + +// Feature-init reliability (FEATURE_HOMOGRAPHY_MIN_INLIERS). +#define DT_HDR_FEATURE_MIN_INLIERS 50 + +// Warp sanity bounds (WARP_SANITY_*). +#define DT_HDR_WARP_MAX_TRANSLATION_DIAG_FRAC 0.30 +#define DT_HDR_WARP_SCALE_MIN 0.5 +#define DT_HDR_WARP_SCALE_MAX 2.0 + +// Below this full-resolution corner motion the warp is treated as a no-op and +// the mosaic is copied through unresampled (avoids needlessly softening a +// frame that did not actually move). +#define DT_HDR_NOOP_MAX_CORNER_PX 0.5 + +// Half-width (px) of the separable same-color resampling tent. For a Bayer +// period-2 sublattice this reduces to exact bilinear; for X-Trans it is a +// smooth same-color weighted average. +#define DT_HDR_CFA_TENT_RADIUS 2 + +// Accepted ranges for the runtime-tunable parameters. These MUST stay in sync +// with the bounds advertised for the matching +// plugins/lighttable/hdr_merge_* keys in data/darktableconfig.xml.in; the +// clamps in dt_hdr_alignment_new() are the safety net behind that UI. +#define DT_HDR_PROXY_SCALE_MIN 0.25 +#define DT_HDR_PROXY_SCALE_MAX 1.0 +#define DT_HDR_FEATURE_GAMMA_MIN 1.0 +#define DT_HDR_FEATURE_GAMMA_MAX 6.0 +#define DT_HDR_CLAHE_CLIP_MIN 0.0 +#define DT_HDR_CLAHE_CLIP_MAX 16.0 +#define DT_HDR_SIFT_KEYPOINTS_MIN 500 +#define DT_HDR_SIFT_KEYPOINTS_MAX 20000 + +/* ------------------------------------------------------------------------- * + * Alignment state. + * ------------------------------------------------------------------------- */ + +struct dt_hdr_align_t +{ + // Runtime-tunable parameters (proxy scale, feature gamma, CLAHE clip, SIFT + // budget), seeded from the compile-time defaults and overridable per run. + dt_hdr_align_params_t params; + + // Per-merge debug-image policy, resolved once in dt_hdr_alignment_new() (env + // override or the preference). NULL => debug images off. `debug_frame` is + // bumped for each aligned moving frame so its visuals share a number. Owned + // here (not a backend global) so concurrent HDR merges cannot interfere. + gchar *debug_dir; + int debug_frame; + + // Cached reference SIFT features (opaque dt_hdr_cv_features handle), computed + // once from the reference proxy in set_reference and matched against by every + // moving frame -- so the reference is detected once per merge, not N-1 times. + // NULL without OpenCV or if the precompute failed. + void *ref_features; + int pw, ph; // proxy dimensions + + // Reference full-resolution geometry / CFA description. + int width, height; + uint32_t filters; + uint8_t xtrans[6][6]; + + gboolean have_reference; +}; + +/* ------------------------------------------------------------------------- * + * Small 3x3 (row-major) homography helpers. + * + * These (and the gate / CFA-warp / logging helpers further down) are only used + * on the alignment path, which needs OpenCV; guard them so a build without + * OpenCV does not trip -Werror=unused-function on the plain static helpers. + * ------------------------------------------------------------------------- */ +#ifdef HAVE_OPENCV + +static void _h_identity(double H[9]) +{ + H[0] = 1.0; H[1] = 0.0; H[2] = 0.0; + H[3] = 0.0; H[4] = 1.0; H[5] = 0.0; + H[6] = 0.0; H[7] = 0.0; H[8] = 1.0; +} + +// Map a point through a row-major homography. +static inline void _h_apply(const double H[9], double x, double y, double *ox, double *oy) +{ + const double w = H[6] * x + H[7] * y + H[8]; + const double iw = (fabs(w) < 1e-12) ? 1e12 : 1.0 / w; + *ox = (H[0] * x + H[1] * y + H[2]) * iw; + *oy = (H[3] * x + H[4] * y + H[5]) * iw; +} + +// Rescale a homography estimated in proxy coordinates to full-resolution +// coordinates: H_full = S * H_proxy * S^-1 with S = diag(sx, sy, 1), where sx/sy +// are the full-res / proxy size ratios (≈ 1 / DT_HDR_PROXY_SCALE). They are +// nearly equal (rounding of the proxy dims) but kept independent so an off-square +// rounding cannot skew the warp. Expanding S*H*S^-1 gives the per-entry factors +// below (translation scales up, perspective scales down, the cross terms by the +// aspect ratio sx/sy). +static void _h_scale_proxy_to_full(double H[9], double sx, double sy) +{ + H[1] *= sx / sy; // b + H[2] *= sx; // translation x + H[3] *= sy / sx; // d + H[5] *= sy; // translation y + H[6] /= sx; // perspective x + H[7] /= sy; // perspective y +} + +#endif // HAVE_OPENCV + +/* ------------------------------------------------------------------------- * + * Percentile bounds (np.percentile(img, (1, 99)) equivalent). + * + * The bounds only gate an 8-bit normalization, so exact order statistics are + * unnecessary: a uniform histogram over [min,max] gives bucket-quantized + * percentiles that are more than accurate enough, in O(n) and fully parallel + * (two reductions, no full sort). This replaces the previous qsort-of-a-copy, + * which was the dominant serial cost on the alignment path. + * ------------------------------------------------------------------------- */ + +#define DT_HDR_PERCENTILE_BINS 4096 + +static void _percentile_bounds(const float *src, size_t n, double low, double high, + float *p_low, float *p_high) +{ + *p_low = 0.0f; + *p_high = 1.0f; + if(n == 0) return; + + // Pass 1: data range (parallel min/max reduction). + float lo = src[0], hi = src[0]; + DT_OMP_FOR(reduction(min : lo) reduction(max : hi)) + for(size_t i = 0; i < n; i++) + { + const float v = src[i]; + if(v < lo) lo = v; + if(v > hi) hi = v; + } + if(!(hi > lo)) // flat (or single-valued) image: nothing to stretch + { + *p_low = lo; + *p_high = hi; + return; + } + + // Pass 2: uniform histogram over [lo,hi] (parallel array reduction). + uint32_t hist[DT_HDR_PERCENTILE_BINS] = { 0 }; + const double sbin = (double)DT_HDR_PERCENTILE_BINS / ((double)hi - (double)lo); + DT_OMP_FOR(reduction(+ : hist[:DT_HDR_PERCENTILE_BINS])) + for(size_t i = 0; i < n; i++) + { + int b = (int)(((double)src[i] - (double)lo) * sbin); + if(b < 0) b = 0; else if(b >= DT_HDR_PERCENTILE_BINS) b = DT_HDR_PERCENTILE_BINS - 1; + hist[b]++; + } + + // Walk the cumulative counts to the requested order statistics (nearest-rank), + // reporting each bin's center as its representative value. + const double bw = ((double)hi - (double)lo) / (double)DT_HDR_PERCENTILE_BINS; + const uint64_t rank_lo = (uint64_t)(low / 100.0 * (double)(n - 1) + 0.5); + const uint64_t rank_hi = (uint64_t)(high / 100.0 * (double)(n - 1) + 0.5); + uint64_t cum = 0; + gboolean got_lo = FALSE; + *p_low = lo; + *p_high = hi; + for(int b = 0; b < DT_HDR_PERCENTILE_BINS; b++) + { + const uint64_t next = cum + hist[b]; + const float center = (float)((double)lo + ((double)b + 0.5) * bw); + if(!got_lo && next > rank_lo) { *p_low = center; got_lo = TRUE; } + if(next > rank_hi) { *p_high = center; break; } + cum = next; + } +} + +/* ------------------------------------------------------------------------- * + * CFA mosaic -> CFA-free luma proxy at a configurable scale. + * + * Luma comes from a stride-1 2x2 window: for *any* 2x2 patch of a Bayer mosaic + * (whatever the phase) the four photosites are exactly one R, one B and two G, + * so 0.25*(sum) is the CFA-free luma (R + 2G + B)/4 at that location -- full + * resolution, no interpolation, no colour bias. (For X-Trans a 2x2 patch is not + * tile-aligned, but averaging still strongly attenuates the CFA modulation, as + * before.) That full-resolution luma is then area-averaged down to the target + * proxy size, so the proxy scale is decoupled from the 0.5x the old 2x2-block + * reduction was locked to. A higher scale keeps more of the structural detail + * SIFT relies on; see DT_HDR_PROXY_SCALE. + * ------------------------------------------------------------------------- */ + +// CFA-free luma at a stride-1 2x2 window with top-left (sx, sy), edge-clamped. +static inline float _luma2x2(const float *m, int width, int height, int sx, int sy) +{ + if(sx < 0) sx = 0; else if(sx > width - 2) sx = width - 2; + if(sy < 0) sy = 0; else if(sy > height - 2) sy = height - 2; + const float a = m[(size_t)sy * width + sx]; + const float b = m[(size_t)sy * width + sx + 1]; + const float c = m[(size_t)(sy + 1) * width + sx]; + const float d = m[(size_t)(sy + 1) * width + sx + 1]; + return 0.25f * (a + b + c + d); +} + +static float *_build_proxy(const float *mosaic, int width, int height, double scale, + int *pw_out, int *ph_out) +{ + if(width < 4 || height < 4) return NULL; + int pw = (int)lround((double)width * scale); + int ph = (int)lround((double)height * scale); + if(pw < 8 || ph < 8) return NULL; + if(pw > width) pw = width; + if(ph > height) ph = height; + + float *proxy = dt_alloc_align_float((size_t)pw * ph); + if(!proxy) return NULL; + + // Source mosaic pixels per proxy pixel (>= ~1; equals 1/scale). + const double fx = (double)width / (double)pw; + const double fy = (double)height / (double)ph; + + DT_OMP_FOR() + for(int ty = 0; ty < ph; ty++) + { + int y0 = (int)(ty * fy); + int y1 = (int)((ty + 1) * fy); + if(y1 <= y0) y1 = y0 + 1; + if(y1 > height - 1) y1 = height - 1; + if(y0 > height - 2) y0 = height - 2; + for(int tx = 0; tx < pw; tx++) + { + int x0 = (int)(tx * fx); + int x1 = (int)((tx + 1) * fx); + if(x1 <= x0) x1 = x0 + 1; + if(x1 > width - 1) x1 = width - 1; + if(x0 > width - 2) x0 = width - 2; + float sum = 0.0f; + int cnt = 0; + for(int sy = y0; sy < y1; sy++) + for(int sx = x0; sx < x1; sx++) + { + sum += _luma2x2(mosaic, width, height, sx, sy); + cnt++; + } + proxy[(size_t)ty * pw + tx] = (cnt > 0) ? sum / (float)cnt : 0.0f; + } + } + + *pw_out = pw; + *ph_out = ph; + return proxy; +} + +// Build the 8-bit SIFT proxy from a float luma proxy: percentile-stretch, apply a +// perceptual (display-gamma) encoding, then scale to [0,255]. The gamma step +// lifts shadow detail out of the linear raw signal into SIFT's operating range +// (see DT_HDR_PROXY_FEATURE_GAMMA); local-contrast enhancement (CLAHE) is applied +// afterwards, in the backend, just before SIFT detection. Mirrors +// _to_feature_uint8() feeding display-referred data. +// +// The stretch and the gamma curve are fused into a single parallel pass (no +// intermediate [0,1] buffer), and the gamma is read from a small LUT so the hot +// loop is a clamp + table lookup instead of a per-pixel powf(). +#define DT_HDR_GAMMA_LUT_SIZE 4096 +static uint8_t *_proxy_to_u8(const float *proxy_f, int pw, int ph, double gamma) +{ + const size_t n = (size_t)pw * ph; + uint8_t *out = dt_alloc_align_uint8(n); + if(!out) return NULL; + + float p_low, p_high; + _percentile_bounds(proxy_f, n, DT_HDR_PERCENTILE_LOW, DT_HDR_PERCENTILE_HIGH, &p_low, &p_high); + const float scale = 1.0f / ((p_high - p_low) + (float)DT_HDR_SMALL_EPS); + + // Precompute the display-gamma response over the normalized [0,1] domain. + const float inv_gamma = 1.0f / (float)(gamma > 0.0 ? gamma : 1.0); + uint8_t lut[DT_HDR_GAMMA_LUT_SIZE]; + for(int k = 0; k < DT_HDR_GAMMA_LUT_SIZE; k++) + { + const float u = (float)k / (float)(DT_HDR_GAMMA_LUT_SIZE - 1); + const float v = powf(u, inv_gamma) * 255.0f; + lut[k] = (uint8_t)(v < 0.0f ? 0.0f : (v > 255.0f ? 255.0f : v)); + } + + DT_OMP_FOR() + for(size_t i = 0; i < n; i++) + { + float nv = (proxy_f[i] - p_low) * scale; // percentile stretch + nv = nv < 0.0f ? 0.0f : (nv > 1.0f ? 1.0f : nv); // clip to [0,1] + const int idx = (int)(nv * (float)(DT_HDR_GAMMA_LUT_SIZE - 1) + 0.5f); + out[i] = lut[idx]; // gamma encode + quantize + } + return out; +} + +/* ------------------------------------------------------------------------- * + * Reliability gates (ported from _warp_sanity_check). + * ------------------------------------------------------------------------- */ +#ifdef HAVE_OPENCV + +static gboolean _warp_is_sane(const double H[9], int w, int h) +{ + const double diag = hypot((double)w, (double)h); + // Translation as a fraction of the image diagonal. + if(hypot(H[2], H[5]) > DT_HDR_WARP_MAX_TRANSLATION_DIAG_FRAC * diag) return FALSE; + // Scale = column norms of the upper-left 2x2. + const double s1 = hypot(H[0], H[3]); + const double s2 = hypot(H[1], H[4]); + if(s1 < DT_HDR_WARP_SCALE_MIN || s1 > DT_HDR_WARP_SCALE_MAX) return FALSE; + if(s2 < DT_HDR_WARP_SCALE_MIN || s2 > DT_HDR_WARP_SCALE_MAX) return FALSE; + return TRUE; +} + +// Per-corner displacement between two warps; fills the 4 distances (px) and +// returns the maximum. Median is derived by the caller. +static double _corner_drift(const double Ha[9], const double Hb[9], int w, int h, double dist[4]) +{ + const double cx[4] = { 0.0, (double)w, 0.0, (double)w }; + const double cy[4] = { 0.0, 0.0, (double)h, (double)h }; + double maxd = 0.0; + for(int i = 0; i < 4; i++) + { + double ax, ay, bx, by; + _h_apply(Ha, cx[i], cy[i], &ax, &ay); + _h_apply(Hb, cx[i], cy[i], &bx, &by); + dist[i] = hypot(ax - bx, ay - by); + if(dist[i] > maxd) maxd = dist[i]; + } + return maxd; +} + +/* ------------------------------------------------------------------------- * + * CFA-aware (same-color) resampling. + * + * Fill output photosite (x, y) -- whose CFA color is _fcol(y, x) -- by mapping + * it through the full-resolution homography into the moving mosaic and + * interpolating ONLY from moving-mosaic photosites of the same color, using a + * separable width-(2*radius) tent. Because we always read color c and write + * it to a cell declared color c, the reference frame's mosaic phase is + * preserved exactly (both frames share the camera's CFA layout). + * ------------------------------------------------------------------------- */ + +static inline float _sample_cfa_same_color(const float *mosaic, int w, int h, + uint32_t filters, const uint8_t (*xtrans)[6], + int c, double sx, double sy) +{ + // Clamp the sampling center to the valid region (clamp-to-edge border). + if(sx < 0.0) sx = 0.0; + if(sy < 0.0) sy = 0.0; + if(sx > (double)(w - 1)) sx = (double)(w - 1); + if(sy > (double)(h - 1)) sy = (double)(h - 1); + + const int x0 = (int)floor(sx) - DT_HDR_CFA_TENT_RADIUS; + const int x1 = (int)floor(sx) + DT_HDR_CFA_TENT_RADIUS; + const int y0 = (int)floor(sy) - DT_HDR_CFA_TENT_RADIUS; + const int y1 = (int)floor(sy) + DT_HDR_CFA_TENT_RADIUS; + const double inv_r = 1.0 / (double)DT_HDR_CFA_TENT_RADIUS; + + double acc = 0.0; + double wsum = 0.0; + for(int j = y0; j <= y1; j++) + { + if(j < 0 || j >= h) continue; + const double wy = 1.0 - fabs(sy - (double)j) * inv_r; + if(wy <= 0.0) continue; + for(int i = x0; i <= x1; i++) + { + if(i < 0 || i >= w) continue; + if(_fcol(j, i, filters, xtrans) != c) continue; + const double wx = 1.0 - fabs(sx - (double)i) * inv_r; + if(wx <= 0.0) continue; + const double wgt = wx * wy; + acc += wgt * (double)mosaic[(size_t)j * w + i]; + wsum += wgt; + } + } + if(wsum > 0.0) return (float)(acc / wsum); + + // Degenerate fallback: no same-color tap had positive tent weight. This only + // happens at extreme X-Trans borders, where the clamped tent window can shrink + // to a patch that carries none of colour c. Search outward in rings for the + // nearest in-bounds photosite that *does* carry colour c, so we never write a + // different-colour sample into a cell declared colour c (which would corrupt + // the mosaic phase). The X-Trans period is 6, so a same-colour site is always + // found within a few rings; the loop is bounded for safety regardless. + int bx = (int)(sx + 0.5); + int by = (int)(sy + 0.5); + if(bx < 0) bx = 0; else if(bx > w - 1) bx = w - 1; + if(by < 0) by = 0; else if(by > h - 1) by = h - 1; + for(int rad = 0; rad <= 8; rad++) + { + float best = 0.0f; + double best_d2 = 0.0; + gboolean found = FALSE; + for(int j = by - rad; j <= by + rad; j++) + { + if(j < 0 || j >= h) continue; + for(int i = bx - rad; i <= bx + rad; i++) + { + if(i < 0 || i >= w) continue; + // ring only: skip interior cells already searched at a smaller radius + if(rad > 0 && abs(i - bx) != rad && abs(j - by) != rad) continue; + if(_fcol(j, i, filters, xtrans) != c) continue; + const double d2 = (sx - i) * (sx - i) + (sy - j) * (sy - j); + if(!found || d2 < best_d2) { best_d2 = d2; best = mosaic[(size_t)j * w + i]; found = TRUE; } + } + } + if(found) return best; + } + // Should be unreachable for a valid CFA; keep a defined result just in case. + return mosaic[(size_t)by * w + bx]; +} + +// Bayer fast path for the same-color sampler. A Bayer CFA is period-2, so in a +// given row the photosites of color c sit on a single column parity (for R/B the +// color occupies only alternate rows; green occupies every row but flips parity). +// We therefore find that parity from the canonical 2x2 cell and step the inner +// loop by 2, visiting exactly the same-color taps with no per-tap _fcol() test. +// The weights and clamp-to-edge fallback are identical to the general sampler, so +// the result is bit-for-bit the same -- only ~1/3 the inner iterations. +static inline float _sample_bayer_same_color(const float *mosaic, int w, int h, + uint32_t filters, int c, double sx, double sy) +{ + if(sx < 0.0) sx = 0.0; + if(sy < 0.0) sy = 0.0; + if(sx > (double)(w - 1)) sx = (double)(w - 1); + if(sy > (double)(h - 1)) sy = (double)(h - 1); + + const int x0 = (int)floor(sx) - DT_HDR_CFA_TENT_RADIUS; + const int x1 = (int)floor(sx) + DT_HDR_CFA_TENT_RADIUS; + const int y0 = (int)floor(sy) - DT_HDR_CFA_TENT_RADIUS; + const int y1 = (int)floor(sy) + DT_HDR_CFA_TENT_RADIUS; + const double inv_r = 1.0 / (double)DT_HDR_CFA_TENT_RADIUS; + + double acc = 0.0; + double wsum = 0.0; + for(int j = y0; j <= y1; j++) + { + if(j < 0 || j >= h) continue; + const double wy = 1.0 - fabs(sy - (double)j) * inv_r; + if(wy <= 0.0) continue; + + // Column parity carrying color c in this row (& 1 is well-defined for the + // negative window edge in two's complement: it yields the 0/1 parity). + const int row_par = j & 1; + int col_par; + if(_fc((size_t)row_par, 0, filters) == c) col_par = 0; + else if(_fc((size_t)row_par, 1, filters) == c) col_par = 1; + else continue; // this row carries no color c (e.g. an R/B row of the wrong parity) + + int istart = x0; + if((istart & 1) != col_par) istart++; + for(int i = istart; i <= x1; i += 2) + { + if(i < 0 || i >= w) continue; + const double wx = 1.0 - fabs(sx - (double)i) * inv_r; + if(wx <= 0.0) continue; + const double wgt = wx * wy; + acc += wgt * (double)mosaic[(size_t)j * w + i]; + wsum += wgt; + } + } + if(wsum > 0.0) return (float)(acc / wsum); + + const int ix = (int)(sx + 0.5); + const int iy = (int)(sy + 0.5); + return mosaic[(size_t)iy * w + ix]; +} + +static void _warp_mosaic_cfa(const float *mosaic, float *out, int width, int height, + uint32_t filters, const uint8_t (*xtrans)[6], const double H[9]) +{ + // Hoist the CFA dispatch out of the per-pixel loop: Bayer uses the period-2 + // fast path, X-Trans the general same-color sampler. + if(filters != 9u) + { + DT_OMP_FOR(collapse(2)) + for(int y = 0; y < height; y++) + for(int x = 0; x < width; x++) + { + const int c = _fcol(y, x, filters, xtrans); + double sx, sy; + _h_apply(H, (double)x, (double)y, &sx, &sy); + out[(size_t)y * width + x] = + _sample_bayer_same_color(mosaic, width, height, filters, c, sx, sy); + } + } + else + { + DT_OMP_FOR(collapse(2)) + for(int y = 0; y < height; y++) + for(int x = 0; x < width; x++) + { + const int c = _fcol(y, x, filters, xtrans); + double sx, sy; + _h_apply(H, (double)x, (double)y, &sx, &sy); + out[(size_t)y * width + x] = + _sample_cfa_same_color(mosaic, width, height, filters, xtrans, c, sx, sy); + } + } +} + +#endif // HAVE_OPENCV + +/* ------------------------------------------------------------------------- * + * Public API. + * ------------------------------------------------------------------------- */ + +#ifdef HAVE_OPENCV +// Log a 3x3 homography (row-major) and its decomposition, matching the Python +// reference's log_warp() / decompose_warp() (MOTION_HOMOGRAPHY) so the C output +// can be compared line-by-line with align_and_blend.py. +static void _log_warp(const char *title, const double H[9]) +{ + dt_print(DT_DEBUG_HDR_MERGE, " %s:", title); + dt_print(DT_DEBUG_HDR_MERGE, " [ %.6f %.6f %.6f ]", H[0], H[1], H[2]); + dt_print(DT_DEBUG_HDR_MERGE, " [ %.6f %.6f %.6f ]", H[3], H[4], H[5]); + dt_print(DT_DEBUG_HDR_MERGE, " [ %.6f %.6f %.6f ]", H[6], H[7], H[8]); + + // Decompose the (homography-normalized) linear part via QR: L = Q*R, with Q a + // pure rotation and R upper-triangular [[sx, shear], [0, sy]]. + double h22 = H[8]; + if(fabs(h22) < DT_HDR_SMALL_EPS) h22 = (h22 >= 0.0 ? DT_HDR_SMALL_EPS : -DT_HDR_SMALL_EPS); + const double tx = H[2] / h22, ty = H[5] / h22; + const double px = H[6] / h22, py = H[7] / h22; + const double a = H[0] / h22, b = H[1] / h22, c = H[3] / h22, d = H[4] / h22; + const double r11 = hypot(a, c); + double rot_deg = 0.0, scale_x = r11, scale_y = 0.0, shear = 0.0; + if(r11 > DT_HDR_SMALL_EPS) + { + const double r12 = (a * b + c * d) / r11; + const double c2x = b - r12 * (a / r11); + const double c2y = d - r12 * (c / r11); + rot_deg = atan2(c, a) * (180.0 / M_PI); // atan2(Q[1,0], Q[0,0]) + scale_x = r11; + scale_y = hypot(c2x, c2y); + shear = r12; + } + dt_print(DT_DEBUG_HDR_MERGE, " Decomposition:"); + dt_print(DT_DEBUG_HDR_MERGE, " translation_x: %.2f", tx); + dt_print(DT_DEBUG_HDR_MERGE, " translation_y: %.2f", ty); + dt_print(DT_DEBUG_HDR_MERGE, " perspective_x: %.2f", px); + dt_print(DT_DEBUG_HDR_MERGE, " perspective_y: %.2f", py); + dt_print(DT_DEBUG_HDR_MERGE, " rotation_deg: %.2f deg", rot_deg); + dt_print(DT_DEBUG_HDR_MERGE, " scale_x: %.2f", scale_x); + dt_print(DT_DEBUG_HDR_MERGE, " scale_y: %.2f", scale_y); + dt_print(DT_DEBUG_HDR_MERGE, " shear: %.2f", shear); +} + +// Log the SIFT+RANSAC initialization metrics block (mirrors log_feature_init_stats). +static void _log_feature_stats(const dt_hdr_cv_feature_stats_t *s) +{ + dt_print(DT_DEBUG_HDR_MERGE, " SIFT+RANSAC initialization metrics:"); + dt_print(DT_DEBUG_HDR_MERGE, " keypoints (template/image): %d/%d", s->kp_template, s->kp_image); + dt_print(DT_DEBUG_HDR_MERGE, " good matches: %d", s->good_matches); + if(s->good_matches > 0) + dt_print(DT_DEBUG_HDR_MERGE, " mutual-consistent matches: %d", s->good_matches); + if(s->good_matches > 0) + dt_print(DT_DEBUG_HDR_MERGE, " RANSAC inliers: %d/%d (%.0f%%)", + s->inliers, s->good_matches, 100.0 * s->inliers / s->good_matches); + else + dt_print(DT_DEBUG_HDR_MERGE, " RANSAC inliers: %d", s->inliers); + if(s->reproj_mean >= 0.0) + dt_print(DT_DEBUG_HDR_MERGE, " reproj mean/median/max: %.2f / %.2f / %.2f px", + s->reproj_mean, s->reproj_median, s->reproj_max); + else + dt_print(DT_DEBUG_HDR_MERGE, " reproj mean/median/max: n/a"); + dt_print(DT_DEBUG_HDR_MERGE, " transform model: %s", + s->used_translation ? "translation (cluster-degraded)" + : (s->used_affine ? "affine-fallback" : "homography")); +} +#endif // HAVE_OPENCV + +void dt_hdr_alignment_default_params(dt_hdr_align_params_t *p) +{ + if(!p) return; + p->proxy_scale = DT_HDR_PROXY_SCALE; + p->feature_gamma = DT_HDR_PROXY_FEATURE_GAMMA; + p->clahe_clip = DT_HDR_CLAHE_CLIP; + p->sift_keypoints = DT_HDR_SIFT_KEYPOINTS; + p->debug_images = 0; +} + +// Clamp user-supplied parameters to the ranges advertised by the matching +// preferences (see the DT_HDR_*_MIN/MAX defines, kept in sync with +// data/darktableconfig.xml.in). +static double _clampd(double v, double lo, double hi) +{ + return v < lo ? lo : (v > hi ? hi : v); +} +static int _clampi(int v, int lo, int hi) +{ + return v < lo ? lo : (v > hi ? hi : v); +} + +// Resolve the per-merge debug-image directory from the env override or the +// preference. Returns a newly-allocated path (free with g_free) or NULL when +// debug images are off / the directory cannot be created. +static gchar *_resolve_debug_dir(int debug_images) +{ + const char *env = getenv("DT_HDR_DEBUG_IMAGE_DIR"); + gchar *dir = (env && env[0]) + ? g_strdup(env) + : (debug_images + ? g_build_filename(g_get_tmp_dir(), "darktable_hdr_align_debug", NULL) + : NULL); + if(!dir) return NULL; + + // Ensure the directory exists (the env override may name a path that does not + // yet exist); otherwise the backend would silently fail to write every dump. + if(g_mkdir_with_parents(dir, 0755) != 0) + { + dt_print(DT_DEBUG_HDR_MERGE, + " [hdr merge] could not create debug image directory '%s'; disabling debug images", + dir); + g_free(dir); + return NULL; + } + dt_print(DT_DEBUG_HDR_MERGE, " debug images -> %s", dir); + return dir; +} + +dt_hdr_align_t *dt_hdr_alignment_new(const dt_hdr_align_params_t *params) +{ + dt_hdr_align_t *a = calloc(1, sizeof(dt_hdr_align_t)); + if(!a) return NULL; + dt_hdr_alignment_default_params(&a->params); + if(params) + { + a->params.proxy_scale = _clampd(params->proxy_scale, + DT_HDR_PROXY_SCALE_MIN, DT_HDR_PROXY_SCALE_MAX); + a->params.feature_gamma = _clampd(params->feature_gamma, + DT_HDR_FEATURE_GAMMA_MIN, DT_HDR_FEATURE_GAMMA_MAX); + a->params.clahe_clip = _clampd(params->clahe_clip, + DT_HDR_CLAHE_CLIP_MIN, DT_HDR_CLAHE_CLIP_MAX); + a->params.sift_keypoints = _clampi(params->sift_keypoints, + DT_HDR_SIFT_KEYPOINTS_MIN, DT_HDR_SIFT_KEYPOINTS_MAX); + a->params.debug_images = params->debug_images ? 1 : 0; + } + a->debug_dir = _resolve_debug_dir(a->params.debug_images); + return a; +} + +void dt_hdr_alignment_free(dt_hdr_align_t *a) +{ + if(!a) return; +#ifdef HAVE_OPENCV + if(a->ref_features) dt_hdr_cv_features_destroy(a->ref_features); +#endif + g_free(a->debug_dir); + free(a); +} + +gboolean dt_hdr_alignment_set_reference(dt_hdr_align_t *a, + const float *mosaic, + int width, + int height, + uint32_t filters, + const uint8_t (*xtrans)[6]) +{ + if(!a || !mosaic) return FALSE; + + int pw = 0, ph = 0; + float *proxy_f = _build_proxy(mosaic, width, height, a->params.proxy_scale, &pw, &ph); + if(!proxy_f) return FALSE; + + // Build the 8-bit SIFT proxy; the float luma is no longer needed afterwards. + uint8_t *proxy_u8 = _proxy_to_u8(proxy_f, pw, ph, a->params.feature_gamma); + dt_free_align(proxy_f); + if(!proxy_u8) return FALSE; + + a->pw = pw; + a->ph = ph; + a->width = width; + a->height = height; + a->filters = filters; + if(xtrans) + memcpy(a->xtrans, xtrans, sizeof(a->xtrans)); + else + memset(a->xtrans, 0, sizeof(a->xtrans)); + +#ifdef HAVE_OPENCV + // Precompute (and cache) the reference's SIFT features once; every moving frame + // then matches against this instead of re-detecting the reference. On failure + // ref_features stays NULL and frames fall through unaligned (graceful). + if(a->ref_features) dt_hdr_cv_features_destroy(a->ref_features); + a->ref_features = dt_hdr_cv_features_create(proxy_u8, pw, ph, + a->params.sift_keypoints, a->params.clahe_clip); +#endif + dt_free_align(proxy_u8); + + a->have_reference = TRUE; + dt_print(DT_DEBUG_HDR_MERGE, + " reference proxy: %dx%d mosaic -> %dx%d luma (%s)", + width, height, pw, ph, filters == 9u ? "X-Trans" : "Bayer"); + return TRUE; +} + +int dt_hdr_alignment_probe_features(const float *mosaic, int width, int height) +{ +#ifndef HAVE_OPENCV + (void)mosaic; + (void)width; + (void)height; + return 0; +#else + if(!mosaic) return 0; + // The probe ranks frames before any state exists, so it uses the default gamma + // (the ranking only needs to be self-consistent). SIFT runs at the probe + // resolution (longest side <= AUTO_REFERENCE_PROBE_DIM), so build the proxy + // directly at that scale instead of building the full DT_HDR_PROXY_SCALE proxy + // and letting the backend downscale it -- for a large raw that avoids building + // (and normalizing) a proxy several times bigger than the probe ever uses. + const int max_side = (width > height) ? width : height; + double probe_scale = (double)DT_HDR_AUTO_REFERENCE_PROBE_DIM / (double)(max_side > 0 ? max_side : 1); + if(probe_scale > DT_HDR_PROXY_SCALE) probe_scale = DT_HDR_PROXY_SCALE; + int pw = 0, ph = 0; + float *proxy = _build_proxy(mosaic, width, height, probe_scale, &pw, &ph); + if(!proxy) return 0; + uint8_t *u8 = _proxy_to_u8(proxy, pw, ph, DT_HDR_PROXY_FEATURE_GAMMA); + dt_free_align(proxy); + if(!u8) return 0; + // Pass the probe dim through as a safety cap; the proxy is already at or below + // it, so this normally does not resize again. + const int n = dt_hdr_cv_count_features(u8, pw, ph, DT_HDR_AUTO_REFERENCE_PROBE_DIM); + dt_free_align(u8); + return n; +#endif +} + +gboolean dt_hdr_alignment_align_frame(dt_hdr_align_t *a, + const float *mosaic, + float *out, + int width, + int height, + uint32_t filters, + const uint8_t (*xtrans)[6], + dt_hdr_align_result_t *info) +{ + if(info) + *info = (dt_hdr_align_result_t){ .status = DT_HDR_ALIGN_IDENTITY, + .feature_inliers = 0, + .corner_drift = 0.0 }; + + // Default behavior on any early-out: pass the frame through unchanged so the + // caller accumulates the unaligned mosaic (current darktable behavior). + if(!a || !a->have_reference || !mosaic || !out || out == mosaic) + { + if(out && mosaic && out != mosaic) + memcpy(out, mosaic, (size_t)width * height * sizeof(float)); + return FALSE; + } + + // Frames must share geometry with the reference (the merge already enforces + // identical size / orientation upstream). + if(width != a->width || height != a->height) + { + memcpy(out, mosaic, (size_t)width * height * sizeof(float)); + return FALSE; + } + +#ifndef HAVE_OPENCV + // Built without OpenCV: registration is unavailable; accumulate unaligned. + // (filters / xtrans are only consumed by the CFA-aware warp below.) + (void)filters; + (void)xtrans; + if(info) info->status = DT_HDR_ALIGN_DISABLED; + memcpy(out, mosaic, (size_t)width * height * sizeof(float)); + return FALSE; +#else + const int pw = a->pw; + const int ph = a->ph; + + // No reference feature cache (precompute failed / too few features): there is + // nothing to align against, so pass the frame through unaligned. + if(!a->ref_features) + { + memcpy(out, mosaic, (size_t)width * height * sizeof(float)); + return FALSE; + } + + // Build the moving frame's proxies. + int mpw = 0, mph = 0; + float *mov_f = _build_proxy(mosaic, width, height, a->params.proxy_scale, &mpw, &mph); + if(!mov_f || mpw != pw || mph != ph) + { + if(mov_f) dt_free_align(mov_f); + memcpy(out, mosaic, (size_t)width * height * sizeof(float)); + return FALSE; + } + uint8_t *mov_u8 = _proxy_to_u8(mov_f, pw, ph, a->params.feature_gamma); + dt_free_align(mov_f); // float proxy is only an intermediate for the u8 proxy + if(!mov_u8) + { + memcpy(out, mosaic, (size_t)width * height * sizeof(float)); + return FALSE; + } + + // The per-merge debug-image directory (env override or preference) was + // resolved once in dt_hdr_alignment_new(); number each aligned moving frame so + // its visuals share an index. Both are passed to the backend -- no global + // debug state -- so concurrent merges stay independent. + const int frame_index = ++a->debug_frame; + + // --- Stage 1: feature initialization (SIFT + RANSAC homography) ---------- + double H_feature[9]; + _h_identity(H_feature); + dt_hdr_cv_feature_stats_t fstats = { 0 }; + const int feature_inliers = dt_hdr_cv_feature_homography( + a->ref_features, mov_u8, pw, ph, a->params.sift_keypoints, a->params.clahe_clip, + a->debug_dir, frame_index, H_feature, &fstats); + const gboolean feature_ok = feature_inliers >= DT_HDR_FEATURE_MIN_INLIERS; + // (the per-stage SIFT / match / inlier lines come from the backend; the + // structured metrics block is emitted below.) + + dt_free_align(mov_u8); + + // --- Stage 2: choose and apply the warp ---------------------------------- + // The SIFT feature-init homography is the final warp when it is reliable and + // geometrically sane; otherwise the frame is accumulated unaligned (never + // worse than the legacy, alignment-free merge). + double H_final[9]; + dt_hdr_align_status_t status; + if(feature_ok && _warp_is_sane(H_feature, pw, ph)) + { + memcpy(H_final, H_feature, sizeof(H_final)); + status = DT_HDR_ALIGN_OK; + } + else + { + _h_identity(H_final); + status = DT_HDR_ALIGN_IDENTITY; + } + + // Rescale proxy -> full-resolution coordinates. + _h_scale_proxy_to_full(H_final, (double)width / pw, (double)height / ph); + + // Report the full-resolution corner motion vs. identity. + double ident[9]; + _h_identity(ident); + double dist_id[4]; + const double corner_motion = _corner_drift(ident, H_final, width, height, dist_id); + + if(info) + { + info->status = status; + info->feature_inliers = feature_inliers; + info->corner_drift = corner_motion; + } + + // --- Structured log ------------------------------------------------------ + // Warps are reported in full-resolution coordinates so they line up with the + // Python reference, which works at full image resolution. + _log_feature_stats(&fstats); + if(status != DT_HDR_ALIGN_IDENTITY) + _log_warp("Final warp matrix (SIFT feature-init)", H_final); + else + dt_print(DT_DEBUG_HDR_MERGE, + " feature init unreliable (%d inliers): frame left unaligned", + feature_inliers); + + if(status == DT_HDR_ALIGN_IDENTITY || corner_motion < DT_HDR_NOOP_MAX_CORNER_PX) + { + // No reliable warp, or motion below the resample threshold: the caller's + // original `mosaic` is already the correct data to accumulate (an identity / + // sub-pixel warp we deliberately do not resample -- see DT_HDR_NOOP_MAX_...). + // Return FALSE so the caller uses its own source buffer and skip the + // full-frame copy into `out` (a static frame is otherwise a pure ~w*h*4-byte + // memcpy per frame). `info->status` still reports the OK/IDENTITY decision. + return FALSE; + } + + _warp_mosaic_cfa(mosaic, out, width, height, filters, xtrans, H_final); + return TRUE; +#endif // HAVE_OPENCV +} + +// clang-format off +// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py +// vim: shiftwidth=2 expandtab tabstop=2 cindent +// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified; +// clang-format on diff --git a/src/common/hdr_alignment.h b/src/common/hdr_alignment.h new file mode 100644 index 000000000000..5f412149c623 --- /dev/null +++ b/src/common/hdr_alignment.h @@ -0,0 +1,210 @@ +/* + This file is part of darktable, + Copyright (C) 2026 darktable developers. + + darktable 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. + + darktable 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 darktable. If not, see . +*/ + +/* Automatic registration of RAW CFA frames for the lighttable "merge HDR" + * feature. This is a C port of the HDR use-case (`--preset hdr`) of the + * reference `align_and_blend.py`: SIFT feature initialization -> RANSAC + * homography, applied to exposure brackets shot on a shaky tripod or handheld. + * + * Unlike the Python tool, which operates on demosaiced RGB, this code operates + * on the undemosaiced Bayer / X-Trans mosaic that the merge job accumulates. + * It therefore (1) builds a reduced-resolution, CFA-free luma proxy for all + * feature work and (2) applies the resulting homography to the full resolution + * mosaic with a CFA-aware (same-color) resampler that preserves the mosaic + * phase. See dev-doc/HDR_Alignment_Design.md for the full design. + * + * The OpenCV-dependent primitives (SIFT, FLANN, findHomography) live behind a + * narrow C-ABI seam implemented in hdr_alignment_cv.cc, because OpenCV 4's API + * is C++ only. When darktable is + * built without OpenCV (HAVE_OPENCV undefined) the public functions degrade to + * no-ops that report "no alignment", and the merge behaves exactly as before. + */ + +#pragma once + +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +/* Outcome of aligning a single frame, for logging / diagnostics. */ +typedef enum dt_hdr_align_status_t +{ + DT_HDR_ALIGN_OK = 0, // a reliable warp was applied + DT_HDR_ALIGN_IDENTITY, // nothing reliable found, frame left unwarped + DT_HDR_ALIGN_DISABLED // built without OpenCV / alignment off +} dt_hdr_align_status_t; + +typedef struct dt_hdr_align_result_t +{ + dt_hdr_align_status_t status; + int feature_inliers; // RANSAC inlier count from the feature stage + double corner_drift; // max corner displacement vs identity, in pixels +} dt_hdr_align_result_t; + +/* Opaque per-merge alignment state. Holds the cached reference proxy. */ +typedef struct dt_hdr_align_t dt_hdr_align_t; + +/* Runtime-tunable alignment parameters. These mirror the compile-time defaults + * (DT_HDR_PROXY_SCALE etc.) but let the merge job override them per run from the + * user's preferences. Pass NULL to dt_hdr_alignment_new() to use the defaults. */ +typedef struct dt_hdr_align_params_t +{ + double proxy_scale; // feature proxy size as a fraction of full res (≈0.5–1.0) + double feature_gamma; // display-gamma applied to the 8-bit SIFT proxy (1.0 = off) + double clahe_clip; // CLAHE clip limit before SIFT (0 = off; on aids extreme DR + // but can cause false matches on repetitive textures) + int sift_keypoints; // per-frame SIFT keypoint budget after spatial balancing + int debug_images; // write per-frame alignment debug visuals (0 = off) +} dt_hdr_align_params_t; + +/* Fill `p` with the built-in defaults (the DT_HDR_PROXY_* compile-time values). + * Callers can then override individual fields from preferences. */ +void dt_hdr_alignment_default_params(dt_hdr_align_params_t *p); + +/* Create / destroy the alignment state for one HDR merge run. `params` may be + * NULL (use defaults); the values are clamped to sane ranges and copied. + * Registration uses a projective (homography) motion model with an affine + * fallback on weak support -- see hdr_alignment_cv.cc. */ +dt_hdr_align_t *dt_hdr_alignment_new(const dt_hdr_align_params_t *params); +void dt_hdr_alignment_free(dt_hdr_align_t *a); + +/* Cache the reference frame. Builds the reduced-resolution 8-bit SIFT luma + * proxy from the CFA mosaic and stores it on `a`. `filters` and `xtrans` + * describe the CFA (filters == 9u => X-Trans, else Bayer). + * Returns TRUE on success. */ +gboolean dt_hdr_alignment_set_reference(dt_hdr_align_t *a, + const float *mosaic, + int width, + int height, + uint32_t filters, + const uint8_t (*xtrans)[6]); + +/* Align one non-reference frame onto the cached reference. + * + * - On success: writes the warped mosaic into `out` (caller-allocated, + * width*height floats) and returns TRUE. + * - Otherwise returns FALSE and the caller must accumulate the original + * `mosaic` (never worse than the current, alignment-free behavior). `out` is + * NOT guaranteed to be populated when FALSE is returned: for a reliable but + * sub-pixel (static) frame it is deliberately left untouched to skip a + * redundant full-frame copy, and on some early rejections it holds a copy of + * `mosaic`. Callers must therefore read `mosaic`, not `out`, whenever FALSE + * is returned. `info->status` (if non-NULL) still carries the alignment + * decision (OK / IDENTITY / DISABLED) regardless of the return value. + * + * `out` may not alias `mosaic`. `info` may be NULL. */ +gboolean dt_hdr_alignment_align_frame(dt_hdr_align_t *a, + const float *mosaic, + float *out, + int width, + int height, + uint32_t filters, + const uint8_t (*xtrans)[6], + dt_hdr_align_result_t *info); + +/* Probe a frame's feature richness for auto-reference selection: builds the luma + * proxy and returns its SIFT keypoint count at the auto-reference probe + * resolution. Returns 0 if built without OpenCV or the proxy cannot be built. + * Mirrors _count_lowres_sift_features() that drives `--auto-reference`. */ +int dt_hdr_alignment_probe_features(const float *mosaic, int width, int height); + +/* ------------------------------------------------------------------------- * + * OpenCV backend seam (implemented in hdr_alignment_cv.cc). + * + * These are the only entry points that touch OpenCV. All image data crosses + * as plain pointers; the homography crosses as a row-major double[9]. Kept in + * this header (rather than an internal one) so both the C and C++ translation + * units agree on the ABI. Not part of the public darktable API. + * ------------------------------------------------------------------------- */ + +typedef struct dt_hdr_cv_feature_stats_t +{ + int kp_template; // keypoints kept in the reference proxy (after scale floor) + int kp_image; // keypoints kept in the moving proxy (after scale floor) + int kp_template_raw; // keypoints before the scale floor (reference) + int kp_image_raw; // keypoints before the scale floor (moving) + int ratio_matches; // matches passing the Lowe ratio test + int good_matches; // mutual-consistent matches fed to RANSAC + int inliers; // RANSAC inliers supporting the returned homography + int used_affine; // 1 if the affine fallback produced the result + int used_translation; // 1 if the result was refit to a pure translation + // (cluster-degradation); takes precedence over used_affine + // inlier reprojection error (pixels) in proxy coords; < 0 if unavailable + double reproj_mean; + double reproj_median; + double reproj_max; +} dt_hdr_cv_feature_stats_t; + +/* Precompute and cache the reference frame's SIFT features (CLAHE -> detect -> + * describe -> scale floor -> spatial balance) from its 8-bit luma proxy, so the + * reference is detected ONCE per merge instead of being re-detected for every + * moving frame. Returns an opaque handle (free with dt_hdr_cv_features_destroy) + * or NULL on failure. `balance_target` / `clahe_clip` are as for + * dt_hdr_cv_feature_homography(). */ +void *dt_hdr_cv_features_create(const uint8_t *proxy, + int width, + int height, + int balance_target, + double clahe_clip); +void dt_hdr_cv_features_destroy(void *features); + +/* Estimate the reference->image homography between the cached reference features + * (`ref_features`, from dt_hdr_cv_features_create) and a moving 8-bit luma proxy + * `img`, using SIFT + FLANN (Lowe ratio + mutual NN) + RANSAC findHomography, + * with an estimateAffine2D fallback on weak homography support. Only the moving + * frame is detected here; the reference comes from the cache. Mirrors + * estimate_initial_warp_feature_ransac() for the HDR path. + * + * Writes the 3x3 row-major homography into H and returns the inlier count + * (0 => estimation failed or ref_features is NULL; H is left as identity). + * `balance_target` caps the moving keypoint set before matching (<= 0 disables + * balancing); `clahe_clip` is the pre-SIFT CLAHE clip limit (<= 0 disables it) + * and should match the value the reference cache was built with. + * + * `debug_dir` (NULL/"" disables) is the directory the per-frame debug visuals + * are written to, and `frame_index` numbers them; both are owned by the caller + * (the per-merge alignment state), so no global debug state is kept here and + * concurrent merges do not interfere. */ +int dt_hdr_cv_feature_homography(const void *ref_features, + const uint8_t *img, + int width, + int height, + int balance_target, + double clahe_clip, + const char *debug_dir, + int frame_index, + double H[9], + dt_hdr_cv_feature_stats_t *stats); + +/* Count SIFT keypoints on a low-resolution copy of `luma` (longest side scaled + * to <= probe_dim). Used by the optional auto-reference pre-pass to pick the + * richest-feature frame. Mirrors _count_lowres_sift_features(). */ +int dt_hdr_cv_count_features(const uint8_t *luma, int width, int height, int probe_dim); + +#ifdef __cplusplus +} +#endif + +// clang-format off +// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py +// vim: shiftwidth=2 expandtab tabstop=2 cindent +// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified; +// clang-format on diff --git a/src/common/hdr_alignment_cv.cc b/src/common/hdr_alignment_cv.cc new file mode 100644 index 000000000000..d9579c632e41 --- /dev/null +++ b/src/common/hdr_alignment_cv.cc @@ -0,0 +1,771 @@ +/* + This file is part of darktable, + Copyright (C) 2026 darktable developers. + + darktable 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. + + darktable 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 darktable. If not, see . +*/ + +/* OpenCV backend for HDR frame registration (the C-ABI seam declared in + * hdr_alignment.h). This is the only translation unit that touches OpenCV. + * It concentrates the feature-based primitives that the HDR path of + * align_and_blend.py relies on: + * + * dt_hdr_cv_feature_homography <- estimate_initial_warp_feature_ransac() + * dt_hdr_cv_count_features <- _count_lowres_sift_features() + * + * All image data crosses the seam as plain pointers operating on the reduced- + * resolution luma proxy built by hdr_alignment.c; the homography crosses as a + * row-major double[9]. No OpenCV type appears in any public header. + * + * Built only when darktable is configured with OpenCV (HAVE_OPENCV). + */ + +#include "common/hdr_alignment.h" + +#ifdef HAVE_OPENCV + +#include "common/darktable.h" // dt_print / DT_DEBUG_HDR_MERGE + +#include +#include +#include +#include +#include // KDTreeIndexParams / SearchParams + +#include +#include +#include // fopen/fwrite (Netpbm debug dump) +#include // getenv +#include +#include + +namespace +{ + +// --- Constants ported from the HDR path of align_and_blend.py -------------- +constexpr double kSiftContrastThreshold = 0.04; +constexpr double kSiftMinScalePx = 6.0; +constexpr double kRatioThreshold = 0.75; +constexpr int kMinFeatureKeypoints = 25; +constexpr int kMinGoodMatches = 25; +constexpr int kMaxMatchesForRansac = 1800; + +constexpr double kRansacReprojThreshold = 2.5; +constexpr int kRansacMaxIters = 5000; +constexpr double kRansacConfidence = 0.995; +constexpr int kHomographyMinInliers = 50; +constexpr double kHomographyMinInlierRatio = 0.40; +constexpr int kSpatialGrid = 6; // FEATURE_SPATIAL_GRID_ROWS/COLS +constexpr int kClusterDegradeMaxCells = 2; // INLIER_CLUSTER_DEGRADE_MAX_CELLS +constexpr double kClusterTranslationMaxMad = 5.0; // INLIER_CLUSTER_TRANSLATION_MAX_MAD_PX +constexpr int kSiftSpatialBalanceTarget = 5000; // SIFT_SPATIAL_BALANCE_TARGET + +// Wrap a borrowed 8-bit buffer as a single-channel cv::Mat header (no copy). +cv::Mat wrapU8(const uint8_t *p, int w, int h) +{ + return cv::Mat(h, w, CV_8U, const_cast(p)); +} + +// Row-major CV_32F 3x3 -> double[9]. +void fromMat3x3(const cv::Mat &m, double H[9]) +{ + cv::Mat md; + m.convertTo(md, CV_64F); + for(int i = 0; i < 9; i++) H[i] = md.at(i / 3, i % 3); +} + +cv::Ptr makeSift() +{ + return cv::SIFT::create(/*nfeatures=*/0, /*nOctaveLayers=*/3, + kSiftContrastThreshold, /*edgeThreshold=*/10.0, + /*sigma=*/1.6); +} + +// Optional local-contrast enhancement (CLAHE) on the 8-bit feature image, with +// clip limit `clip` (<= 0 disables it). CLAHE recovers features in extreme +// dynamic range, but -- as align_and_blend.py warns and defaults off +// (SIFT_USE_CLAHE = False) -- on *repetitive* textures it changes descriptor +// signatures between frames and manufactures false matches, so on a periodic +// scene it makes the matcher consense on a period-aliased shift. We therefore +// leave it off by default and rely on the display-gamma proxy (which already +// mimics the prototype's tone-mapped JPEG input); it stays available as a knob +// for genuinely feature-starved / extreme-DR brackets. +void applyClahe(cv::Mat &img, double clip) +{ + if(clip <= 0.0) return; + cv::Ptr clahe = cv::createCLAHE(clip, cv::Size(8, 8)); + clahe->apply(img, img); +} + +// Indices of keypoints at or above the scale floor (SIFT_MIN_SCALE_PX): the +// smallest octaves are descriptor-ambiguous noise on low-contrast scenes. +// Returns an index list (rather than mutating in place) so the caller can apply +// the same selection to the parallel descriptor matrix via gatherKpDesc(). +std::vector scaleFloorKeep(const std::vector &kps) +{ + std::vector keep; + keep.reserve(kps.size()); + for(int idx = 0; idx < (int)kps.size(); idx++) + if(kSiftMinScalePx <= 0.0 || kps[idx].size >= kSiftMinScalePx) keep.push_back(idx); + return keep; +} + +// Reorder keypoints and (when present) their descriptor rows by `keep`, in +// keep-order, so the two stay aligned after any filter/balance step. des may be +// empty (no descriptors computed yet), in which case only kps is reordered. +void gatherKpDesc(std::vector &kps, cv::Mat &des, const std::vector &keep) +{ + std::vector nk; + nk.reserve(keep.size()); + const bool has_desc = !des.empty() && des.rows == (int)kps.size(); + cv::Mat nd; + if(has_desc) nd.create((int)keep.size(), des.cols, des.type()); + for(size_t k = 0; k < keep.size(); k++) + { + nk.push_back(kps[keep[k]]); + if(has_desc) des.row(keep[k]).copyTo(nd.row((int)k)); + } + kps.swap(nk); + if(has_desc) des = nd; +} + +// Spatially balance keypoints down to `target`, keeping the strongest (by SIFT +// response) in each cell of a frame-spanning grid. Mirrors align_and_blend.py's +// "SIFT spatial balance" step. Two effects matter for the HDR raw path: +// (1) it caps a feature-dense reference (we have seen 12888 vs 5075 between the +// two frames) so the descriptor matcher is not flooded with near-duplicate +// candidates that pass the ratio test at the *wrong* instance of a periodic +// structure -- the failure mode where RANSAC locks onto a translation that +// is one structure-period short of the true (large) camera motion; +// (2) it makes the template/image keypoint budgets symmetric, which is what +// mutual-nearest-neighbour matching assumes. +// The grid has ~target cells (≈ 1 kept keypoint per cell), proportioned to the +// image aspect ratio, and we round-robin the per-cell response-sorted lists so +// coverage stays uniform if some cells are sparse. Returns the kept indices +// (into kps) for gatherKpDesc() to apply to keypoints and descriptors together. +std::vector spatialBalanceKeep(const std::vector &kps, int width, int height, int target) +{ + std::vector all(kps.size()); + for(int idx = 0; idx < (int)kps.size(); idx++) all[idx] = idx; + if(target <= 0 || (int)kps.size() <= target) return all; + const double aspect = (double)std::max(1, width) / (double)std::max(1, height); + int cols = std::max(1, (int)std::lround(std::sqrt((double)target * aspect))); + int rows = std::max(1, (int)std::lround((double)target / cols)); + std::vector> cells(rows * cols); + for(int idx = 0; idx < (int)kps.size(); idx++) + { + int c = (int)(kps[idx].pt.x / std::max(1, width) * cols); + int r = (int)(kps[idx].pt.y / std::max(1, height) * rows); + c = std::min(std::max(c, 0), cols - 1); + r = std::min(std::max(r, 0), rows - 1); + cells[r * cols + c].push_back(idx); + } + for(auto &cell : cells) + std::sort(cell.begin(), cell.end(), + [&](int a, int b) { return kps[a].response > kps[b].response; }); + std::vector out; + out.reserve(target); + size_t maxlen = 0; + for(const auto &cell : cells) maxlen = std::max(maxlen, cell.size()); + for(size_t k = 0; k < maxlen && (int)out.size() < target; k++) + for(const auto &cell : cells) + if(k < cell.size()) + { + out.push_back(cell[k]); + if((int)out.size() >= target) break; + } + return out; +} + +// Inlier reprojection-error stats (mean/median/max, in pixels) of warp H mapping +// template `src` -> image `dst`, restricted to RANSAC inliers. Mirrors the +// "reproj mean/median/max" line of log_feature_init_stats(). +void reprojStats(const std::vector &src, const std::vector &dst, + const cv::Mat &inliers, const cv::Mat &H, dt_hdr_cv_feature_stats_t *stats) +{ + if(!stats) return; + stats->reproj_mean = stats->reproj_median = stats->reproj_max = -1.0; + if(H.empty() || src.empty()) return; + std::vector proj; + cv::perspectiveTransform(src, proj, H); + std::vector errs; + errs.reserve(src.size()); + for(size_t k = 0; k < src.size(); k++) + { + if(!inliers.empty() && inliers.at((int)k) == 0) continue; + const double dx = (double)proj[k].x - dst[k].x; + const double dy = (double)proj[k].y - dst[k].y; + errs.push_back(std::sqrt(dx * dx + dy * dy)); + } + if(errs.empty()) return; + double sum = 0.0, mx = 0.0; + for(double e : errs) { sum += e; mx = std::max(mx, e); } + std::sort(errs.begin(), errs.end()); + stats->reproj_mean = sum / (double)errs.size(); + stats->reproj_median = errs[errs.size() / 2]; + stats->reproj_max = mx; +} + +// Subsample matches to an even spatial distribution over a kSpatialGrid^2 grid of +// the image, so RANSAC is constrained by correspondences from the whole frame +// instead of one dense, well-exposed region. Mirrors _spatially_uniform_subsample(). +std::vector spatialSubsample(const std::vector &matches, + const std::vector &kp_img, + int width, int height, int target) +{ + if((int)matches.size() <= target) return matches; + const int gr = kSpatialGrid, gc = kSpatialGrid; + std::vector> cells(gr * gc); + for(int idx = 0; idx < (int)matches.size(); idx++) + { + const cv::Point2f &p = kp_img[matches[idx].queryIdx].pt; + int r = (int)(p.y / std::max(1, height) * gr); + int c = (int)(p.x / std::max(1, width) * gc); + r = std::min(std::max(r, 0), gr - 1); + c = std::min(std::max(c, 0), gc - 1); + cells[r * gc + c].push_back(idx); + } + for(auto &cell : cells) + std::sort(cell.begin(), cell.end(), + [&](int a, int b) { return matches[a].distance < matches[b].distance; }); + std::vector out; + out.reserve(target); + size_t maxlen = 0; + for(const auto &cell : cells) maxlen = std::max(maxlen, cell.size()); + for(size_t k = 0; k < maxlen && (int)out.size() < target; k++) + for(const auto &cell : cells) + { + if(k < cell.size()) + { + out.push_back(matches[cell[k]]); + if((int)out.size() >= target) break; + } + } + return out; +} + +// When RANSAC inliers cluster into <= kClusterDegradeMaxCells grid cells, an 8-DOF +// homography overfits scale/shear/perspective to a tiny region and extrapolates +// wildly. If the inlier displacements are consistent (low MAD), refit as a pure +// translation from their median. Mirrors the INLIER_CLUSTER_DEGRADE_* logic. +bool degradeClusteredToTranslation(const std::vector &src, + const std::vector &dst, + const cv::Mat &inliers, int width, int height, + double H[9]) +{ + const int gr = kSpatialGrid, gc = kSpatialGrid; + bool occ[kSpatialGrid * kSpatialGrid] = { false }; + int n_cells = 0; + std::vector dxs, dys; + for(size_t k = 0; k < src.size(); k++) + { + if(!inliers.empty() && inliers.at((int)k) == 0) continue; + int r = (int)(src[k].y / std::max(1, height) * gr); + int c = (int)(src[k].x / std::max(1, width) * gc); + r = std::min(std::max(r, 0), gr - 1); + c = std::min(std::max(c, 0), gc - 1); + if(!occ[r * gc + c]) { occ[r * gc + c] = true; n_cells++; } + dxs.push_back((double)dst[k].x - src[k].x); + dys.push_back((double)dst[k].y - src[k].y); + } + if((int)dxs.size() < 4 || n_cells > kClusterDegradeMaxCells) return false; + + std::vector sx = dxs, sy = dys; + std::sort(sx.begin(), sx.end()); + std::sort(sy.begin(), sy.end()); + const double tx = sx[sx.size() / 2], ty = sy[sy.size() / 2]; + std::vector ax, ay; + ax.reserve(dxs.size()); + ay.reserve(dys.size()); + for(double v : dxs) ax.push_back(std::fabs(v - tx)); + for(double v : dys) ay.push_back(std::fabs(v - ty)); + std::sort(ax.begin(), ax.end()); + std::sort(ay.begin(), ay.end()); + if(std::max(ax[ax.size() / 2], ay[ay.size() / 2]) > kClusterTranslationMaxMad) return false; + + H[0] = 1; H[1] = 0; H[2] = tx; + H[3] = 0; H[4] = 1; H[5] = ty; + H[6] = 0; H[7] = 0; H[8] = 1; + return true; +} + +/* ------------------------------------------------------------------------- * + * Optional debug-image dump. + * + * When the caller passes a directory (resolved in the C layer from the + * preference or the DT_HDR_DEBUG_IMAGE_DIR override), each aligned (moving) + * frame writes a numbered set of Netpbm images (PGM/PPM) to it: the CLAHE'd + * SIFT input, the detected keypoints, and the colour-coded match visualisation + * (green = inlier, red = outlier). Mirrors align_and_blend.py's "Saved feature + * debug visuals". Entirely diagnostic: no effect unless a directory is given. + * ------------------------------------------------------------------------- */ + +// Longest side of a dumped image; large proxies are scaled down to stay viewable +// and to keep the (uncompressed) files reasonable. +constexpr int kDebugMaxSide = 1600; + +// Write `img` (8-bit 1- or 3-channel, or float/other normalised to [0,255]) as a +// Netpbm image -- PGM (P5) for grayscale, PPM (P6) for colour -- scaling it down +// if larger than kDebugMaxSide. Netpbm needs no image library (keeps OpenCV to +// its minimal module set); the files open in any common image viewer. Never +// throws into the caller. +void writeDebugImage(const char *dir, int frame, const char *name, const cv::Mat &img) +{ + try + { + if(img.empty()) return; + cv::Mat u8; + if(img.depth() == CV_8U) + u8 = img; + else + { + cv::normalize(img, u8, 0, 255, cv::NORM_MINMAX); + u8.convertTo(u8, img.channels() == 3 ? CV_8UC3 : CV_8U); + } + const int side = std::max(u8.cols, u8.rows); + if(side > kDebugMaxSide) + { + const double f = (double)kDebugMaxSide / side; + cv::resize(u8, u8, cv::Size(), f, f, cv::INTER_AREA); + } + const int ch = u8.channels(); + char path[1024]; + snprintf(path, sizeof(path), "%s/hdr_frame%02d_%s.%s", dir, frame, name, + ch == 3 ? "ppm" : "pgm"); + FILE *fp = fopen(path, "wb"); + if(!fp) return; + fprintf(fp, "P%c\n%d %d\n255\n", ch == 3 ? '6' : '5', u8.cols, u8.rows); + if(ch == 3) + { + // OpenCV is BGR; Netpbm wants RGB, so swap per pixel, one row at a time. + std::vector row((size_t)u8.cols * 3); + for(int y = 0; y < u8.rows; y++) + { + const unsigned char *p = u8.ptr(y); + for(int x = 0; x < u8.cols; x++) + { + row[x * 3 + 0] = p[x * 3 + 2]; + row[x * 3 + 1] = p[x * 3 + 1]; + row[x * 3 + 2] = p[x * 3 + 0]; + } + fwrite(row.data(), 1, row.size(), fp); + } + } + else + { + for(int y = 0; y < u8.rows; y++) + fwrite(u8.ptr(y), 1, (size_t)u8.cols, fp); + } + fclose(fp); + dt_print(DT_DEBUG_HDR_MERGE, " debug image: %s", path); + } + catch(const cv::Exception &e) + { + dt_print(DT_DEBUG_HDR_MERGE, " debug image '%s' failed: %s", name, e.what()); + } +} + +// Single SIFT pass on one 8-bit proxy: detectAndCompute builds the Gaussian +// pyramid once (a separate detect + compute would build it twice), then the +// scale floor and spatial balance prune keypoints and their descriptor rows +// together via gatherKpDesc. On the raw HDR proxy the detection count is modest, +// so the descriptors computed for later-pruned keypoints cost little next to the +// saved second pyramid; with very high detection counts (e.g. CLAHE on) that +// overhead grows but the single pyramid still dominates. `out_img` receives the +// CLAHE'd image (the SIFT input) for the optional debug visuals. +void detectDescribe(const uint8_t *proxy, int width, int height, int balance_target, + double clahe_clip, std::vector &kp, cv::Mat &des, + cv::Mat &out_img, int &raw, int &after_floor) +{ + out_img = wrapU8(proxy, width, height).clone(); // clone: CLAHE writes in place + applyClahe(out_img, clahe_clip); + cv::Ptr sift = makeSift(); + sift->detectAndCompute(out_img, cv::noArray(), kp, des); + raw = (int)kp.size(); + gatherKpDesc(kp, des, scaleFloorKeep(kp)); + after_floor = (int)kp.size(); + // Spatially balance to a common budget so a feature-dense frame cannot flood + // the matcher with aliased candidates. Mirrors align_and_blend.py. + gatherKpDesc(kp, des, spatialBalanceKeep(kp, width, height, balance_target)); +} + +// Build a FLANN KD-tree matcher (5 randomized trees, 50 checks) -- the same +// parameters used for the per-frame matching. Factored out so the reference +// index and the per-frame moving index are built identically. +cv::Ptr makeFlann() +{ + return cv::makePtr(cv::makePtr(5), + cv::makePtr(50)); +} + +// Cached reference features: detected once per merge (dt_hdr_cv_features_create) +// and matched against by every moving frame, so the reference SIFT pass is not +// repeated for all N-1 moving frames. Holds the CLAHE'd reference image too, for +// the optional debug visuals. +struct RefFeatures +{ + cv::Mat image; // CLAHE'd 8-bit reference proxy (debug visuals) + std::vector kp; // balanced reference keypoints + cv::Mat des; // matching descriptor rows + int kp_raw = 0; // detections before the scale floor + int kp_floor = 0; // kept after the scale floor (balance input) + // FLANN index over `des`, trained once. Every moving frame's image->template + // match reuses this instead of rebuilding the reference KD-tree N-1 times. + cv::Ptr matcher; +}; + +} // namespace + +extern "C" { + +void *dt_hdr_cv_features_create(const uint8_t *proxy, int width, int height, + int balance_target, double clahe_clip) +{ + if(balance_target <= 0) balance_target = kSiftSpatialBalanceTarget; + try + { + RefFeatures *f = new RefFeatures(); + detectDescribe(proxy, width, height, balance_target, clahe_clip, + f->kp, f->des, f->image, f->kp_raw, f->kp_floor); + // Train the reference-side FLANN index once; every moving frame reuses it. + if(!f->des.empty()) + { + f->matcher = makeFlann(); + f->matcher->add(std::vector{ f->des }); + f->matcher->train(); + } + return f; + } + catch(const std::exception &e) + { + dt_print(DT_DEBUG_HDR_MERGE, "HDR ref feature precompute raised: %s", e.what()); + return nullptr; + } +} + +void dt_hdr_cv_features_destroy(void *features) +{ + delete static_cast(features); +} + +int dt_hdr_cv_feature_homography(const void *ref_features, const uint8_t *img, + int width, int height, int balance_target, + double clahe_clip, const char *debug_dir, + int frame_index, double H[9], + dt_hdr_cv_feature_stats_t *stats) +{ + if(balance_target <= 0) balance_target = kSiftSpatialBalanceTarget; + // Identity by default. + H[0] = 1; H[1] = 0; H[2] = 0; H[3] = 0; H[4] = 1; H[5] = 0; H[6] = 0; H[7] = 0; H[8] = 1; + if(stats) + { + *stats = dt_hdr_cv_feature_stats_t{}; + stats->reproj_mean = stats->reproj_median = stats->reproj_max = -1.0; + } + + const RefFeatures *ref = static_cast(ref_features); + if(!ref || ref->kp.empty() || ref->des.empty()) return 0; + + try + { + // The reference (template) comes from the per-merge cache; only the moving + // frame is detected here. Bind t / kp_t / des_t to the cache so the matching + // and debug code below is identical to the symmetric (uncached) version. + const cv::Mat &t = ref->image; + const std::vector &kp_t = ref->kp; + const cv::Mat &des_t = ref->des; + + cv::Mat mov_img, des_i; + std::vector kp_i; + int ki_raw = 0, ki_floor = 0; + detectDescribe(img, width, height, balance_target, clahe_clip, + kp_i, des_i, mov_img, ki_raw, ki_floor); + + if(stats) { stats->kp_template_raw = ref->kp_raw; stats->kp_image_raw = ki_raw; } + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT scale floor: kept template %d/%d, image %d/%d keypoints with kp.size >= %.1fpx", + ref->kp_floor, ref->kp_raw, ki_floor, ki_raw, kSiftMinScalePx); + if((int)kp_i.size() < ki_floor) + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT spatial balance: image %d->%d keypoints", ki_floor, (int)kp_i.size()); + + if(stats) { stats->kp_template = (int)kp_t.size(); stats->kp_image = (int)kp_i.size(); } + if(des_t.empty() || des_i.empty() + || (int)kp_t.size() < kMinFeatureKeypoints || (int)kp_i.size() < kMinFeatureKeypoints) + { + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT init: insufficient features (need >= %d), using identity initial warp", + kMinFeatureKeypoints); + return 0; + } + + // FLANN kNN match with Lowe's ratio test, in both directions, then keep + // only the mutually-best correspondences (FEATURE_REQUIRE_MUTUAL_CONSISTENCY). + // it (image->template): query = moving (des_i), train = the cached + // reference index -> queryIdx in image, trainIdx in template (the Python + // convention). Reuses ref->matcher, so the reference KD-tree is built + // once per merge, not once per moving frame. + // ti (template->image): query = reference (des_t), train = a fresh index + // over the moving descriptors (rebuilt every frame, unavoidable). + // The two directions are independent, so they run in two parallel sections; + // the team is capped to 2 threads (only two sections do work). + std::vector> knn_it, knn_ti; + bool match_failed = false; + std::string match_err; + auto capture_err = [&](const cv::Exception &e) + { +#ifdef _OPENMP +#pragma omp critical(hdr_match_err) +#endif + { + match_failed = true; + match_err = e.what(); + } + }; + +#ifdef _OPENMP +#pragma omp parallel sections num_threads(2) +#endif + { +#ifdef _OPENMP +#pragma omp section +#endif + { + try + { + // Reuse the trained reference index; fall back to a fresh one only if + // the cache is somehow absent (defensive -- des_t is non-empty here). + if(ref->matcher) + ref->matcher->knnMatch(des_i, knn_it, 2); + else + makeFlann()->knnMatch(des_i, des_t, knn_it, 2); + } + catch(const cv::Exception &e) { capture_err(e); } + } +#ifdef _OPENMP +#pragma omp section +#endif + { + try + { + makeFlann()->knnMatch(des_t, des_i, knn_ti, 2); + } + catch(const cv::Exception &e) { capture_err(e); } + } + } + if(match_failed) + { + dt_print(DT_DEBUG_HDR_MERGE, "SIFT init match raised: %s", match_err.c_str()); + return 0; + } + + // Reverse best (template kp -> image kp) after the ratio test. + std::vector rev_best(kp_t.size(), -1); + for(const auto &m : knn_ti) + if(m.size() >= 2 && m[0].distance < kRatioThreshold * m[1].distance) + rev_best[m[0].queryIdx] = m[0].trainIdx; + + int ratio_count = 0; + std::vector good; + good.reserve(knn_it.size()); + for(const auto &m : knn_it) + { + if(m.size() < 2) continue; + if(m[0].distance >= kRatioThreshold * m[1].distance) continue; // Lowe ratio + ratio_count++; + // Mutual-consistency: the template kp must point back to this image kp. + if(rev_best[m[0].trainIdx] != m[0].queryIdx) continue; + good.push_back(m[0]); + } + if(stats) { stats->ratio_matches = ratio_count; stats->good_matches = (int)good.size(); } + if((int)good.size() > 0) + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT init: %d/%d ratio-test matches are mutual-consistent", + (int)good.size(), ratio_count); + if((int)good.size() < kMinGoodMatches) + { + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT init: only %d good matches, using identity initial warp", + (int)good.size()); + return 0; + } + + // Subsample matches to an even spatial distribution, so RANSAC is constrained + // by correspondences across the whole frame rather than one dense region. + { + const int before = (int)good.size(); + good = spatialSubsample(good, kp_i, width, height, kMaxMatchesForRansac); + if((int)good.size() < before) + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT init: spatially-subsampled to %d matches (%dx%d grid) for RANSAC", + (int)good.size(), kSpatialGrid, kSpatialGrid); + } + + // src = template coords, dst = image coords => warp maps template->image, + // consistent with WARP_INVERSE_MAP sampling. + std::vector src, dst; + src.reserve(good.size()); + dst.reserve(good.size()); + for(const auto &m : good) + { + src.push_back(kp_t[m.trainIdx].pt); + dst.push_back(kp_i[m.queryIdx].pt); + } + + cv::Mat inliers; + cv::Mat Hh = cv::findHomography(src, dst, cv::RANSAC, kRansacReprojThreshold, inliers, + kRansacMaxIters, kRansacConfidence); + int n_in = inliers.empty() ? 0 : cv::countNonZero(inliers); + const int homog_inliers = n_in; + const double homog_ratio = good.empty() ? 0.0 : (double)n_in / (double)good.size(); + const bool weak = Hh.empty() || n_in < kHomographyMinInliers + || homog_ratio < kHomographyMinInlierRatio; + + // Choose the model: homography, or estimateAffine2D fallback on weak support. + cv::Mat chosen, chosen_inl; + bool used_affine = false; + if(weak) + { + cv::Mat ainl; + cv::Mat A = cv::estimateAffine2D(src, dst, ainl, cv::RANSAC, kRansacReprojThreshold, + kRansacMaxIters, kRansacConfidence); + if(!A.empty()) + { + chosen = cv::Mat::eye(3, 3, CV_64F); + A.copyTo(chosen(cv::Rect(0, 0, 3, 2))); + chosen_inl = ainl; + used_affine = true; + n_in = ainl.empty() ? 0 : cv::countNonZero(ainl); + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT init: weak homography support (%d/%d inliers, %.0f%%); using affine fallback", + homog_inliers, (int)good.size(), homog_ratio * 100.0); + } + } + if(chosen.empty()) + { + if(Hh.empty()) return 0; + chosen = Hh; + chosen_inl = inliers; + } + fromMat3x3(chosen, H); + + // Cluster-degradation: a homography (or affine) fit to inliers crammed into a + // couple of grid cells overfits scale/shear and extrapolates wildly across the + // frame; if the inlier displacements agree, refit as a robust translation. + double Htrans[9]; + cv::Mat reproj_H = chosen; + bool used_translation = false; + if(degradeClusteredToTranslation(src, dst, chosen_inl, width, height, Htrans)) + { + for(int k = 0; k < 9; k++) H[k] = Htrans[k]; + used_translation = true; + reproj_H = (cv::Mat_(3, 3) << 1, 0, Htrans[2], 0, 1, Htrans[5], 0, 0, 1); + dt_print(DT_DEBUG_HDR_MERGE, + "SIFT init: inliers cluster in <= %d cells; refitting as translation-only " + "tx=%.2f ty=%.2f", + kClusterDegradeMaxCells, Htrans[2], Htrans[5]); + } + + if(stats) + { + stats->inliers = n_in; + stats->used_affine = used_affine ? 1 : 0; + stats->used_translation = used_translation ? 1 : 0; + } + reprojStats(src, dst, chosen_inl, reproj_H, stats); + + // Optional debug visuals: the (CLAHE'd) SIFT input, the kept keypoints, and + // the colour-coded matches -- the quickest way to see *why* the match set + // consensed where it did (e.g. a periodic structure matching one period off). + // `debug_dir` / `frame_index` are supplied by the (per-merge) caller. + if(debug_dir && debug_dir[0]) + { + const char *dir = debug_dir; + const int fr = frame_index; + // The reference (template) visuals are identical for every moving frame, + // so write them once (on the first aligned frame) instead of per frame. + if(frame_index <= 1) + { + writeDebugImage(dir, fr, "1_template_sift_input", t); + cv::Mat kt_vis; + cv::drawKeypoints(t, kp_t, kt_vis, cv::Scalar::all(-1), + cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS); + writeDebugImage(dir, fr, "3_template_keypoints", kt_vis); + } + writeDebugImage(dir, fr, "2_image_sift_input", mov_img); + cv::Mat ki_vis; + cv::drawKeypoints(mov_img, kp_i, ki_vis, cv::Scalar::all(-1), + cv::DrawMatchesFlags::DRAW_RICH_KEYPOINTS); + writeDebugImage(dir, fr, "4_image_keypoints", ki_vis); + // All matches, colour-coded: red = RANSAC outlier, green = inlier. `good` + // matches have queryIdx -> kp_i (image), trainIdx -> kp_t (template), and + // chosen_inl is the per-match inlier mask over `good`. Draw the outliers + // first, then the inliers over the top so the consensus stands out. + std::vector inlier_mask(good.size(), 0), outlier_mask(good.size(), 1); + for(size_t k = 0; k < good.size(); k++) + { + const bool in = !chosen_inl.empty() && chosen_inl.at((int)k); + inlier_mask[k] = in ? 1 : 0; + outlier_mask[k] = in ? 0 : 1; + } + const cv::Scalar red(0, 0, 255), green(0, 255, 0); // BGR + cv::Mat match_vis; + cv::drawMatches(mov_img, kp_i, t, kp_t, good, match_vis, red, red, outlier_mask, + cv::DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS); + cv::drawMatches(mov_img, kp_i, t, kp_t, good, match_vis, green, green, inlier_mask, + cv::DrawMatchesFlags::DRAW_OVER_OUTIMG + | cv::DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS); + writeDebugImage(dir, fr, "5_all_matches_green_inlier_red_outlier", match_vis); + } + return n_in; + } + catch(const cv::Exception &e) + { + dt_print(DT_DEBUG_HDR_MERGE, "SIFT init raised: %s", e.what()); + return 0; + } +} + +int dt_hdr_cv_count_features(const uint8_t *luma, int width, int height, int probe_dim) +{ + try + { + cv::Mat in = wrapU8(luma, width, height); + cv::Mat resized; + const int max_side = std::max(width, height); + if(probe_dim > 0 && max_side > probe_dim) + { + const double s = (double)probe_dim / (double)max_side; + cv::resize(in, resized, cv::Size(), s, s, cv::INTER_AREA); + } + // The probe ranks the gamma-encoded proxies as detection sees them (CLAHE + // off by default); the ranking only needs to be self-consistent. Detect + // directly on the wrapped (or resized) buffer -- SIFT does not modify its + // input, so no defensive copy is needed. + const cv::Mat &probe = resized.empty() ? in : resized; + + cv::Ptr sift = makeSift(); + std::vector kps; + sift->detect(probe, kps); + return (int)kps.size(); + } + catch(const cv::Exception &) + { + return 0; + } +} + +} // extern "C" + +#endif // HAVE_OPENCV diff --git a/src/control/jobs/control_jobs.c b/src/control/jobs/control_jobs.c index 9c8fc02e4a63..f0102cb265fa 100644 --- a/src/control/jobs/control_jobs.c +++ b/src/control/jobs/control_jobs.c @@ -23,6 +23,7 @@ #include "common/exif.h" #include "common/film.h" #include "common/gpx.h" +#include "common/hdr_alignment.h" #include "common/history.h" #include "common/history_snapshot.h" #include "common/image.h" @@ -404,7 +405,7 @@ static int32_t _control_write_sidecar_files_job_run(dt_job_t *job) typedef struct dt_control_merge_hdr_t { - uint32_t first_imgid; + dt_imgid_t first_imgid; uint32_t first_filter; uint8_t first_xtrans[6][6]; @@ -420,6 +421,18 @@ typedef struct dt_control_merge_hdr_t float adobe_XYZ_to_CAM[4][3]; char camera_makermodel[128]; + // auto-alignment of exposure brackets (handheld / shaky tripod). When + // enabled, the first frame seeds the reference and every later frame is + // registered onto it before accumulation. NULL/disabled => legacy behavior. + dt_hdr_align_t *align; + gboolean align_enabled; + float *aligned_buf; // scratch holding the warped moving frame (wd*ht) + + // auto-reference pre-pass: pick the richest-feature frame as the reference. + gboolean probe_mode; // set during the feature-probe pre-pass + int probe_best_count; // best SIFT feature count seen so far + dt_imgid_t probe_best_imgid; // imgid of the richest frame + // 0 - ok; 1 - errors, abort gboolean abort; } dt_control_merge_hdr_t; @@ -486,6 +499,27 @@ static int _control_merge_hdr_process(dt_imageio_module_data_t *datai, const dt_image_t image = *img; dt_image_cache_read_release(img); + // Auto-reference pre-pass: just probe each frame's feature richness and keep + // the richest as the reference; do not accumulate anything. + if(d->probe_mode) + { + if(image.buf_dsc.filters != 0u + && image.buf_dsc.channels == 1 + && image.buf_dsc.datatype == TYPE_UINT16) + { + const int n = dt_hdr_alignment_probe_features((const float *)ivoid, + datai->width, datai->height); + dt_print(DT_DEBUG_HDR_MERGE, + "[hdr merge] auto-reference probe: imgid %d -> %d features", imgid, n); + if(n > d->probe_best_count) + { + d->probe_best_count = n; + d->probe_best_imgid = imgid; + } + } + return 0; + } + if(!d->pixels) { d->first_imgid = imgid; @@ -532,6 +566,14 @@ static int _control_merge_hdr_process(dt_imageio_module_data_t *datai, return 1; } + // Authoritative per-frame guard on the *decoded* buffer. + // _control_merge_hdr_validate() gives a friendly early failure, but it can only + // see the cached image flags, which are derived from metadata and can lag what + // rawspeed actually decodes here: DT_IMAGE_MONOCHROME / DT_IMAGE_HDR and the CFA + // layout are only settled during the load. A monochrome raw in particular + // carries DT_IMAGE_RAW yet decodes to filters == 0. The accumulation loop below + // indexes the raw CFA mosaic and the result is written as a CFA DNG, so a + // non-CFA / non-uint16 buffer here must abort rather than corrupt the output. if(image.buf_dsc.filters == 0u || image.buf_dsc.channels != 1 || image.buf_dsc.datatype != TYPE_UINT16) @@ -540,9 +582,11 @@ static int _control_merge_hdr_process(dt_imageio_module_data_t *datai, d->abort = TRUE; return 1; } - else if(datai->width != d->wd - || datai->height != d->ht - || d->orientation != image.orientation) + // Re-check size / orientation too: the accumulation loop indexes the current + // frame's buffer with the first frame's stride. + if(datai->width != d->wd + || datai->height != d->ht + || d->orientation != image.orientation) { dt_control_log(_("images have to be of same size and orientation!")); d->abort = TRUE; @@ -562,13 +606,39 @@ static int _control_merge_hdr_process(dt_imageio_module_data_t *datai, const float photoncnt = 100.0f * aperture * exp / iso; float saturation = 1.0f; d->whitelevel = fmaxf(d->whitelevel, saturation * cal); + + // Auto-align this frame onto the reference before accumulation. The first + // frame becomes the reference; every later frame is registered onto it and + // accumulated from a warped scratch buffer. On any failure (or when + // alignment is disabled / built without OpenCV) we fall back to the raw + // unaligned mosaic, i.e. the legacy behavior. + const float *in_buf = (const float *)ivoid; + if(d->align_enabled && d->align) + { + dt_print(DT_DEBUG_HDR_MERGE, "Processing image %d/%d: imgid %d", num, total, imgid); + if(imgid == d->first_imgid) + { + dt_print(DT_DEBUG_HDR_MERGE, " Reference image (no alignment needed)"); + if(dt_hdr_alignment_set_reference(d->align, (const float *)ivoid, d->wd, d->ht, + d->first_filter, d->first_xtrans)) + d->aligned_buf = dt_alloc_align_float((size_t)d->wd * d->ht); + } + else if(d->aligned_buf) + { + if(dt_hdr_alignment_align_frame(d->align, (const float *)ivoid, d->aligned_buf, + d->wd, d->ht, d->first_filter, d->first_xtrans, NULL)) + in_buf = d->aligned_buf; + } + } + DT_OMP_FOR(collapse(2)) for(int y = 0; y < d->ht; y++) for(int x = 0; x < d->wd; x++) { // read unclamped raw value with subtracted black and rescaled - // to 1.0 saturation. this is the output of the rawprepare iop. - const float in = ((float *)ivoid)[x + d->wd * y]; + // to 1.0 saturation. this is the output of the rawprepare iop + // (optionally registered onto the reference frame, see in_buf above). + const float in = in_buf[x + d->wd * y]; // weights based on siggraph 12 poster zijian zhu, zhengguo li, // susanto rahardja, pasi fraenti 2d denoising factor for high // dynamic range imaging @@ -588,8 +658,8 @@ static int _control_merge_hdr_process(dt_imageio_module_data_t *datai, for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++) { - M = MAX(M, ((float *)ivoid)[xx + i + d->wd * (yy + j)]); - m = MIN(m, ((float *)ivoid)[xx + i + d->wd * (yy + j)]); + M = MAX(M, in_buf[xx + i + d->wd * (yy + j)]); + m = MIN(m, in_buf[xx + i + d->wd * (yy + j)]); } // move envelope a little to allow non-zero weight even for // clipped regions. this is because even if the 2x2 block is @@ -631,17 +701,102 @@ static int _control_merge_hdr_process(dt_imageio_module_data_t *datai, return 0; } +// Validate that the selected images form a consistent raw exposure bracket. +// Reject non-raw inputs, already-merged HDR DNGs, or images with mismatched +// resolution/orientation before expensive decode/align work. Returns TRUE if +// mergeable, otherwise logs the odd file and returns FALSE. +static gboolean _control_merge_hdr_validate(GList *images) +{ + gboolean have_ref = FALSE; + int32_t ref_width = 0, ref_height = 0; + dt_image_orientation_t ref_orientation = ORIENTATION_NONE; + + for(GList *l = images; l; l = g_list_next(l)) + { + const dt_imgid_t imgid = GPOINTER_TO_INT(l->data); + const dt_image_t *img = dt_image_cache_get(imgid, 'r'); + if(!img) continue; + + const int32_t flags = img->flags; + const int32_t width = img->width; + const int32_t height = img->height; + const dt_image_orientation_t orientation = img->orientation; + char filename[DT_MAX_FILENAME_LEN]; + g_strlcpy(filename, img->filename, sizeof(filename)); + dt_image_cache_read_release(img); + + // exposure bracketing operates on the raw CFA mosaic; reject anything that + // is not raw, reject already-merged HDR DNGs (floating point raws carry the + // DT_IMAGE_HDR flag) which look like raws but are not exposure brackets, and + // reject monochrome raws (no CFA -> filters == 0, cannot be merged as a CFA + // mosaic). These flags are extension/metadata derived and may not yet be set + // for a not-yet-decoded frame; _control_merge_hdr_process() re-checks the + // decoded buffer as the authoritative guard. + if(!(flags & DT_IMAGE_RAW) || (flags & DT_IMAGE_HDR) || (flags & DT_IMAGE_MONOCHROME)) + { + dt_control_log(_("HDR merge aborted: `%s' is not a raw exposure" + " (already merged HDR, monochrome or non-raw image)"), filename); + return FALSE; + } + + if(!have_ref) + { + ref_width = width; + ref_height = height; + ref_orientation = orientation; + have_ref = TRUE; + } + else if(width != ref_width + || height != ref_height + || orientation != ref_orientation) + { + dt_control_log(_("HDR merge aborted: `%s' has a different size" + " or orientation"), filename); + return FALSE; + } + } + return TRUE; +} + static int32_t _control_merge_hdr_job_run(dt_job_t *job) { dt_control_image_enumerator_t *params = dt_control_job_get_params(job); GList *t = params->index; const guint total = g_list_length(t); double fraction = 0; + + // reject mismatched selections (already-merged HDR DNGs, differing resolution + // or orientation) before doing any expensive decode / alignment work + if(!_control_merge_hdr_validate(t)) + return 0; + dt_control_job_set_progress_message(job, ngettext("merging %d image", "merging %d images", total), total); dt_control_merge_hdr_t d = (dt_control_merge_hdr_t){.epsw = 1e-8f, .abort = FALSE }; + // Opt-in auto-alignment of the brackets (only when built with OpenCV). +#ifdef HAVE_OPENCV + d.align_enabled = dt_conf_get_bool("plugins/lighttable/hdr_merge_auto_align"); + if(d.align_enabled) + { + // Seed the runtime parameters from preferences (falling back to the built-in + // defaults for any value the user has not changed). + dt_hdr_align_params_t align_params; + dt_hdr_alignment_default_params(&align_params); + align_params.proxy_scale = dt_conf_get_float("plugins/lighttable/hdr_merge_proxy_scale"); + align_params.feature_gamma = dt_conf_get_float("plugins/lighttable/hdr_merge_feature_gamma"); + align_params.clahe_clip = dt_conf_get_float("plugins/lighttable/hdr_merge_clahe_clip"); + align_params.sift_keypoints = dt_conf_get_int("plugins/lighttable/hdr_merge_sift_keypoints"); + align_params.debug_images = dt_conf_get_bool("plugins/lighttable/hdr_merge_debug_images"); + d.align = dt_hdr_alignment_new(&align_params); + if(d.align) + dt_control_log(_("merging HDR: auto-aligning exposure brackets")); + else + d.align_enabled = FALSE; + } +#endif + dt_imageio_module_format_t buf = (dt_imageio_module_format_t) {.mime = _control_merge_hdr_mime, .levels = _control_merge_hdr_levels, @@ -651,6 +806,51 @@ static int32_t _control_merge_hdr_job_run(dt_job_t *job) dt_control_merge_hdr_format_t dat = (dt_control_merge_hdr_format_t){.parent = { 0 }, .d = &d }; + // Auto-reference pre-pass: probe every frame's feature richness and move the + // richest one to the front, so it becomes both the geometry seed and the + // alignment reference (mirrors align_and_blend.py --auto-reference). Costs one + // extra decode pass, hence opt-in. +#ifdef HAVE_OPENCV + if(d.align_enabled && total > 1 + && dt_conf_get_bool("plugins/lighttable/hdr_merge_auto_reference")) + { + dt_control_job_set_progress_message(job, _("examining frames for HDR auto-reference")); + d.probe_mode = TRUE; + d.probe_best_count = -1; + d.probe_best_imgid = NO_IMGID; + int pnum = 1; + for(GList *p = t; p && !d.abort + && dt_control_job_get_state(job) != DT_JOB_STATE_CANCELLED; + p = g_list_next(p)) + { + const dt_imgid_t pimgid = GPOINTER_TO_INT(p->data); + dt_imageio_export_with_flags(pimgid, "unused", &buf, (dt_imageio_module_data_t *)&dat, + TRUE, FALSE, TRUE, TRUE, FALSE, 1.0, + FALSE, "pre:rawprepare", FALSE, + FALSE, DT_COLORSPACE_NONE, NULL, DT_INTENT_LAST, NULL, + NULL, pnum, total, NULL, -1); + // the probe pre-pass is a full extra decode sweep; advance the bar over it + // so a large bracket does not look hung, and honor a cancel between frames + dt_control_job_set_progress(job, (double)pnum / (double)total); + pnum++; + } + d.probe_mode = FALSE; + if(d.abort || dt_control_job_get_state(job) == DT_JOB_STATE_CANCELLED) + goto end; + if(dt_is_valid_imgid(d.probe_best_imgid)) + { + params->index = g_list_remove(params->index, GINT_TO_POINTER(d.probe_best_imgid)); + params->index = g_list_prepend(params->index, GINT_TO_POINTER(d.probe_best_imgid)); + t = params->index; + dt_print(DT_DEBUG_HDR_MERGE, + "[hdr merge] auto-reference: imgid %d selected (%d features)", + d.probe_best_imgid, d.probe_best_count); + } + dt_control_job_set_progress_message(job, ngettext("merging %d image", + "merging %d images", total), total); + } +#endif + int num = 1; while(t) { @@ -692,7 +892,15 @@ static int32_t _control_merge_hdr_job_run(dt_job_t *job) const int exif_len = dt_exif_read_blob(&exif, pathname, d.first_imgid, FALSE, d.wd, d.ht, TRUE); char *c = pathname + strlen(pathname); while(*c != '.' && c > pathname) c--; - g_strlcpy(c, "-hdr.dng", sizeof(pathname) - (c - pathname)); + const size_t free_space = sizeof(pathname) - (c - pathname); + + // don't overwrite an existing merge: append a numeric suffix + // (-hdr_01.dng, -hdr_02.dng, ...) until we find a free filename, + // matching the behavior of export and raw denoise + g_strlcpy(c, "-hdr.dng", free_space); + for(int seq = 1; g_file_test(pathname, G_FILE_TEST_EXISTS); seq++) + snprintf(c, free_space, "-hdr_%02d.dng", seq); + dt_imageio_dng_write_float(pathname, d.pixels, d.wd, @@ -728,6 +936,8 @@ static int32_t _control_merge_hdr_job_run(dt_job_t *job) end: free(d.pixels); free(d.weight); + if(d.aligned_buf) dt_free_align(d.aligned_buf); + if(d.align) dt_hdr_alignment_free(d.align); return 0; } diff --git a/src/tests/unittests/CMakeLists.txt b/src/tests/unittests/CMakeLists.txt index 239c0a0ebd3f..44ee0c4a7622 100644 --- a/src/tests/unittests/CMakeLists.txt +++ b/src/tests/unittests/CMakeLists.txt @@ -13,3 +13,25 @@ if(WIN32) target_link_libraries(test_sample PRIVATE lib_darktable) _copy_required_library(test_sample lib_darktable) endif(WIN32) + +# HDR exposure-bracket auto-alignment. Calls the public dt_hdr_alignment_* API +# exported from lib_darktable (which carries the OpenCV-linked objects); the test +# itself adapts at runtime to OpenCV-present / -absent builds. +add_cmocka_test(test_hdr_alignment + SOURCES test_hdr_alignment.c + LINK_LIBRARIES lib_darktable cmocka) + +if(WIN32) + _copy_required_library(test_hdr_alignment lib_darktable) +endif(WIN32) + +# White-box tests for the internal (static) helpers of the HDR alignment code. +# Following test_filmicrgb.c, the .c is #included directly for static access; +# lib_darktable supplies what it links against (OpenCV seam, glib, ...). +add_cmocka_test(test_hdr_alignment_internal + SOURCES test_hdr_alignment_internal.c + LINK_LIBRARIES lib_darktable cmocka) + +if(WIN32) + _copy_required_library(test_hdr_alignment_internal lib_darktable) +endif(WIN32) diff --git a/src/tests/unittests/test_hdr_alignment.c b/src/tests/unittests/test_hdr_alignment.c new file mode 100644 index 000000000000..bb28b838e80a --- /dev/null +++ b/src/tests/unittests/test_hdr_alignment.c @@ -0,0 +1,358 @@ +/* + This file is part of darktable, + Copyright (C) 2026 darktable developers. + + darktable 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. + + darktable 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 darktable. If not, see . +*/ + +/* Unit tests for the HDR exposure-bracket auto-alignment (common/hdr_alignment). + * + * The tests synthesise a Bayer mosaic, warp it by a known homography, and check + * that the aligner recovers the motion and reduces the misalignment error; they + * also check that the auto-reference probe ranks frames by feature richness. + * + * The assertions are robust to a build *without* OpenCV: in that case the public + * API reports DT_HDR_ALIGN_DISABLED / a zero probe count and the tests verify the + * documented graceful fall-through instead. + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "common/hdr_alignment.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#ifdef _WIN32 +#include "win/main_wrapper.h" +#endif + +#define TEST_W 512 +#define TEST_H 512 +// RGGB-style Bayer encoding compatible with the FC() formula; the scene is +// grayscale so the exact phase only exercises the same-color sampler. +#define TEST_FILTERS 0x94949494u + +static unsigned int _rng = 12345u; +static float _frand(void) +{ + _rng = _rng * 1664525u + 1013904223u; + return (float)((_rng >> 8) & 0xFFFFFF) / (float)0xFFFFFF; +} + +// Deterministic textured scene: gentle gradients + scattered bright blobs and +// rectangles, which give SIFT plenty of distinctive keypoints. +static void _make_scene(float *s, int rich) +{ + _rng = 12345u; + for(int y = 0; y < TEST_H; y++) + for(int x = 0; x < TEST_W; x++) + s[y * TEST_W + x] = 0.30f + 0.08f * sinf(x * 0.011f + y * 0.004f) + + 0.08f * cosf(y * 0.013f - x * 0.006f); + if(!rich) return; + for(int k = 0; k < 120; k++) + { + int rx = (int)(_frand() * TEST_W), ry = (int)(_frand() * TEST_H); + int rw = 12 + (int)(_frand() * 26), rh = 12 + (int)(_frand() * 26); + float ri = 0.5f + 0.4f * _frand(); + for(int y = ry; y < ry + rh && y < TEST_H; y++) + for(int x = rx; x < rx + rw && x < TEST_W; x++) s[y * TEST_W + x] = ri; + } + for(int k = 0; k < 240; k++) + { + float cx = _frand() * TEST_W, cy = _frand() * TEST_H, r = 10.0f + _frand() * 18.0f; + float inten = 0.4f + 0.5f * _frand(); + for(int y = (int)(cy - r - 1); y <= (int)(cy + r + 1); y++) + for(int x = (int)(cx - r - 1); x <= (int)(cx + r + 1); x++) + { + if(x < 0 || x >= TEST_W || y < 0 || y >= TEST_H) continue; + const float d = hypotf(x - cx, y - cy); + if(d <= r) + { + const float v = inten * (1.0f - d / r); + if(v > s[y * TEST_W + x]) s[y * TEST_W + x] = v; + } + } + } +} + +static float _bilin(const float *im, float x, float y) +{ + if(x < 0) x = 0; + if(y < 0) y = 0; + if(x > TEST_W - 1) x = TEST_W - 1; + if(y > TEST_H - 1) y = TEST_H - 1; + const int x0 = (int)x, y0 = (int)y; + const int x1 = x0 < TEST_W - 1 ? x0 + 1 : x0, y1 = y0 < TEST_H - 1 ? y0 + 1 : y0; + const float fx = x - x0, fy = y - y0; + return (1 - fx) * (1 - fy) * im[y0 * TEST_W + x0] + fx * (1 - fy) * im[y0 * TEST_W + x1] + + (1 - fx) * fy * im[y1 * TEST_W + x0] + fx * fy * im[y1 * TEST_W + x1]; +} + +// Alignment of a frame displaced by a known homography must either reduce the +// misalignment error (OpenCV build) or pass the frame through unchanged +// (no-OpenCV build). +static void test_hdr_align_reduces_misalignment(void **state) +{ + (void)state; + float *ref = malloc(sizeof(float) * TEST_W * TEST_H); + float *mov = malloc(sizeof(float) * TEST_W * TEST_H); + float *out = malloc(sizeof(float) * TEST_W * TEST_H); + assert_non_null(ref); + assert_non_null(mov); + assert_non_null(out); + _make_scene(ref, 1); + + // Ground-truth ref->moving homography: rotation about the centre + translation. + const double th = 1.2 * M_PI / 180.0, cx = TEST_W / 2.0, cy = TEST_H / 2.0; + const double c = cos(th), s = sin(th), tx = 5.0, ty = -3.0; + const double H0 = c, H1 = -s, H2 = cx - c * cx + s * cy + tx; + const double H3 = s, H4 = c, H5 = cy - s * cx - c * cy + ty; + const double det = H0 * H4 - H1 * H3; + // inverse affine, to synthesise moving(u) = ref(Hinv * u) + const double i0 = H4 / det, i1 = -H1 / det, i3 = -H3 / det, i4 = H0 / det; + const double i2 = -(i0 * H2 + i1 * H5), i5 = -(i3 * H2 + i4 * H5); + for(int y = 0; y < TEST_H; y++) + for(int x = 0; x < TEST_W; x++) + mov[y * TEST_W + x] = _bilin(ref, (float)(i0 * x + i1 * y + i2), (float)(i3 * x + i4 * y + i5)); + + dt_hdr_align_t *a = dt_hdr_alignment_new(NULL); + assert_non_null(a); + assert_true(dt_hdr_alignment_set_reference(a, ref, TEST_W, TEST_H, TEST_FILTERS, NULL)); + + dt_hdr_align_result_t info; + const gboolean ok = + dt_hdr_alignment_align_frame(a, mov, out, TEST_W, TEST_H, TEST_FILTERS, NULL, &info); + + if(info.status == DT_HDR_ALIGN_DISABLED) + { + // Built without OpenCV: the frame must be passed through untouched. + assert_false(ok); + assert_memory_equal(out, mov, sizeof(float) * TEST_W * TEST_H); + } + else + { + // Built with OpenCV: alignment must meaningfully reduce the misalignment. + double e_before = 0.0, e_after = 0.0; + for(int y = 48; y < TEST_H - 48; y++) + for(int x = 48; x < TEST_W - 48; x++) + { + e_before += fabs(ref[y * TEST_W + x] - mov[y * TEST_W + x]); + e_after += fabs(ref[y * TEST_W + x] - out[y * TEST_W + x]); + } + assert_true(ok); + assert_true(e_after < 0.7 * e_before); + } + + dt_hdr_alignment_free(a); + free(ref); + free(mov); + free(out); +} + +// CFA color of photosite (y, x) for TEST_FILTERS (RGGB), matching the FC() +// formula in hdr_alignment.c. 0 = R, 1 = G, 2 = B. +static int _fc_rggb(int y, int x) +{ + static const int rggb[4] = { 0, 1, 1, 2 }; // (0,0)R (0,1)G (1,0)G (1,1)B + return rggb[((y & 1) << 1) | (x & 1)]; +} + +// Headline property: aligning a genuinely CFA-*modulated* mosaic must warp each +// colour channel without cross-talk. We build a real Bayer mosaic (per-channel +// gains, so adjacent photosites differ by colour), displace it by a known +// homography, and require the aligner to recover the reference per photosite. +// A colour-blind (bilinear) warp would blend the strongly-different R/G/B taps +// and fail to reduce the error; the same-colour resampler must not. +static void test_hdr_align_cfa_modulated(void **state) +{ + (void)state; + float *scene = malloc(sizeof(float) * TEST_W * TEST_H); + float *ref = malloc(sizeof(float) * TEST_W * TEST_H); + float *mov = malloc(sizeof(float) * TEST_W * TEST_H); + float *out = malloc(sizeof(float) * TEST_W * TEST_H); + assert_non_null(scene); + assert_non_null(ref); + assert_non_null(mov); + assert_non_null(out); + _make_scene(scene, 1); + + // Per-channel gains: make the mosaic strongly colour-modulated so any + // channel cross-talk in the warp shows up as a large per-photosite error. + const float gain[3] = { 1.6f, 1.0f, 0.55f }; // R, G, B + + // Ground-truth ref->moving homography: rotation about the centre + translation. + const double th = 1.2 * M_PI / 180.0, cx = TEST_W / 2.0, cy = TEST_H / 2.0; + const double c = cos(th), s = sin(th), tx = 5.0, ty = -3.0; + const double H0 = c, H1 = -s, H2 = cx - c * cx + s * cy + tx; + const double H3 = s, H4 = c, H5 = cy - s * cx - c * cy + ty; + const double det = H0 * H4 - H1 * H3; + const double i0 = H4 / det, i1 = -H1 / det, i3 = -H3 / det, i4 = H0 / det; + const double i2 = -(i0 * H2 + i1 * H5), i5 = -(i3 * H2 + i4 * H5); + + for(int y = 0; y < TEST_H; y++) + for(int x = 0; x < TEST_W; x++) + { + const float g = gain[_fc_rggb(y, x)]; + // reference mosaic: scene sampled at (x,y), modulated by CFA colour + ref[y * TEST_W + x] = scene[y * TEST_W + x] * g; + // moving mosaic: the shifted scene, sampled at the same photosite / colour + const float sv = _bilin(scene, (float)(i0 * x + i1 * y + i2), (float)(i3 * x + i4 * y + i5)); + mov[y * TEST_W + x] = sv * g; + } + + dt_hdr_align_t *a = dt_hdr_alignment_new(NULL); + assert_non_null(a); + assert_true(dt_hdr_alignment_set_reference(a, ref, TEST_W, TEST_H, TEST_FILTERS, NULL)); + + dt_hdr_align_result_t info; + const gboolean ok = + dt_hdr_alignment_align_frame(a, mov, out, TEST_W, TEST_H, TEST_FILTERS, NULL, &info); + + if(info.status == DT_HDR_ALIGN_DISABLED) + { + assert_false(ok); + assert_memory_equal(out, mov, sizeof(float) * TEST_W * TEST_H); + } + else + { + double e_before = 0.0, e_after = 0.0; + for(int y = 48; y < TEST_H - 48; y++) + for(int x = 48; x < TEST_W - 48; x++) + { + e_before += fabs(ref[y * TEST_W + x] - mov[y * TEST_W + x]); + e_after += fabs(ref[y * TEST_W + x] - out[y * TEST_W + x]); + } + assert_true(ok); + assert_true(e_after < 0.7 * e_before); + } + + dt_hdr_alignment_free(a); + free(scene); + free(ref); + free(mov); + free(out); +} + +// The auto-reference probe must rank a feature-rich frame above a flat one. +static void test_hdr_probe_ranks_richness(void **state) +{ + (void)state; + float *flat = malloc(sizeof(float) * TEST_W * TEST_H); + float *rich = malloc(sizeof(float) * TEST_W * TEST_H); + assert_non_null(flat); + assert_non_null(rich); + for(int i = 0; i < TEST_W * TEST_H; i++) flat[i] = 0.3f; + _make_scene(rich, 1); + + const int n_flat = dt_hdr_alignment_probe_features(flat, TEST_W, TEST_H); + const int n_rich = dt_hdr_alignment_probe_features(rich, TEST_W, TEST_H); + + if(n_flat == 0 && n_rich == 0) + { + // Built without OpenCV: probing is unavailable; nothing to rank. + } + else + { + assert_true(n_rich > n_flat); + } + + free(flat); + free(rich); +} + +// One reference, several differently-shifted moving frames. The reference +// features are detected once (in set_reference) and reused by every align_frame, +// so this exercises the per-merge feature cache and confirms it is not corrupted +// across frames: every frame must still be corrected (OpenCV build) or passed +// through (no-OpenCV build). +static void test_hdr_align_reference_reused_across_frames(void **state) +{ + (void)state; + float *ref = malloc(sizeof(float) * TEST_W * TEST_H); + float *mov = malloc(sizeof(float) * TEST_W * TEST_H); + float *out = malloc(sizeof(float) * TEST_W * TEST_H); + assert_non_null(ref); + assert_non_null(mov); + assert_non_null(out); + _make_scene(ref, 1); + + dt_hdr_align_t *a = dt_hdr_alignment_new(NULL); + assert_non_null(a); + assert_true(dt_hdr_alignment_set_reference(a, ref, TEST_W, TEST_H, TEST_FILTERS, NULL)); + + const double shifts[3][2] = { { 4.0, -2.0 }, { -6.0, 3.0 }, { 2.0, 5.0 } }; + for(int t = 0; t < 3; t++) + { + const double tx = shifts[t][0], ty = shifts[t][1]; + for(int y = 0; y < TEST_H; y++) + for(int x = 0; x < TEST_W; x++) + mov[y * TEST_W + x] = _bilin(ref, (float)(x - tx), (float)(y - ty)); + + dt_hdr_align_result_t info; + const gboolean ok = + dt_hdr_alignment_align_frame(a, mov, out, TEST_W, TEST_H, TEST_FILTERS, NULL, &info); + + if(info.status == DT_HDR_ALIGN_DISABLED) + { + assert_false(ok); + assert_memory_equal(out, mov, sizeof(float) * TEST_W * TEST_H); + } + else + { + double e_before = 0.0, e_after = 0.0; + for(int y = 48; y < TEST_H - 48; y++) + for(int x = 48; x < TEST_W - 48; x++) + { + e_before += fabs(ref[y * TEST_W + x] - mov[y * TEST_W + x]); + e_after += fabs(ref[y * TEST_W + x] - out[y * TEST_W + x]); + } + assert_true(ok); + assert_true(e_after < 0.7 * e_before); + } + } + + dt_hdr_alignment_free(a); + free(ref); + free(mov); + free(out); +} + +int main(int argc, char *argv[]) +{ + (void)argc; + (void)argv; + const struct CMUnitTest tests[] = { + cmocka_unit_test(test_hdr_align_reduces_misalignment), + cmocka_unit_test(test_hdr_align_cfa_modulated), + cmocka_unit_test(test_hdr_align_reference_reused_across_frames), + cmocka_unit_test(test_hdr_probe_ranks_richness), + }; + return cmocka_run_group_tests(tests, NULL, NULL); +} +// clang-format off +// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py +// vim: shiftwidth=2 expandtab tabstop=2 cindent +// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified; +// clang-format on diff --git a/src/tests/unittests/test_hdr_alignment_internal.c b/src/tests/unittests/test_hdr_alignment_internal.c new file mode 100644 index 000000000000..c83562048d42 --- /dev/null +++ b/src/tests/unittests/test_hdr_alignment_internal.c @@ -0,0 +1,196 @@ +/* + This file is part of darktable, + Copyright (C) 2026 darktable developers. + + darktable 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. + + darktable 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 darktable. If not, see . +*/ + +/* White-box unit tests for the internal (static) helpers of the HDR alignment + * code. Following the test_filmicrgb.c pattern, the translation unit is + * #included directly so the static helpers are reachable; lib_darktable supplies + * everything hdr_alignment.c links against (the OpenCV seam, glib, ...). + * + * Complements test_hdr_alignment.c, which exercises the public API end-to-end. */ + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "common/hdr_alignment.c" + +#ifdef _WIN32 +#include "win/main_wrapper.h" +#endif + +static int _cmp_f(const void *a, const void *b) +{ + const float x = *(const float *)a, y = *(const float *)b; + return (x > y) - (x < y); +} + +// The histogram-based percentile must track the exact (sorted) order statistics +// to well under one part in a hundred of the data span. +static void test_percentile_bounds_accuracy(void **state) +{ + (void)state; + const size_t n = 300000; + float *v = malloc(n * sizeof(float)); + float *c = malloc(n * sizeof(float)); + assert_non_null(v); + assert_non_null(c); + unsigned s = 4242u; + for(size_t i = 0; i < n; i++) + { + s = s * 1103515245u + 12345u; + v[i] = ((s >> 9) & 0x7fff) / 32767.0f * 8.0f - 1.0f; // ~[-1, 7] + } + memcpy(c, v, n * sizeof(float)); + qsort(c, n, sizeof(float), _cmp_f); + const float ex_lo = c[(size_t)(0.01 * (n - 1) + 0.5)]; + const float ex_hi = c[(size_t)(0.99 * (n - 1) + 0.5)]; + + float p_lo, p_hi; + _percentile_bounds(v, n, DT_HDR_PERCENTILE_LOW, DT_HDR_PERCENTILE_HIGH, &p_lo, &p_hi); + + const float span = ex_hi - ex_lo; + assert_true(fabsf(p_lo - ex_lo) < 0.01f * span); + assert_true(fabsf(p_hi - ex_hi) < 0.01f * span); + free(v); + free(c); +} + +// A flat image has no spread: the bounds must collapse to the constant value so +// the downstream stretch is a well-defined no-op (no division blow-up). +static void test_percentile_bounds_flat(void **state) +{ + (void)state; + float v[1024]; + for(int i = 0; i < 1024; i++) v[i] = 0.5f; + float p_lo = -1.0f, p_hi = -1.0f; + _percentile_bounds(v, 1024, DT_HDR_PERCENTILE_LOW, DT_HDR_PERCENTILE_HIGH, &p_lo, &p_hi); + assert_true(p_lo == 0.5f); + assert_true(p_hi == 0.5f); +} + +// _proxy_to_u8 must be monotone (preserves ordering) and use the full code range +// on a ramp, so SIFT sees the intended contrast. +static void test_proxy_to_u8_monotone(void **state) +{ + (void)state; + const int pw = 256, ph = 1; + float ramp[256]; + for(int i = 0; i < pw; i++) ramp[i] = (float)i; + uint8_t *u8 = _proxy_to_u8(ramp, pw, ph, 2.2); + assert_non_null(u8); + for(int i = 1; i < pw; i++) assert_true(u8[i] >= u8[i - 1]); + assert_int_equal(u8[0], 0); + assert_int_equal(u8[pw - 1], 255); + dt_free_align(u8); +} + +#ifdef HAVE_OPENCV +// The Bayer same-color fast path must be bit-for-bit identical to the general +// X-Trans-capable sampler it replaces, across the whole sampling domain +// (including out-of-bounds edges) and every CFA colour. A divergence here would +// silently corrupt every aligned Bayer frame. +static void test_bayer_sampler_bit_exact(void **state) +{ + (void)state; + const int W = 71, H = 53; // odd dims exercise the edges + const uint32_t filters = 0x94949494u; // RGGB + float *m = malloc(sizeof(float) * W * H); + assert_non_null(m); + for(int y = 0; y < H; y++) + for(int x = 0; x < W; x++) + m[y * W + x] = (float)((x * 7 + y * 13) % 101) / 101.0f + + 0.013f * (float)((x * x + 3 * y) % 17); + + for(double sy = -2.0; sy <= H + 2.0; sy += 0.31) + for(double sx = -2.0; sx <= W + 2.0; sx += 0.27) + for(int c = 0; c < 3; c++) // R=0, G=1, B=2 + { + const float fast = _sample_bayer_same_color(m, W, H, filters, c, sx, sy); + const float general = _sample_cfa_same_color(m, W, H, filters, NULL, c, sx, sy); + assert_true(fast == general); + } + free(m); +} + +// The general (X-Trans-capable) same-color sampler must never mix colours: with +// a mosaic where only one colour's photosites are non-zero, sampling that colour +// anywhere in-bounds (including the borders that exercise the outward-search +// fallback) returns that value, and sampling a different colour returns zero. +// This is the invariant that keeps a warped mosaic a valid CFA in the reference +// frame's phase. +static void test_xtrans_sampler_same_color_only(void **state) +{ + (void)state; + const int W = 48, H = 42; + // A valid X-Trans 6x6 layout (8 R, 20 G, 8 B per tile). + static const uint8_t xtrans[6][6] = { + { 1, 1, 0, 1, 1, 2 }, + { 0, 2, 1, 2, 0, 1 }, + { 1, 1, 2, 1, 1, 0 }, + { 1, 1, 0, 1, 1, 2 }, + { 2, 0, 1, 0, 2, 1 }, + { 1, 1, 2, 1, 1, 0 }, + }; + float *m = malloc(sizeof(float) * W * H); + assert_non_null(m); + + // Only green (colour 1) photosites carry signal; R and B are zero. + const int target = 1; + for(int y = 0; y < H; y++) + for(int x = 0; x < W; x++) + m[y * W + x] = (_fcol(y, x, 9u, xtrans) == target) ? 1.0f : 0.0f; + + for(double sy = -1.0; sy <= H + 1.0; sy += 0.5) + for(double sx = -1.0; sx <= W + 1.0; sx += 0.5) + { + // sampling green sees only green (all 1.0) -> exactly 1.0 + const float g = _sample_cfa_same_color(m, W, H, 9u, xtrans, target, sx, sy); + assert_true(g > 0.999f && g < 1.001f); + // sampling red sees only red photosites, all of which are 0.0 + const float r = _sample_cfa_same_color(m, W, H, 9u, xtrans, 0, sx, sy); + assert_true(r >= 0.0f && r < 0.001f); + } + free(m); +} +#endif // HAVE_OPENCV + +int main(void) +{ + const struct CMUnitTest tests[] = { + cmocka_unit_test(test_percentile_bounds_accuracy), + cmocka_unit_test(test_percentile_bounds_flat), + cmocka_unit_test(test_proxy_to_u8_monotone), +#ifdef HAVE_OPENCV + cmocka_unit_test(test_bayer_sampler_bit_exact), + cmocka_unit_test(test_xtrans_sampler_same_color_only), +#endif + }; + return cmocka_run_group_tests(tests, NULL, NULL); +} + +// clang-format off +// modelines: These editor modelines have been set for all relevant files by tools/update_modelines.py +// vim: shiftwidth=2 expandtab tabstop=2 cindent +// kate: tab-indents: off; indent-width 2; replace-tabs on; indent-mode cstyle; remove-trailing-spaces modified; +// clang-format on