在 MATLAB 中对图像进行颜色校正

Question

我有 2 张图片 im1 和 im2 如下所示。 im2 图片与im1 相同，但它们之间的唯一区别是颜色。 im1 每个颜色通道的 RGB 范围为 (0-255, 0-255, 0-255)，而 im2 的 RGB 范围为 (201-255, 126-255, 140-255)。我的练习是反转添加的效果，这样我就可以尽可能地将im2 恢复为im1。我有两个想法。首先是匹配它们的直方图，使它们具有相同的颜色。我使用histeq 进行了尝试，但它只恢复了图像的一部分。有什么方法可以将im2 的直方图更改为与im1 完全相同？第二种方法只是将每个像素值从im1 复制到im2，但这是错误的，因为它不会恢复原始图像状态。有什么恢复图片的建议吗？

Answer 1

下面的@sepdek 几乎建议了@NKN 提到的方法，但我将提供另一种方法。我可以建议的另一种选择是基于最小均方解进行颜色校正。这意味着我们可以假设将像素从im2 转换为im1 需要权重的线性组合。换句话说，给定一个 RGB 像素，其中它的红色、绿色和蓝色分量从损坏的图像 (im2) 中形成一个 3 x 1 向量，存在一些线性变换以在干净图像中获得其等效像素 (@ 987654331@)。换句话说，我们有这样的关系：

[R_im1]       [R_im2]
[G_im1] = A * [G_im2]
[B_im1]       [B_im2]

  Y     = A *   X

A 在这种情况下将是一个 3 x 3 矩阵。这实质上是执行矩阵乘法以获得输出校正像素。 im2 的输入 RGB 像素将是 X，im1 的输出 RGB 像素将是 Y。我们可以将其扩展到任意数量的像素，其中来自im1 和im2 的像素对将沿Y 和X 建立列。一般来说，这将进一步扩展X 和Y 到3 x N 矩阵。要找到矩阵A，您将找到最小均方误差解决方案。我不会深入探讨，但要找到A 的最优矩阵，这需要找到伪逆矩阵。在我们这里的例子中，A 将等于：

找到此矩阵A 后，您需要获取图像中的每个像素，对其进行整形以使其成为 3 x 1 向量，然后将A 与该向量相乘，就像上面的方法一样。您可能会问自己的一件事是，我需要从两张图像中获取什么样的像素才能使上述方法起作用？您必须遵守的一项准则是，您需要确保从两个图像之间的相同空间位置进行采样。因此，如果我们要在...比如说...第 4 行第 9 列抓取一个像素，您需要确保来自 im1 和 im2 的两个像素都来自同一行和同一列，并且它们被放置在X 和Y 的相同对应列中。

使用这种方法的另一个小提示是，您需要确保对图像中的大量像素进行采样以获得良好的解决方案，并且您还需要确保您的采样分布在整个图像上。如果我们将采样定位在一个小区域内，那么您将无法获得足够好的颜色分布，因此输出看起来不会很好。为问题选择多少像素取决于您，但根据经验，您会达到输出开始平稳且看不到任何差异的地步。出于演示目的，我在整个图像的随机位置选择了 2000 个像素。

因此，这就是代码的样子。我使用randperm 生成从1 到M 的随机排列，其中M 是图像中的像素总数。这些生成线性索引，以便我们可以从图像中采样并构建我们的矩阵。然后我们应用上面的等式找到A，然后取每个像素并与A 进行矩阵乘法以获得输出。废话不多说：

close all;
clear all;
im1 = imread('http://i.stack.imgur.com/GtgHU.jpg');
im2 = imread('http://i.stack.imgur.com/wHW50.jpg');

rng(123); %// Set seed for reproducibility
num_colours = 2000;
ind = randperm(numel(im1) / size(im1,3), num_colours);

%// Grab colours from original image
red_out = im1(:,:,1);
green_out = im1(:,:,2);
blue_out = im1(:,:,3);

%// Grab colours from corrupted image
red_in = im2(:,:,1);
green_in = im2(:,:,2);
blue_in = im2(:,:,3);

%// Create 3 x N matrices
X = double([red_in(ind); green_in(ind); blue_in(ind)]);
Y = double([red_out(ind); green_out(ind); blue_out(ind)]);

%// Find A
A = Y*(X.')/(X*X.');

%// Cast im2 to double for precision
im2_double = double(im2);

%// Apply matrix multiplication    
out = cast(reshape((A*reshape(permute(im2_double, [3 1 2]), 3, [])).', ...
          [size(im2_double,1) size(im2_double,2), 3]), class(im2));

---

让我们慢慢看一下这段代码。我正在直接从 StackOverflow 读取您的图像。之后，我使用rng 设置种子，以便您最终可以重现相同的结果。设置种子很有用，因为它允许您重现我所做的随机像素选择。我们生成这些线性索引，然后为im1 和im2 创建我们的3 x N 矩阵。查找A 正是我所描述的，但您可能不习惯rdivide / / 运算符。 rdivide 在运算符右侧找到逆，然后将其与左侧的任何内容相乘。这是一种更有效的计算方式，而不是分别计算右侧的倒数，然后在完成后与左侧相乘。事实上，MATLAB 会给你一个警告，说明避免单独计算逆，而你应该使用除法运算符。接下来，我将im2 转换为double 以确保精度，因为A 很可能是浮点值，然后将每个像素与A 相乘以计算结果。最后一行代码看起来很吓人，但如果你想知道我是如何得出这个的，我用它来创建复古风格的照片，这也需要像这种方法一样的矩阵乘法，你可以在这里阅读：@987654323 @ 。 out 存储我们的最终图像。运行这段代码并显示out 的样子后，我们得到了：

现在，输出看起来完全乱码，但颜色分布或多或少模仿了输入原始图像的样子。我有一些解释为什么会这样：

1. 存在量化噪声。如果你看一下最终的图像，到处都是各种白点。这可能是由于压缩图像时引入的量化误差。由于量化，应该在图像之间映射到相同颜色的像素会有轻微的变化，这给了我们斑点
2. im2 中不止一种颜色映射到im1。如果从im2 映射到im1 的颜色不止一种，则在给定单个像素的情况下，与矩阵A 的线性乘法不可能为im1 生成一种以上的颜色在im2。相反，最小均方解决方案将尝试生成一种颜色，使误差最小化，并为您提供可能的最佳颜色。这可能是由于这个确切原因，图像的面部和其他精细细节被遮挡的方式。
3. 图像有噪点。您的im2 不完全干净。我还可以在所有频道上看到各种椒盐噪声。这种方法的一个不好的地方是，如果你的图像受到噪声的影响，那么这种方法将无法正确正确地重建原始图像。您的图像只能被错误的颜色映射损坏。如果引入了任何其他类型的图像噪声，那么当您尝试基于噪声图像重建原始图像时，此方法肯定不起作用。嘈杂图像中有一些像素从未出现在原始图像中，因此您将无法将其恢复到以前的状态！

如果你想看一下原始图像和输出图像之间每个通道的直方图，这就是我们得到的：

我用来生成上图的代码是：

names = {'Red', 'Green', 'Blue'};
figure;
for idx = 1 : 3
    subplot(3,2,2*idx - 1);
    imhist(im1(:,:,idx));
    title([names{idx} ': Image 1']);
end

for idx = 1 : 3
    subplot(3,2,2*idx);
    imhist(out(:,:,idx));
    title([names{idx} ': Output']);
end

左侧显示原始图像的红色、绿色和蓝色直方图，而右侧显示重建图像的相同直方图。您可以看到总体形状或多或少模仿了原始图像，但始终存在一些尖峰 - 很可能归因于量化噪声和两个图像颜色之间的非唯一映射。

---

总而言之，这是我能做的最好的事情，但我认为这就是练习的全部意义......表明这是不可能的。

有关如何执行色彩校正的更多信息，请查看Richard Alan Peters' II Digital Image Processing slides on colour correction。这是我开始的，如何计算A 的推导可以在他的幻灯片中找到。也许你可以在你未来的工作中使用他所说的一些内容。

祝你好运！

Answer 2

您似乎需要一个缩放函数来将im2 的值映射到im1 的值。 这相当简单，您可以编写一个缩放函数以使其适用于任何此类情况。 一个基本的缩放映射将按如下方式工作：

out_value = min_output + (in_value - min_input) * (outrange / inrange)

假设有一个输入值in_value 在值inrange=max_input-min_input 的范围内，并且映射导致输出值out_value 在outrange=max_output-min_output 的范围内。我们还需要考虑最小输入和输出范围界限（min_input 和 min_output）才能获得正确的映射。
例如，请参见以下缩放函数代码：

%
% scale the values of a matrix using a set of limits
% possible ways to use:
% y = scale( x, in_range, out_range)  --> ex. y = scale( x, [8 230], [0 255])
% y = scale( x, out_range)            --> ex. y = scale( x, [0 1])
% 
function y = scale( x, varargin );
    if nargin<2,
        error([upper(mfilename),':: Syntax: y=',mfilename,'(x[,in_range],out_range)']);
    end;
    if nargin==2,
        inrange=[min(x(:)) max(x(:))];      % compute the limits of the input variable
        outrange=varargin{1};               % get the output limits from the arguments
    else
        inrange=varargin{1};                % get the input limits from the arguments
        outrange=varargin{2};               % get the output limits from the arguments
    end;

    if diff(inrange)==0,                    % row or column vector matrix or scalar
        % just do a clipping...
        if x>=outrange(2), 
            y=outrange(2); 
        elseif x<=outrange(1), 
            y=outrange(1); 
        else 
            y=x; 
        end;
    else
        % actually scale the data
        % using: out = min_output + (x-min_input) * (outrange / inrange)
        y = outrange(1) + (x-inrange(1))*abs(diff(outrange))/abs(diff(inrange));
    end;

此函数获取值矩阵并将它们缩放到所需范围。
在您的情况下，它可以按如下方式使用（变量img 是缩放的im2）：

for i=1:size(im1,3),        % for each of the input/output image channels
    output_range = [min(min(im1(:,:,i))) max(max(im1(:,:,i)))];
    img(:,:,i) = scale( im2(:,:,i), output_range);
end;

这样im2 被缩放到im1 的值范围，一次一个通道。输出变量img 应该是所需的变量。