Hi,

I need to add photos taken with different cameras and have the sum be
physically meaningful. Every image must therefore land in the same linear XYZ space, with a
saturation mask, a clipped photosite records a lower bound, not a measurement, so those
pixels have to be excluded.

I decode with LibRaw in what I hope is a strictly linear mode, gamma=(1,1),
no_auto_bright=True, user_wb=[1,1,1,1],
output_color=raw, half_size=True,
adjust_maximum_thr=0, then apply the DNG colour matrix myself, interpolated
between ColorMatrix1/2 at the as-shot white point.
adjust_maximum_thr=0 is deliberate: letting maximum follow each
image's own content would put every file on a different scale, which is fatal when the
whole point is to add two files together.

My mask was bayer >= black + 0.99 * (white_level - black) and it came out
empty on most images, which seemed wrong for 1500x2000 photos with obviously blown skies.
So I measured white_level against the actual plateau in the raw data (the code
where tens of thousands of photosites pile up), over 120 FiveK DNGs, all Adobe DNG
Converter output:

camera                white_level   plateau   clipped files
Canon EOS 5D                 3692      3692        10/20
Nikon D700                  15892     16383         9/11
Canon EOS 40D               13600     13824          2/3
Nikon D70s                   4095      3120         4/13
Canon EOS 20D                4095      3594          3/9

1. What is the right saturation threshold? The data goes above
   white_level on some bodies and plateaus well below it on others. My
   reading is that color.maximum is a linearity limit rather than a clipping
   level, that values above it are legal in a DNG and get flattened by CLIP() in
   scale_colors, and that on bodies where it is merely the ADC range (4095) the
   analog saturation is lower and nothing reports it. Is
   min(white_level, measured_plateau) the correct threshold, and do I have to
   detect the plateau from the histogram myself and build a per-model table?
   (color.linear_max equals maximum on all 23 models in my sample,
   so it gives me nothing.)
2. What is CameraWhite for? The DNG SDK computes
   CameraWhite = ColorMatrix * XYtoXYZ(whiteXY), normalised so its max entry is
   1, about (0.45, 1.0, 0.55) on my files. I used it to clamp the decoded image with
   cam = min(cam, CameraWhite), which clamps red to 0.45 everywhere,
   destroying saturated reds that were never clipped, and breaks additivity since
   min(t,w) + min(r,w) != min(t+r,w). Is it right that it belongs only to
   highlight reconstruction, telling me a clipped channel's true value was above the
   ceiling, and has no place in a linear decode?
3. Being honest with you, I'm a bit lost about what to do and I do not find good ressources online, AI is not very helpful too. I Would be grateful for any enlightment!!

Here's my code for more details:

I follow a pythonized version of the adobe DNG SNK to get the color_matrix and the whiteXY point.

def read_raw_to_xyz(path, half_size=True):
    with rawpy.imread(str(path)) as raw:
        bayer = raw.raw_image_visible
        sat_b = bayer >= black + 0.99 * (raw.white_level - black)          # <- my question is here
        H, W  = (sat_b.shape[0]//2)*2, (sat_b.shape[1]//2)*2
        saturated = sat_b[:H, :W].reshape(H//2, 2, W//2, 2).any(axis=(1, 3))
 
        cam = raw.postprocess(
            gamma=(1, 1), no_auto_bright=True,    # no tone curve, no auto exposure
            output_bps=16,
            adjust_maximum_thr=0,                 # keep `maximum` independent of image content
            use_camera_wb=False, use_auto_wb=False,
            user_wb=[1.0, 1.0, 1.0, 1.0],         # no WB: I do it myself in XYZ
            output_color=rawpy.ColorSpace.raw,    # no colour conversion: I apply the DNG matrix
            half_size=half_size,                  # avoids the non-linear demosaic algorithms
        )
 
    cam_to_xyz, xy = camera_to_xyz(path)          # inv(AB · CC · CM(xy))
    cam = cam.astype(np.float32) / 65535.0
    xyz = np.einsum("ij,hwj->hwi", cam_to_xyz, cam)
    return xyz, xy_to_xyz(xy), saturated # image in the xyz space, illuminance, saturation mask

Thanks.