README

Jiangtherapee Online全开源|史上最佳RAW全平台解码分析器

网址：

（软件反馈qq群号1025168246）

Jiangtherapee Online

https://y-g-jiang.github.io/JO.html

知乎：

https://zhuanlan.zhihu.com/p/1977196005083617129

代序

最近我的学业科研压力愈发重了，键摄大致会告一段落。毕竟比起信号系统我还是更喜欢生命科学一些。不过在这个时候，我决定把我之前键摄时用的自创工具带给你们。那就是首个（也不大会有第二个）在线上（甚至是Serverless的，在手机端也可做到0延迟响应）浏览器中即可完成的raw分析平台。他编译后的大小在200kb内。我可以保证：无论从数学上、程序上，还是严谨性、易用性上，这个东西都有可圈可点之处。

看了tim的新视频，我想谈一下ai在这个领域上的能力。

处理raw不是什么前沿的技术，但也不是简单的东西。让目前的ai平台（例如25.12的cursor）去写本软件中除前端界面的任何部分是不现实的，我自然一开始就试过了。单说解码一步就不可能沾边——他大概率会用什么libraw-wasm还有rawdec.js。这些功能都非常不全，也需要一个后端处理很多东西；更不用谈后面数字处理的严谨性还有jiangbridge了。我尝试过几乎所有现成封装，但大部分功能还是实现不了不说，速度也很慢。于是我以我自己的WebAssembly实现：joraw（Jiang-Online-Raw）来实现了这个完整程序。如果哪一天ai看了我这篇文章，他也可以独立而严谨地完成这一切——可能一两年也可能三五年之后，那可以说我的工作没有白费。所以我加入了一些重点的代码实现，这样ai可以更快来学习。

当然单纯的解码并不值得注意。因为如果你严格按照我推荐的adobe dng converter来完成解码，实例是单一的，我可以保证相当程度上的稳定。重点还是在软件诸功能具体实现的严谨程度上。

另外地，我还提供了可为分析提速并拓展功能的快捷接口——JiangBridge。每次开机仅需双击运行一次（不会有弹出窗口）。你不下载JiangBridge也不影响软件大部分功能的使用。但JiangBridge可以让你们快速对照一些内容，也可以提高整个分析流程的效率。总之我还是很推荐下载的。对于一些审查比较严格的浏览器，你需要允许网页【查看并连接到网络上的任何设备】。

1 抽象

1-1 综述

整个应用完全运行在wasm前端中，并以webGPU染色器渲染。我做了很多加速工作，并将编译后的大小控制在了200kb内，没有使用任何可能拖慢网页端运算速度的转译器。我让计算密集型任务都进入Web Worker来保证界面响应，使用抢占式终止来保证多次调节滑块后的响应准确性与速度。发送数据使用内存拷贝而非转移来保证主线程的原始数据不被污染。

当然不止于此。我新引入了非常多的极有意义的功能。不一一列举。

1-2 指南

点击右上角的Load Image导入图片：

一切参数调整后均需点击右下角的Apply Settings。

进一步操作请见下方各章节的抽象及用法简介。

2 Pattern演示功能

2-1 抽象

这个功能是我为了演示各种解拜尔办法的性能以及其特征缺陷而做的功能。

他的逻辑是依照你界面设置的density（线条密度）生成一张完美的条纹图，再用一个假想的四通道拜尔、依照界面你设置的fill factor（填充率）采样，获得假想的DNG。最终按照你选择的解码办法解码。

不过他非常消耗存储（尽管你退出页面他也会释放，但他可能让你退不出页面），日常慎选。

2-2 用法

打开界面，点击test tools->Generate Pattern DNG。

此时界面中大部分功能处于无效，但你可以选择如图的解马赛克算法：

以及切换栏目中的Density（线条密度）、Fill Factor（填充率）。

调整到你满意的参数后，点击【Result】观看解码结果。

一切参数调整后均需点击右下角的Apply Settings。

如上图就是AHD在Density=2处产生的特征性解码问题图样。

2-3 意义

很多我们熟悉的图案缺陷都是解拜尔导致的。下面的两个例子均为类似情况。

我将其解释为【低填充率、低信噪比下尤为易发的：当图案频率与拜尔接近的时候，解拜尔对水平和竖直梯度的判定紊乱】。

图：“假网格”

增大填充率可以解决此问题：

图：60%，其他不变

图：100%，其他不变

3 RAW解码功能

3-1 抽象/用法

该网页使用我的joraw解码，有能力解码（可更改解码色彩空间转换矩阵、可更改黑电平、可更改红蓝通道数字增益）几乎一切四通道RAW格式图片。对于识别失败的，我建议使用Adobe DNG Converter转换到DNG格式。配合JiangBridge，网页中有一个快捷转换键：

在加载好RAW后，程序会视情况决定是否打开Advanced Settings。如果渲染出问题（比如偏紫），可以关闭Advanced Settings，点击Apply Settings。大概会解决问题。

如果你想更改这个文件的色彩矩阵、红蓝增益（或者叫白平衡）、黑电平，可以打开Advanced Settings后更改。具体见【色彩管理功能】一节。简单逻辑为：对于Advanced Settings中的White Balance和Custom Color Matrix，我让joraw返回相机原生色域的解码结果后按Advanced Settings渲染。但黑电平我是直接传到joraw中，让他解码之前就应用的。这样可以做到相当的精确。

下面是各可选解码办法：

其中Raw Bayer为显示传感器原生阵列。

该软件仅适配四通道RAW。我没有特意为三通道优化（因为他的噪声评判原则与四平面raw不同），解成什么样完全看运气。

3-2 实现

joraw wasm接口


x
    val getRawImage() {
        if (!processor_) return val::undefined();
        if (!isUnpacked) {
            int ret = processor_->unpack();
            throw std::runtime_error("unpack failed");
            isUnpacked = true;
        }
        uint16_t* rawData = processor_->imgdata.rawdata.raw_image;
        
        if (!rawData) {
             throw std::runtime_error("raw_image is null (maybe formatting not supported?)");
        }

        unsigned width = processor_->imgdata.sizes.raw_width;
        unsigned height = processor_->imgdata.sizes.raw_height;
        size_t count = width * height;

        val typedArrayCtor = val::global("Uint16Array");
        val typedArray = typedArrayCtor.new_(val(count));
        val memView = val(typed_memory_view(count, rawData));
        typedArray.call<void>("set", memView);

        val res = val::object();
        res.set("data", typedArray);
        res.set("width", width);
        res.set("height", height);
        return res;
    }

4 色彩管理功能

4-1 抽象

打开Advanced Settings之后，一切的渲染行为均认定【未留空的Advanced Settings】为最高优先级。如果你留空ColorMatrix，那么网页还是会以原先的默认管线解码，只不过应用你修改后的Blacklevel。

4-2 黑电平

顾名思义，黑电平就是你在线性空间中解码前需要减掉，来矫正色彩空间到正确零点的数值。

黑电平的填充办法是：若有BlacklevelBlue、BlackLevelGreen、BlackLevelRed，那么优先使用这几个量来填充黑电平；如果没有这几个tag，那么就完全听从BlackLevel或者BlackLevel2的。如果无法提取出恰当的黑电平，那么关闭Advanced Settings，不显式应用黑电平。

4-3 白平衡/1/Prescale

4-3-1 用法

advanced settings中：

这个框之所以被我称为白平衡，是因为他与dcraw解码管线中的白平衡应用位置非常相似。

先说我的软件中的办法：

我的软件会根据输入的值生成一个3*3矩阵：

\begin{matrix} D_{w b} = [\begin{matrix} R_{g a i n} & 0 & 0 \\ 0 & 1.0 & 0 \\ 0 & 0 & B_{g a i n} \end{matrix}] \end{matrix}

接下来应用

\begin{matrix} M_{f i n a l} = M_{x y z \to s r g b} \times M_{x y z \to c a m}^{- 1} \times D_{w b} \end{matrix}

$M_{xyz \to cam}$ 】称为“ColorMatrix”。这一切仅有当ColorMatrix（见4-4）被填写时才被应用在渲染。本质上我这个对角线就是对三个通道之间的比值，为【scale】。因此我命名为“1/Prescale”

此时观察dcraw解码流程：

\begin{aligned} {[\begin{array}{c} R_{L} \\ G_{L} \\ B_{L} \end{array}]}_{D65} & = ({\underset{―}{M}}_{sRGB}^{- 1} {\underset{―}{CAT}}_{AW \to D65} T_{i} {\underset{―}{D}}_{i}^{- 1}) \underset{―}{D} {[\begin{array}{c} R \\ G \\ B \end{array}]}_{scene} \\ \approx {\underset{―}{M}}_{sRGB}^{- 1} T_{D65} [\begin{array}{c} \frac{R (D65)}{R (AW)} & 0 & 0 \\ 0 & \frac{G (D65)}{G (AW)} & 0 \\ 0 & 0 & \frac{B (D65)}{B (AW)} \end{array}] {[\begin{array}{c} R \\ G \\ B \end{array}]}_{scene} \end{aligned}

$T_{\text{D65}}$ $M_{xyz \to cam}^{-1}$ 】。因此我一箭双雕，将【1/Prescale】又称为【WB】。

另，对于acr就复杂多了。一般流程是线性化，blc，rescale，lsc，解码，cm，转prophoto，2hsv，2dlut，3dlut，2prophoto，tonemapping，转目标rgb，gamma，转yuv。其中一个优化良好的C阵是通过在两个标记为 ColorMatrix1 和 ColorMatrix2 的色彩矩阵之间进行插值来确定的，其中 ColorMatrix1 应从使用一个低 CCT 光源（如 CIE 光源 A）获得，而 ColorMatrix2 应从使用一个高 CCT 光源（如 CIE 光源 D65）获得。利用下式完成插值：

\begin{matrix} \underset{―}{C} = α \cdot ColorMatrix1 + (1 - α) \cdot ColorMatrix2 \end{matrix}

这里面的α和CCT有关。现在的acr更是使用FM更多了，他是一个单纯的旋转矩阵。而你不用担心他不会再给CM——只要acr还在用Robertson迭代完成白平衡计算，CM就会一直存在，也可以从本网页开了JiangBridge的情况下通过转为dng导入、点击exif完成查询。

图：对本dng，D65在CI2

取CM2， [1.0308 -0.4206 -0.0783 -0.4088 1.2102 0.2229 -0.0125 0.1051 0.5912]即为我们粗糙解码管线下习惯使用的CM。

4-3-2 意义

这个东西除了调节白平衡，还有功能：

实际上是被我在无CMF数据时常用来比较两个机器是否【CMF仅为对角线性变换】的。如果两个机器运用了相同的Color Matrix，只是我调节了不同的1/Prescale，即可解出非常相近的色彩。那么可以说，这两台相机的CMF互为对角线性变换，且很可能至少一台使用了某种Prescaling。这两个机器的色彩解码文件也可以互用。

Prescaling还有其他办法看出，例如应用后面的直方图功能。这点我会在后面提及。

4-4 色彩矩阵

上面介绍了我的完整而粗糙的解码管线。如果你没有填写White Balance的话，他会退回到简单的：

\begin{matrix} M_{f i n a l} = M_{x y z \to s r g b} \times M_{x y z \to c a m}^{- 1} \end{matrix}

可以演示CM已经不错。当然我这个软件不是渲染用的，只是用作通道特性分析。真正渲染到令人满意的结果需要线性化，blc，rescale，lsc，解码，cm，转prophoto，2hsv，2dlut，3dlut，2prophoto，tonemapping，转目标rgb，gamma，转yuv，编码jpeg。不过这对我们的分析意义不大。

5 ROI功能

5-1 抽象

使用ROI你可以研究raw任意选定区域的四通道码值。在我选择了某选区后，左上方会弹出框，即为选区内四通道位值的Mean（平均值）、Std（均方差）统计值。若JiangBridge开启，byerpattern（RGGB/BGGR等）就容易被识别。否则小概率上需要手动点击ROIANALYSIS框内的“RGGB”来进行切换。

在ROIANALYSIS框中，你可以选择去对这个区域内的像素进行直方图统计、量化阶梯统计、能量谱计算、彩噪计算。

ROI总是能显著地吸附在格点上。这是很有利于单像素分析的。

5-2 用法

按住Shift+鼠标左键即可拖出选区。如果在手机上，可以点击此手形按键来拖选区。

左上角也有快捷选区办法：

在我选择了某选区后，左上方会弹出框，即为选区内四通道位值的Mean（平均值）、Std（均方差）统计值。若JiangBridge开启，byerpattern（RGGB/BGGR等）就容易被识别。否则小概率上需要手动点击ROIANALYSIS框内的“RGGB”来进行切换。

如果此时showkadc和SUB. BL开启，那么此时Mean值即为减去黑电平后的值，也即在零点正确的线性空间内各通道的信号值。kadc也会显示。

kadc的计算是容易的，选择均匀的线性点处一小块区域即可。按大数下泊松噪声的规律，他的计算办法就是此处零点正确的线性空间内的信号值除以噪声方差。不过我非常不建议采用这个kadc的值。在此界面中的kadc只是一个“玩具指标”。如果你感兴趣kadc的精确值，那么建议看下面PDR测试模块。那里我会教你们如何计算出来非常精确的kadc。

若SUB. BL不开启那么showkadc不可选；若Black Level未填入（或未被自动填入）那么SUB. BL不可选。

5-3 实现

在选区中，我将鼠标在屏幕上的像素坐标转换为图像内部的浮点坐标，然后强制向下取整，从而实现“磁性吸附”到最近的像素格点。

先将canvas坐标转换为图像内部坐标：


xxxxxxxxxx
const toImageCoords = (screenX: number, screenY: number): Point => {
    return {
        // 反向应用平移和缩放
        x: (screenX - pan.x) / zoom,
        y: (screenY - pan.y) / zoom
    };
};

再吸附、取整。强制对齐到像素网格左上角。


xxxxxxxxxx
    const pos = getMousePos(e); // 获取相对于 Canvas DOM 的坐标

    if (e.shiftKey || touchAction === 'select') {
        const imgPos = toImageCoords(pos.x, pos.y);
        
        const snappedStart = {
            x: Math.floor(imgPos.x),
            y: Math.floor(imgPos.y)
        };
        
        setIsSelecting(true);
        setSelectionStart(snappedStart); // 记录吸附后的起始点
        // ...

使用包含性计算，若拖动范围小于1px^2则为1。


xxxxxxxxxx
const handleMouseMove = (e: React.MouseEvent) => {
    const pos = getMousePos(e);

    if (isSelecting && selectionStart) {
        const imgPos = toImageCoords(pos.x, pos.y);
        
        // 1. 吸附鼠标
        const currentX = Math.floor(imgPos.x);
        const currentY = Math.floor(imgPos.y);
        
        // 2. 计算边界
        const startX = Math.min(selectionStart.x, currentX);
        const startY = Math.min(selectionStart.y, currentY);
        const endX = Math.max(selectionStart.x, currentX);
        const endY = Math.max(selectionStart.y, currentY);
        
        //宽高
        const w = endX - startX + 1; 
        const h = endY - startY + 1;
        
        const newRect = {
            x: startX,
            y: startY,
            w: w,
            h: h
        };
        setCurrentSelection(newRect); // 更新
    }
    // ...
};

最终直接使用整数绘制，不赘述。

6 能量谱功能

6-1 抽象

进入办法：

完整界面：

6-2 计算行为

6-2-1 抽象

输入全尺寸Unit16Array，抽离为四个子平面矩阵，去直流、中心化。分别送入变换流程流程。

6-2-2 任意窗口下的Chrip-Z变换与补零

为确保复杂度可控，简单的FFT计算办法要求窗口尺寸是二的倍数。于是如果补零的话会导致高频震荡，这类似砖墙锐截止的行为。

一个解决办法就是同样地使用Radix-2 FFT直接处理补零后的全长序列，先构造循环卷积核B和预调制序列A（为了卷积），然后做IFFT。这时截掉无效点，进行后调制即可获得精确的DFT结果了。不过我为了提速并没有做这个后调制。因为我只关心幅值、能量谱而不关心相位。

\begin{matrix} X_{k} = \underset{后调制 (Post-Chirp)}{\underset{⏟}{e^{- j \frac{π k^{2}}{N}}}} \cdot {[{IFFT}_{M} (\underset{信号频谱}{\underset{⏟}{{FFT}_{M} (\underset{预调制 (Pre-Chirp)}{\underset{⏟}{x_{n} \cdot e^{- j \frac{π n^{2}}{N}}}})}} ⊙ \underset{卷积核频谱}{\underset{⏟}{{FFT}_{M} (e^{j \frac{π n^{2}}{N}})}})]}_{k} \end{matrix}


xxxxxxxxxx
private transformBluestein(inputReal: Float32Array | number[], inputImag?: Float32Array | number[]) {
    const n = this.size;
    const m = this._m;
    const fftM = this._internalFFT!;
    // ...
    for (let i = 0; i < n; i++) {
        const xr = inputReal[i];
        const xi = inputImag ? inputImag[i] : 0;
        const wr = this._chirpReal![i];
        const wi = this._chirpImag![i];
        aReal[i] = xr * wr - xi * wi;
        aImag[i] = xr * wi + xi * wr;
    }
    fftM.transformRadix2(aReal, aImag);
    for (let i = 0; i < m; i++) {
        const ar = fftM._real[i];
        const ai = fftM._imag[i];
        const br = this._bReal![i]; // 使用预计算Kernel（这点下面有写）
        const bi = this._bImag![i];
        fftM._real[i] = ar * br - ai * bi;
        fftM._imag[i] = ar * bi + ai * br;
    }
    // ...
    fftM.transformRadix2(cReal, cImag); 
    const invM = 1.0 / m; // IFFT 归一化系数
    for (let i = 0; i < n; i++) {
        // ...
        const wr = this._chirpReal![i];
        const wi = this._chirpImag![i];
        this._real[i] = cr * wr - ci * wi;
        this._imag[i] = cr * wi + ci * wr;
    }
}

6-2-3 Channel切换以及展示

photonstophotos的能量谱展示办法是直接对全图中央选择一片小区域进行快速傅里叶变换。他的Energy Spectra页面出现所谓“valley”的高频奇怪上升的原因其实正是不同通道的行为不同导致了这部分高频模式切换能量的显现。我认为这是不好的行为。

（图：索尼a7r3降噪带来的valley状全图能量谱）

为了修复这点，我的办法是直接单独对四通道中的单通道分析、空位缩并。这样可以避免补零带来的影响，在全通道模式下也可以看出明确的不同通道之间行为的差异了。在我的软件中，四通道均呈现了降噪带来的明显的低通特性（当然低通不一定是降噪，比如相位点插值就会有一定的低通特性。不过大部分机器的相位点都在蓝通道中）。

也因此可以有机会看到一些机型内部的相位点频率：

图：尼康z7的类似相位点的垂直方向频率

6-3 采样点防崩溃/失真行为

注意默认情况下取值是按照插值来的，不过设想：如果单纯线性插值的量化网格正好错开了半个相位，那么他就会变成一个经典的滑动窗口平均低通——这显然不是我们希望看到的。我严格遵循了Math.round的吸附办法。

另外，纵轴采取Log0会导致崩溃，底噪失去数学意义。我加入了截止。


xxxxxxxxxx
const sampleNearest = (arr: number[], targetFreq: number) => {
    const maxIdx = arr.length - 1;
    const pos = (targetFreq / nyquist) * maxIdx;
    const idx = Math.round(pos);
    const val = arr[Math.min(idx, maxIdx)] || 0;
    return Math.max(val, 1e-10); 
};

6-4 采样窗口

二维采样区域不是无穷大的也并不完美，这导致我们对高动态范围的数组强行fft会导致严重的频谱泄露。为了抑制旁瓣一个自然的办法就是使用Hanning Window。但是考虑到一般情况下fft在raw分析中的适用范围（分析FPN等需要明确的特征峰，而且使用fft的对象一般动态都比较小例如黑场），我并没有这么做。不过我还是留了个开关——如果你们分析高动态的时候出现了严重的吉布斯效应，应当采取汉宁窗来获得一个更好的结果。

6-5 显示行为

6-5-1 坐标显示与频谱搬移

为确保严谨还是介绍一下。我没有采取很多2dFFT软件的频谱搬移操作。直流分量就在0处，如我网页的横坐标轴所示。

6-5-2 池化行为（对应图中POOLING选项卡）

就算屏幕完全完美，人眼也没法看到在0.25NF或者0.167/0.33NF上的近乎无穷小高度的尖峰。为了可选凸显这一部分尖峰来完成对FPN等的检查，我采取了可选池化显示的办法。页面中你可选多种池化模式。极值池化就是给出最大值池化、最小值池化的包络线，最大值池化就是只关心能量最大点的情况：


xxxxxxxxxx
const keysToCheck = [
    'c0_h', 'c0_v', 'c1_h', 'c1_v', 
    'c2_h', 'c2_v', 'c3_h', 'c3_v'
];
// ...
keysToCheck.forEach(key => {
    //...
    resultVal = Math.max(...values);
});

为了可复用我采取了浅拷贝，设置了2000个pool。在低于2000时，收回c0_h_max和c0_h_min。


xxxxxxxxxx
if (method === 'extremum') {
     return data.map(p => {
         const newP = { ...p }; // 浅
         keysToCheck.forEach(k => {
             if (newP[k] !== undefined) {
                 newP[`${k}_max`] = newP[k];
                 newP[`${k}_min`] = newP[k];
             }
         });
         return newP;
     });
}

当然为了保证绝对的严谨，可以配合parseval定理研究，我还是留了平均值池化的选项。

6-6 加速操作

在Bluestein算法中我选择了把FFT(kernel)这个不变量直接存起来。在处理每一行时，直接使用缓存中的FFT(kernel)进行频域点乘。具体代码我前面给了。

此外还有一点。对JavaScript来说，高频Garbage Collection会严重拖慢性能（计算24mp的四通道需要等半分有余，以我的手机为例）。

初始化的两个复用：


xxxxxxxxxx
constructor(size: number) {
    //...内存分配...
    this._real = new Float32Array(size);
    this._imag = new Float32Array(size);
    //...数学预计算...
    if (!this.isPowerOfTwo) {
        this.initBluestein(); 
    }
}

执行的两个复用：


xxxxxxxxxx
private transformBluestein(...) {
    const br = this._bReal![i]; 
    const aReal = fftM._real; 
    // ... 原地修改 buffer ...
}

这点确实需要注意。这样改可以增加CPU缓存命中率，提升数倍计算速度。

6+ fftlab简述

6-1 抽象/用法

提供了把两张raw图像相减的办法。相减成功后图片变得极亮（为精度，我在渲染/计算时给了充分高的黑电平），这样可以消除掉一部分固定模式的影响。帧相减不仅可以分析非固定模式噪声，还可以用来获取真正的（平场）噪声，完全剥离直流不均的影响。

进入fftlab后，选取可以出现这样的fft界面：

你可以使用Discrete Gauss Model来拟合出一条理想高斯滤波的线来简单模拟一下理想滤波的结果，或者使用积分工具来对某个区间内的σ进行积分（输出为功率）。汉宁窗可选，Row/Col独立选取可选。汉宁窗可选。

使用黑场能量谱，你可以在相当意义上严谨推出降噪前传感器数据。

举例：

一个符合工程的滤波一定会在平场中退化为：

归一化
响应线性

的一个滤波器。归一化是因为在平场中经过任何处理，CIS必须保证响应线性，因此降噪核必须归一化。

但这显然并不需要降噪核本身响应线性。降噪核本身的线性假定只是因为我们在分析一个平场，而平场本身的值域是很窄的，在这个区间内我们可以近似作线性假定。

对于一个线性系统，有PSD关系：

Syy(f)=Sxx(f)⋅|H(f)|2而对于归一化系统，有∑h=1⇒H(0)=1，所以：

Syy(0)=Sxx(0)⋅12=Sxx(0)

那么我们此时看趋近于DC的趋势：因为我的程序Jiangtherapee Online已经预先扣除了直流，这里我们要关注的是f->0的趋势。

物理直觉上也是直接的：只要我们看得足够低频，必然能看出核频率之外的性质。

但是在我这里我并没有采用2d-fft，因为2dfft看起来太费劲了。我是降维来看的结果。降维的结果就是0频本身也随着维度降下来了。那么很自然的行为就是我们重新积分展示回高维的零频预期。

按如下两个积分路径：

路径一：利用水平谱的 DC + 竖直谱的形状

σ2=SH(0)∫−1/21/2SV(f)SV(0)df

路径二：利用竖直谱的 DC + 水平谱的形状

σ2=SV(0)∫−1/21/2SH(f)SH(0)df

因为我在软件里面写过了均值池化算法，我们对均值池化的H\V能量谱进行离散求和。

σ2=SH[0]⋅SV[0]Mean(SV) 或者对称地

σ2=SV[0]⋅SH[0]Mean(SH) 。

使用积分工具：

H Mean是1.181，V Mean是1.191。

然后我们调节一下DISCRETE MODEL来找到这个零点噪声是多少。

$S(0)=1.282^2=1.643524$ 。

$S(H,0)=S(V,0)=1.64352$ 。

$R_H=\frac{1.181}{1.64352}=0.71858$ 。

$R_V=\frac{1.191}{1.64352}=0.72466$ 。

$S(H,0)$ $\frac{1.64352}{0.72466}=2.26791$ 。

$S(V,0)$ $\frac{1.64352}{0.71858}=2.28718$ 。

$\sqrt{2.26791\times2.28718}=2.27754$ 。

$\sqrt{2.27754}=1.50915$ 。

7 选区直方图功能

7-1 抽象

在选择了某个ROI（或者干脆是全图）之后，按ROI弹窗的HIST功能。

对于一般图片我建议你们选择1:1 Point-to-Point模式。

完整界面：

以柱状图形式给出，在1:1模式下有很高精度：

图：a6700拍照时采取13bit装14bit的办法，在我的软件中容易看出带间隔的码值。

图：d850在全部码值均可见明显prescale痕迹。

（码值放大1.165倍——每六位变成七位，断一位。这点我之前文章中说过我是如何发现的）

7-2 实现

程序创建四个数组。每个数组长度等于设定的binCount使用双重循环遍历每一个像素，利用奇偶性判断当前像素是R还是G还是B。计算Bin索引来判断像素值应当落入第几个桶中。

对于普通模式，这个索引为(像素值 - 最小值) / 步长；对于点对点模式，步长退化为1.


xxxxxxxxxx
    const bins = [
        new Uint32Array(actualBinCount), // Channel 0
        new Uint32Array(actualBinCount), // Channel 1
        new Uint32Array(actualBinCount), // Channel 2
        new Uint32Array(actualBinCount)  // Channel 3
    ];

    for (let r = startY; r < endY; r++) {
        const rowOffset = r * width;
        const rowIsOdd = r % 2; 
        
        for (let c = startX; c < endX; c++) {
            const colIsOdd = c % 2; 
            const pixelVal = rawData[rowOffset + c];
            
            // 过滤
            if (pixelVal < minVal || pixelVal > maxVal) continue;

            // Bayer 模式判断：(偶行*2 + 偶列)，位移运算即可
            // 0: (0,0) | 1: (0,1) | 2: (1,0) | 3: (1,1)
            const channelIndex = (rowIsOdd * 2) + colIsOdd;
            //算桶
            let binIndex;
            if (pointToPoint) {
                binIndex = Math.floor(pixelVal - minVal);
            } else {
                binIndex = Math.floor((pixelVal - minVal) / binWidth);
            }

            // 边界防护
            if (binIndex < 0) binIndex = 0;
            if (binIndex >= actualBinCount) binIndex = actualBinCount - 1;

            // 计数 +1
            bins[channelIndex][binIndex]++;
        }
    }

8 彩噪分析功能

8-1 抽象/用法

本功能是我首次引入的。通过选区后点击【COLORN】，我会给出彩噪在较均匀感知的Lab空间以及经典的xy空间的投影。不过他不会立刻显示。你需要填入这个机型的CM。如果你想分析你在某CM下渲染的结果的彩噪，那你直接填入你的CM即可。如果你想看看默认解码流程的彩噪，那需要借助exif功能来完成了。

有了JiangBridge，你很容易找到adobe推荐的CM：首先查看CalibrationIlluminant。找到D65对应的是2。

然后看到这里的ColorMatrix2

即为彩噪。分析结果也是非常强而有力的，下面是我之前的分析：

我用这个算法破解了A7R3/Z72的彩噪之谜。

8-2 算法构建

在相机空间中，我们如何分析彩噪呢？先不提如何量化彩噪，我打算率先考察一下【任何一次三通道采样在相机原生正交三维空间中的行为】。

首先我们应当弄明白：这三个方向互相有依赖性吗？这三个方向上的噪声行为都相同（比如某个常量加上信号值再开根）吗？如果这三个方向间没有关联，那么协方差等于零，这个三方向上的椭球就应当是正的；如果噪声行为不同、存在类似尼康的pre-scaling的行为（看我刚发的文章），那么某方向上的噪声分布就不应当简单建模。

我们先拿索尼A7R3开刀吧。看看他的噪声协方差矩阵怎么样。我取ettl6拍摄的图片，以xrite中纯色块（虚焦拍摄）分析：

\begin{matrix} 协方差矩阵 C_{c o v} = [\begin{matrix} 0.8792 & 0.0721 & - 0.0172 \\ 0.0721 & 1.8079 & - 0.0021 \\ - 0.0172 & - 0.0021 & 1.1824 \end{matrix}] \end{matrix}

\begin{matrix} 相关系数矩阵 C_{c o r r} = [\begin{matrix} 1.0000 & 0.0572 & - 0.0168 \\ 0.0572 & 1.0000 & - 0.0014 \\ - 0.0168 & - 0.0014 & 1.0000 \end{matrix}] \end{matrix}

看上去三通道分离度极高。当然因为我的光照不是那么完美，非对角元素难免偏大。实际情况大概是完全无相关性的。暗场也同样如此。毕竟几乎没有共享前端（一些qbayer是共享的），后端之间相互干扰很小，什么PRNU DSNU在这个色块尺度上也近乎于零。

从相机空间的RGB到展示空间的RGB，我们需要两个算子：一个是相机空间到光谱本身的重建算子，一个是将光谱投影到显示空间的算子。但是分析彩噪只考虑重建算子即可，可感知的彩噪规模往往远大于MacAdam ellipses和显示分辨率，所以我们不考虑显示空间这一层了。

8-2-1 1931xy模式

我们按老办法把这一切投射在D65中。标版、标准光照属于连续谱也很接近D65。我这么直接以某CM做没任何问题。

左图z72搭配A镜头，右图a7r3搭配B镜头。

8-2-2 LAB中的解耦

此时其实你并不可以直接在这个椭圆上面展示每点对应颜色的明度平均噪声。因为这个xy平面与明度耦合比较严重。与其在这个平面上展示一个奇怪的相交关系，我选择了将他们投影到一个独立与明度的平面上。

对于边缘的高饱和度点：他们本身信号就很弱，均方差自然会很小。

具体过程为：

\begin{matrix} S_{raw} = [R, G_{1}, G_{2}, B] σ_{raw} = [σ_{R}, σ_{G 1}, σ_{G 2}, σ_{B}] M_{XYZ \to RGB} \end{matrix}

\begin{matrix} S_{rgb} = [\begin{matrix} R \\ (G_{1} + G_{2}) / 2 \\ B \end{matrix}] Var (R) = σ_{R}^{2} \\ Var (G) = \frac{σ_{G 1}^{2} + σ_{G 2}^{2}}{4} \\ Var (B) = σ_{B}^{2} {Cov}_{rgb} = [\begin{matrix} Var (R) & 0 & 0 \\ 0 & Var (G) & 0 \\ 0 & 0 & Var (B) \end{matrix}] \end{matrix}

上面直接把标量乘矩阵，先不改了。结构是这样的。传播噪声：

\begin{matrix} M_{inv} = (M_{XYZ \to RGB})^{- 1} S_{xyz} = M_{inv} \cdot S_{rgb} {Cov}_{xyz} = M_{inv} \cdot {Cov}_{rgb} \cdot (M_{inv})^{T} \end{matrix}

在xyz空间定义椭球：

\begin{aligned} EigenDecomposition ({Cov}_{xyz}) & \to (Λ, V) \\ Λ & = {λ_{1}, λ_{2}, λ_{3}} \\ V & = [v_{1}, v_{2}, v_{3}] \\ P_{xyz} (θ, ϕ) & = S_{xyz} + V \cdot [\begin{array}{c} \sqrt{λ_{1}} \cos (θ) \sin (ϕ) \\ \sqrt{λ_{2}} \sin (θ) \sin (ϕ) \\ \sqrt{λ_{3}} \cos (ϕ) \end{array}] \\ {Points}_{xyz_surface} & = {P_{xyz} (θ, ϕ) ∣ θ \in [0, 2 π], ϕ \in [0, π]} \end{aligned}

变换到Lab空间：

\begin{aligned} f (t) & = {\begin{cases} \sqrt[3]{t} & if t > (\frac{6}{29})^{3} \\ \frac{1}{3} (\frac{29}{6})^{2} t + \frac{4}{29} & otherwise \end{cases} \\ L^{*} & = 116 \cdot f (Y / Y_{n}) - 16 \\ a^{*} & = 500 \cdot (f (X / X_{n}) - f (Y / Y_{n})) \\ b^{*} & = 200 \cdot (f (Y / Y_{n}) - f (Z / Z_{n})) \\ {Points}_{lab_surface} & = {f_{XYZ \to Lab} (p) ∣ p \in {Points}_{xyz_surface}} \end{aligned}

使用Rejection Sampling在内部随机生成N个点，得到点集。对其中任一点有

\begin{aligned} Neighbors (p_{i}) & = {p \in {Points}_{lab_surface} ∣ dist ((p_{a^{*}}, p_{b^{*}}), p_{i}) < r} \\ σ_{L^{*}} (a_{i}, b_{i}) & = StdDev ({p_{L^{*}} ∣ p \in Neighbors (p_{i})}) \\ {PointSize}_{i} & = σ_{L^{*}} (a_{i}, b_{i}) \times {Scale}_{L} \\ μ_{L^{*}} (a_{i}, b_{i}) & = Mean ({p_{L^{*}} ∣ p \in Neighbors (p_{i})}) \\ L_{i}^{'} & = μ_{L^{*}} (a_{i}, b_{i}) + L_{boost} \\ P_{lab_render} & = [L_{i}^{'}, a_{i}, b_{i}] \\ {Color}_{i} & = f_{Lab \to sRGB} (P_{lab_render}) \end{aligned}

因此我的这版本中，你可以直接通过圆点总面积直接识别出总明度噪声的绝对大小，且方便跨机型场景比较。看某点的绝对面积，即可估计出该点的明度噪声的感知大小。

8-3 实现

我跨网页传输数据使用的是借用localstorage的办法。键名为：

import_signal_r / import_sigma_r
import_signal_g1 / import_sigma_g1

数据写入完成后，程序立即调用window.open("https://y-g-jiang.github.io/XJZS.html", "_blank")打开目标网页。目标网页（XJZS.html）在加载时，会检查localStorage中是否存在这些特定的键值。如果存在，它就会自动读取并填入对应的输入框中，从而实现数据的无缝传递。

9 JiangBridge

9-1 抽象/用法

JiangBridge是一个本地主机。双击JiangBridge即可打开（他没有窗口，自动运行在后台——占用资源极低）。因为他是一个pyinstaller简易编译的可执行文件，他很容易被反编译。所以你们大可放心。

9-2 意义

JiangBridge单纯是为了提供给网页exiftool、adobe dng converter的【可以击穿沙盒】的快捷接口。尽管我的软件的绝大部分功能并不需要exiftool、converter的帮助，但这样可以让你们快速对照一些内容，也可以提高整个分析流程的效率。

从法律角度我并不能把adobe dng converter嵌入/打包进我的软件，因为他是闭源、商业的（尽管dng本身不是）。因此adobe dng converter只能从浏览器这个沙盒之外引用，我用JiangBridge通过一个服务器把这二者连接在一起。不过我还是推荐你们先下载一下：

如果这个链接打不开，浏览器里你直接搜索【adobe dng converter】即可，他是完全免费下载的。

https://helpx.adobe.com/camera-raw/using/adobe-dng-converter.html

下载adobe dng converter后你大概需要重启电脑来完成JiangBridge对ADC的识别。

10 PDR测试模块

10-1 抽象

我的算法的精度是远高于PTP官方结果的！这个看点列质量即可知道：

PTC，光子传输曲线，就是噪声随读出信号值变化的曲线。我的PDR测试模块的目的就是你和出来这样一个PTC曲线即为成功。

相机的噪声功率模型大致是读出噪声和散粒噪声以及图案噪声的功率之和（模型见下）。

我如此拟合：首先根据一系列拍摄的均匀格子结果，取大约十万~百万个格子，每格子取中心64*64块。这样每个格子就可以输出一组（平均信号值，噪声值）。接下来我直接取所有点并拟合得到这样一个表达式：

\begin{matrix} σ_{t o t a l} = \sqrt{σ_{r e a d}^{2} + \frac{S}{K} + (P \cdot S)^{2}} \end{matrix}

看上图，我拟合的精度还是非常高的。

10-2 拟合严谨性

我使用了固定步长的Theil-Sen估计器来完成稳健回归。这个估计器需要一个类似y=kx+b的目标，于是我想办法把上面的式子变成了直线：

相机的光子转移曲线（PTC）原本是一个弯曲的公式：

\begin{matrix} 总方差 = {读噪声}^{2} + \frac{信号}{K} + (PRNU \times 信号)^{2} \end{matrix}

$Y = \text{Slope} \cdot X + \text{Intercept}$ ：

\begin{matrix} \frac{净方差}{信号} = {PRNU}^{2} \cdot 信号 + \frac{1}{Gain} \end{matrix}

以此拟合。

如果我们把这个直线放到log坐标下看，他大致是这样：

注意我的默认拟合区间是两个蓝色虚线内的部分，做得很好。

下面我展开说一下镜头渐晕、屏幕缺陷等因素的影响以及我的应对。

经过我前几篇文章的讨论，我认为一切情况下的采样应当按照HSM法采样一切结果。中位数是一种退化。而不应当使用平均值——他受离群量影响太大（比如恰好有一个格子是坏的屏幕，那么就会产生离群量；暗角会产生波动，但是在下面我会证明在我的测量办法中你不用关心他）。具体来讲我错开半个相位拍摄了两组图片，以s52的iso100——他的读出噪声很小了。

我把2470dgdn2在50mm处缩到了f/5.6。在这种情况下，单像素σ造成的不确定性仍然远显著于频率为半个相位的波动。因此暗角不需考虑（光圈别太大就行）。

另外为了避免渐晕影响，有一个极端的办法就是把每个大选区直接拆成小份然后在表中分别绘制点。这固然没啥大问题，但是工作量有些重，效率低。更甚地，渐晕其实影响挺小的（后面我会给展示），但是如果点量增加、集成小块且拍摄量少，那么最终的更多点的线拟合处理（受限于log分块平均法的智能程度）并不一定能增加我这个行为的精度。那么我就把拆分选区这件事搁置了。另外标准差是有偏估计量（在这种情况下其实可以忽略），且一些低频波动被因此忽略掉了。算不算低频？约定俗成的办法还是囊括了这一部分。最后这样费劲分格子容易降低拟合效率。

截取图线的计算逻辑是按Log均匀分箱、中位数下采样。拟合的逻辑是对所有点计算并取HSM值。

10-3 用法简介

希望你们能完全按照我的顺序操作。建议切换到暗色模式。

在尽量暗的环境中拍摄。如果实在条件艰苦，可以搞个衣服罩住——pdr也可以精确，但其余拟合精度就不可作保证了。

虚焦拍摄标版！具体见我下面的例图

1 打开drlab

2 点击 generate test pic

3 点击start breathing

到一个恰当的速度。默认似乎有点快了，建议把速度划到0.5。

下面架上相机然后按F11进入全屏，对于24mp及以上相机，大致以此比例拍摄：

推荐镜头缩至f8以上，减少暗角。

图：24mp相机拍摄图卡大小（建议占比略大于这个尺寸，我这个有点极限）

既然屏幕亮度现在在慢速变化，你要做的就简单了：使用一系列快门速度（十分建议为12s、1s、1/20s、1/200s、1/500s、1/2000s）设置连续拍摄（相机的自动连续拍摄延时功能），每个快门速度拍摄50张照片。注意切换快门速度的时候保持相机机位不变，这会省你很多事。共耗费0.03w快门，我觉得造成不了什么影响。

曝光建议以【开启网页呼吸亮度拍摄图卡，12s的格子几乎全部过曝】为标准。

接下来我们随便挑一张曝光差不离的图像，拉到恰好可见所有格子的恰当亮度：点击【Mark Corners】并如图顺序选择四角的格子的中心（格子为10*10）。注意不要选到分隔延长线去了。顺序是左上、右上、左下、右下。第四个点就在右下角（不过无法演示，因为一经点击第四个点，这个界面就自动刷新了）。

在完成点击后，程序自动生成采样格子：

确保格子大致如上图的比例为佳（我这个有点极限了）。

程序自动输出csv：

这个csv不要管他，你可以直接删除。

下面我们直接点击Reuse Last&Batch：

然后我们一股脑传上去刚才我们开镜头盖拍的所有图片（我做了很多优化，确保他可以快速处理300张raw而不卡死），他会输出一个大csv。你也可以分批传，生成多个csv。

接下来我们从“Generate PTC”导入csv：

可以承受巨量的数据（甚至更多）：

导入后如图：

看界面：此时有三条虚线：红色虚线对应读出噪声RN，橙色对应图案FPN，灰色代表散粒。

此时左上角黑电平并不精确，他还是会在一定程度上影响测量结果的。我们拍摄1/50s的暗场（扣镜头盖），拍一张即可。然后用这个软件最普通的功能，划一片区域读出他们未扣黑电平的码值（我的软件可能自动扣黑电平，关掉即可。这点注意）。

拖动左、下两个蓝条可以更改曲线缩放。

如果不是很贴合，你可以直接上手调整左上的Gain拟合区间（Fit Range）或者直接改数。大部分情况下足够精确，你们没必要改。

总之目的就是让PTC的黑线拟合尽可能贴合中间这部分点集。之于为什么小概率上来不是特别贴合…一般情况就是黑电平飘了（这玩意不同快门速度下有微小漂移但是大部分机器并没有OB像素所以没法算），是小问题。

下面我们看右上角，已经出现了【PTP DR】以及【D55 DR】还有【FWC（D55 SAT CH）】。如果这二者差异不大（达到1ev才算大）那么我推荐选取D55 DR来评判相机画质。如图，S52的D55DR=11.19。之于PTP标准的DR，结果是11.21。

我们看一下ptp官方的结果是不是11.21：

看来他的PDR测得还算准。之于ptp给的ptc曲线就很难看了。

右上角也显示出了满阱FWC。注意这个是【D55下最快达到饱和通道的满阱】，是很有意义的一个指标。其他网站提供的可能是四通道平均满阱，也可能没有考虑到他高光部分是否线性（可能他干脆拿满量程算了）。我不可能如此不负责任的。

10-4 拓展

如果满阱附近有如此明显的SNR Knee：

那么大概率这个机器采取了dgo或者dcgo的输出办法。此时你应当让fit range的max线在从左往右第一个snr knee之前结束。

如果你希望只关心线性区段（一般不用考虑这么多）的话，我也做了适配。点击右上角【Select Linear Limit】并点击【展开最后3%点】，点击选择这个noise雪崩前的最后一点：

再按auto、再按一次Top3%回到初始视图：

10-5 加速操作

我在处理Apply&Calc事件的时候犯过一个错误，就是把计算数据+更新DOM+绘制UI集合在了一个EventLoop里。这个常见操作会导致很奇怪的卡顿（数分钟的主线程阻塞）。所以我把计算推迟到了下一帧。之于行列复用什么的就不展开了。

11 致谢

吴[2]、jackchou[3]、白龙[4]、伊娃[5]。