Hello all,
I have been trying to process raw images in python (with libraw) such that the output looks the same as if I applied white/black point, daylight white balance, demosaic, and output to linear rec. 2020 in darktable (or even XYZ).
The goal is to do some processing on the raw image (possibly even demosaic) and then import it in darktable (or any image processor), select the right input color profile (eg XYZ or lin. rec. 2020), and develop it exactly as if it were raw (post-demosaic).
Unfortunately the result is always far from what it should be. I’ve tried various combinations and tricks and at best I get an off white balance. Most of the time I get something very red/purple/wrong.
My basic process is to apply the black/white point, multiply raw values by the normalized (s.t. green = x1) “daylight” white balance, debayer, then do profiledRGB_image = debayered_image @ invert(xyz_to_cam) @ xyz_to_profiledRGB.
What is obviously wrong with my process?
Below is the sample source code which will reproduce this and download a sample image:
import os
import requests
import numpy as np
import rawpy
import cv2
SAMPLE_RAW_URL = (
"https://nc.trougnouf.com/index.php/s/zMA8QfgPoNoByex/download/DSC01568.ARW"
)
SAMPLE_RAW_FPATH = os.path.join("data", SAMPLE_RAW_URL.rpartition("/")[-1])
def raw_fpath_to_mono_img_and_metadata(fpath: str) -> tuple[np.ndarray, dict]:
"""Convert from raw fpath to (tuple) [1,h,w] numpy array and metadata dictionary."""
# step 1: get raw data and metadata
rawpy_img = rawpy.imread(fpath)
metadata = dict()
metadata["black_level_per_channel"] = rawpy_img.black_level_per_channel
metadata["white_level"] = rawpy_img.white_level # imgdata.rawdata.color.maximum
metadata["daylight_whitebalance"] = rawpy_img.daylight_whitebalance
metadata["rgb_xyz_matrix"] = rawpy_img.rgb_xyz_matrix
metadata["daylight_whitebalance_norm"] = np.array(
metadata["daylight_whitebalance"], dtype=np.float32
)
metadata["daylight_whitebalance_norm"][3] = metadata[
"daylight_whitebalance_norm"
][1]
metadata["daylight_whitebalance_norm"] /= metadata[
"daylight_whitebalance_norm"
][1]
mono_img = rawpy_img.raw_image
mono_img = np.expand_dims(mono_img, axis=0)
return mono_img, metadata
def scale_to_bw_points(img: np.ndarray, metadata: dict) -> None:
"""
Scale image to black/white points described in metadata.
compat: bool: use the conservative (higher) "white_level" value as in darktable.
"""
scaled_img = img.astype(np.float32)
for ch in range(img.shape[-3]):
scaled_img[ch] -= metadata["black_level_per_channel"][ch]
# step 4: normalize s.t. white level is 1: divide each value by (white level - black level)
# darktable uses global max only, aka "white_level", which is higher.
vrange = metadata["white_level"] - metadata["black_level_per_channel"][ch]
scaled_img[ch] /= vrange
return scaled_img
def apply_whitebalance(
img: np.ndarray, metadata: dict, in_place: bool = True
) -> None:
"""
Apply camera white balance.
"""
# step 5: apply camera reference white balance
if not in_place:
img = img.copy()
img[0, 0::2, 0::2] *= metadata["daylight_whitebalance_norm"][0] # R
img[0, 0::2, 1::2] *= metadata["daylight_whitebalance_norm"][1] # G
img[0, 1::2, 1::2] *= metadata["daylight_whitebalance_norm"][2] # B
img[0, 1::2, 0::2] *= metadata["daylight_whitebalance_norm"][3] # G
assert img.max() <= 1.0, img.max()
return img
def demosaic(mono_img: np.ndarray):
"""
Transform mono image to camRGB colors.
"""
# convert to uint16 and scale to ensure we don't lose negative / large values
black_offset = 0 if mono_img.min() >= 0 else -mono_img.min()
mono_img += black_offset
max_value = 1 if mono_img.max() <= 1 else mono_img.max()
mono_img /= max_value
# mono_img = mono_img.clip(0, 1)
assert mono_img.min() >= 0 and mono_img.max() <= 1.0, (
f"demosaic: image is out of bound; destructive operation. {mono_img.min()=}, "
f"{mono_img.max()=}"
)
mono_img *= 65535
mono_img = mono_img.astype(np.uint16).reshape(mono_img.shape[1:] + (1,))
rgb_img = cv2.demosaicing(mono_img, cv2.COLOR_BayerRGGB2RGB)
rgb_img = rgb_img.transpose(2, 0, 1) # opencv h, w, ch to numpy ch, h, w
rgb_img = rgb_img.astype(np.float32) / 65535.0 * max_value - black_offset
return rgb_img
def get_XYZ_to_profiledRGB_matrix(profile: str) -> np.ndarray:
if profile == "lin_rec2020":
return np.array(
[
[1.71666343, -0.35567332, -0.25336809],
[-0.66667384, 1.61645574, 0.0157683],
[0.01764248, -0.04277698, 0.94224328],
],
dtype=np.float32,
)
else:
raise NotImplementedError(
f"get_std_profile_matrix: {profile} not lin_rec2020."
)
def get_camRGB_to_profiledRGB_matrix(
metadata: dict, output_color_profile: str
) -> np.ndarray:
"""Get conversion matrix from camRGB to a given color profile."""
cam_to_xyzd65 = np.linalg.inv(metadata["rgb_xyz_matrix"][:3])
if output_color_profile.lower() == 'xyz':
return cam_to_xyzd65
xyz_to_profiledRGB = get_XYZ_to_profiledRGB_matrix(output_color_profile)
color_matrix = cam_to_xyzd65 @ xyz_to_profiledRGB
return color_matrix
def camRGB_to_profiledRGB(
camRGB_img: np.ndarray, metadata: dict, output_color_profile: str
) -> np.ndarray:
"""Convert camRGB debayered image to a given RGB color profile."""
color_matrix = get_camRGB_to_profiledRGB_matrix(metadata, output_color_profile)
orig_dims = camRGB_img.shape
#breakpoint()
profiledRGB_img = (camRGB_img.swapaxes(0,2) @ color_matrix).swapaxes(0,2)
#profiledRGB_img = (color_matrix @ camRGB_img.reshape(3, -1)).reshape(orig_dims)
return profiledRGB_img
def rgb_img_to_file(img: np.ndarray, fpath: str) -> None:
"""Save (c,h,w) RGB image to file."""
hwc_img = img.transpose(1, 2, 0) # numpy ch, h, w to opencv h, w, ch
hwc_img = cv2.cvtColor(hwc_img, cv2.COLOR_RGB2BGR)
hwc_img = (hwc_img * 65535).clip(0,65535).astype(np.uint16) # while debugging
cv2.imwrite(fpath, hwc_img)#.clip(0, 1))
def get_sample_raw_file(url: str = SAMPLE_RAW_URL) -> str:
"""Get a testing image online."""
fn = url.split("/")[-1]
fpath = os.path.join("data", fn)
if not os.path.exists(fpath):
os.makedirs("data", exist_ok=True)
r = requests.get(url, allow_redirects=True, verify=False)
open(fpath, "wb").write(r.content)
return fpath
if __name__ == "__main__":
raw_fpath = get_sample_raw_file(url=SAMPLE_RAW_URL)
out_base_path = os.path.join("tests", f"raw.py.main")
mono_img, metadata = raw_fpath_to_mono_img_and_metadata(raw_fpath)
print(f"raw.py: opened {raw_fpath} with {metadata=}")
# prepare (common)
mono_img = scale_to_bw_points(mono_img, metadata)
mono_img_wb = apply_whitebalance(mono_img, metadata, in_place=False)
camRGB_img = demosaic(mono_img_wb)
# output representations
xyz_img = camRGB_to_profiledRGB(camRGB_img, metadata, 'XYZ')
lin_rec2020_img = camRGB_to_profiledRGB(camRGB_img, metadata, "lin_rec2020")
os.makedirs("tests", exist_ok=True)
rgb_img_to_file(lin_rec2020_img, out_base_path+'.lin_rec2020.tiff')
rgb_img_to_file(
img=camRGB_img,
fpath=out_base_path + ".camRGB.tiff",
)
rgb_img_to_file(
img=xyz_img,
fpath=out_base_path + ".XYZ.tiff",
)
print(f"raw.py: output images saved to {out_base_path}.*")
but it seems to be all there
