在 Matlab 中删除图像偏移（2D 基线）

Posted 2023-02-23

【中文标题】在 Matlab 中删除图像偏移（2D 基线）

【英文标题】：Image Offset (2D Baseline) Remove in Matlab

【发布时间】：2021-01-23 00:22:33

【问题描述】：

我有一张类似的图片：

并希望删除其底层基线，使其看起来像：

图像总是不同的，通常有一些峰，有一个总的绝对偏移量和一个倾斜/弯曲/非线性的基础表面。

我正在考虑使用 1D 基线拟合和减法技术来处理常见信号光谱，并创建 2D 基线图像，然后以数字方式从另一个中减去。但在 2D 中我无法完全理解它。

这是我之前提出的一个改进的问题，但这个问题应该更清楚。

【问题讨论】：

只是头脑风暴，当我们计算一维基线时，您考虑采用哪种拟合/近似值？

我认为解决方案比我想象的要简单。它可以简单地拟合沿 y 的所有 x 线。然后可能沿着 x 的所有 y 线。表面贴合？我无法在 2D 中理解这一点。

这看起来可以通过对图像应用傅里叶变换并在反转之前去除低频来解决

【参考方案1】：

在我看来，我们可以应用某种高通滤波器来解决您的问题。这样做的一种方法是使用模糊滤镜（某种低通滤镜），并从原始滤镜中减去模糊部分（称为“锐化遮罩”）。因此，对于低通滤波，您可以使用带有高斯的卷积或仅使用普通的箱形滤波器。或者，您也可以使用我在这里所做的中值滤波器：

%% setup
v = 0:0.01:1;
[x,y] = meshgrid(v, v);
z0 = cos(pi*x).*cos(pi*y);z = z0; % "baseline" surface

pks = [1,1; 3,3; 7,5; 2,8; 7, 3]/10;% add 5 peaks
for p=pks';
   z = z + exp(-((x-p(1)).^2 + (y-p(2)).^2)/0.02.^2);
end

subplot(221);mesh(x,y,z);title('data');
%% recover "baseline"
z0_ = medfilt2(z, [1,1]*20, 'symmetric'); % low pass filter of your choice

subplot(222);mesh(x,y,z0_);title('recovered baseline');
subplot(223);mesh(x,y,z0_-z0);title('error');

%% subtract recovered baseline
subplot(224);mesh(x,y,z-z0_);title('recovered baseline removed');

【参考方案2】：

以前的答案提出了一些有趣的数学方法来删除基线。但我想这个问题是您previousquestions 的延续，而“图像”是指您的数据实际上是图像。如果是这样，您可以使用图像处理技术来找到峰值并将其周围的区域变平。

1.预处理

在应用不同的过滤器之前，最好将像素值映射到一定范围内。这样我们可以更好地控制过滤器所需参数的值。 首先，我们将图像数据类型转换为双精度，用于像素值为整数的情况。

I = double(I);

然后，通过应用平均滤波器，我们减少了图像中的噪声。

SI = imfilter(I,fspecial('disk',40),'replicate');

最后，我们将所有像素的值映射到从零到一的范围内。

NI = SI-min(SI(:));
NI = NI/max(NI(:));

2。细分

准备好图像后，我们可以识别每个峰所在的部分。为此，我们首先计算图像梯度。

G = imgradient(NI,'sobel');

然后我们识别出图像中斜率较高的部分。由于“高坡度”在不同的图像中可能有不同的含义，因此我们使用graythresh函数将图像分为低坡度和高坡度两部分。

SA = im2bw(G, graythresh(G));

上一步的分割区域可能有几个问题：

被归类为高坡度区域一部分的小连续分量可能仅由噪声引起。因此，应移除面积小于阈值的组件。 由于斜率在峰的顶部达到零，因此在上一步中发现的组件中可能存在孔。 沿边界的峰的斜率不一定相同，发现的区域可能具有不规则的形状。一种解决方案是通过将它们替换为凸形大厅来扩展它们。

    [L, nPeaks] = bwlabel(SA);
    minArea = 0.03*numel(I);
    P = false(size(I));
    for i=1:nPeaks
        P_i = bwconvhull(L==i);
        area = sum(P_i(:));
        if (area>minArea)
            P = P|P_i;
        end
    end

3。删除基线

在上一步计算的P 矩阵在峰值处包含值 1，在其他区域包含值为零。到目前为止，我们可以通过在主图像中乘以这个矩阵来删除基线。但最好先软化发现区域的边缘，以免峰的边缘突然降为零。

FC = imfilter(double(P),fspecial('disk',50),'replicate');
F = I.*FC;

您还可以移动边缘图像最少的峰。

E = bwmorph(P, 'remove');
o = min(I(E));
T = max(0, F-o);

上述所有步骤都在一个函数中

function [hlink, T] = removeBaseline(I, demoSteps)
% converting image to double
I = double(I);
% smoothing image to reduce noise
SI = imfilter(I,fspecial('disk',40),'replicate');
% normalizing image in [0..1] range
NI = SI-min(SI(:));
NI = NI/max(NI(:));
% computng image gradient
G = imgradient(NI,'sobel');
% finding steep areas of the image
SA = im2bw(G, graythresh(G));
% segmenting image to find regions covered by each peak
[L, nPeaks] = bwlabel(SA);
% defining a threshold for minimum area covered by each peak
minArea = 0.03*numel(I);
% filling each of the regions, and eliminating small ones
P = false(size(I));
for i=1:nPeaks
    % finding convex hull of the region
    P_i = bwconvhull(L==i);
    % computing area of the filled region
    area = sum(P_i(:));
    if (area>minArea)
        % adding the region to peaks mask
        P = P|P_i;
    end
end
% applying the average filter on peaks mask to compute coefficients
FC = imfilter(double(P),fspecial('disk',50),'replicate');
% Removing baseline by multiplying the coefficients
F = I.*FC;
% finding edge of peaks
E = bwmorph(P, 'remove');
% finding minimum value of edges in the image 
o = min(I(E));
% shifting the flattened image
T = max(0, F-o);

if demoSteps
    figure
    subplot 231, imshow(I, []); title('Original Image');
    subplot 232, imshow(SI, []); title('Smoothed Image');
    subplot 233, imshow(NI); title('Normalized in [0..1]');
    subplot 234, imshow(G, []); title('Gradient of Image');
    subplot 235, imshow(SA); title('Steep Areas');
    subplot 236, imshow(P); title('Peaks');
    
    figure;
    subplot 221, imshow(FC); title('Flattening Coefficients');
    subplot 222, imshow(F, []); title('Base Line Removed');
    subplot 223, imshow(E); title('Peak Edge');
    subplot 224, imshow(T, []); title('Final Result');
    
    figure
    h1 = subplot(1, 3, 1);
    surf(I, 'edgecolor', 'none'); hold on;
    contour3(I, 'k', 'levellist', o, 'linewidth', 2)
    h2 = subplot(1, 3, 2);
    surf(F, 'edgecolor', 'none'); hold on;
    contour3(F, 'k', 'levellist', o, 'linewidth', 2)
    h3 = subplot(1, 3, 3);
    surf(T, 'edgecolor', 'none');
    
    hlink = linkprop([h1 h2 h3],'CameraPosition','CameraUpVector', 'xlim', 'ylim', 'zlim', 'clim');
    set(h1, 'zlim', [0 max(I(:))])
    set(h1, 'ylim', [0 size(I, 1)])
    set(h1, 'xlim', [0 size(I, 2)])
    set(h1, 'clim', [0 max(I(:))])
end
end

使用包含多个噪声峰值的图像执行该函数：

close all; clc; clear variables;
I = abs(peaks(1200));
J1 = imnoise(ones(size(I))*0.5,'salt & pepper', 0.05);
J1 = imfilter(double(J1),fspecial('disk',20),'replicate');
[X, Y] = meshgrid(linspace(0, 1, size(I, 2)), linspace(0, 1, size(I, 1)));
J2 = X.^3+Y.^2;
I = max(I, 2*J2) + 5*J1;
lp3 = removeBaseline(I, true);

调用从文件中读取图像的函数：

I = rgb2gray(imread('imagefile.jpg'));
[~, I2] = removeBaseline(I, true);

先前问题中提供的图像的结果：

【讨论】：

非常感谢！我尝试运行您的代码，但我对由于“demoSteps”而得到的错误感到困惑。 >> removeBaseline("2020_10_5_12_57_7.png") 输入参数不足。如果 demoSteps.removeBaseline（第 40 行）出错。 ___ 如何为我的图像应用代码？

@nolimits，它需要两个参数。首先是 灰度 图像，其次是布尔值，指示您是否希望函数绘制算法步骤。在我的帖子中添加了一个示例。

谢谢！ I2 是一个 1x1 LinkProp。如何取回图像文件？ >> imshow(I2) 使用 imageDisplayValidateParams 时出错 预期输入数字 1，I，是以下类型之一：数字、逻辑

@nolimits，对不起，又是我的错。修改了样本。但我真的不鼓励你在不知道它是如何工作的情况下使用发布在 Internet 上的一段代码。 The LinkProp is returned so that in the last figure, you can rotate all three plots together.

非常感谢！我不擅长编码，你的代码看起来很复杂。我需要时间来解决它。第一步是让它运行，然后通过反复试验查看每个步骤的作用。 :)

【参考方案3】：

我在 Python 中有一个解决方案，但猜想将其传输到 MATLAB 不会太复杂。

它适用于一堆峰。不过，我做了一些假设，比如有几个峰值。它适用于一个，但如果有几个峰值会更好。峰可能重叠。主要假设当然是背景的形状，但如果存在其他模型，则可以对其进行修改。 主要思想是减去一个函数，但禁止负值。这是通过额外的成本函数完成的，为了最小化，我保持可微分。因此，接近零的值可能存在问题。这种情况可以通过迭代额外成本的急剧程度来处理。可以从大约 1 的适度斜率开始，然后以更陡峭的斜率和前一次拟合的起始值重新进行拟合。我已经在类似的问题上做到了，它工作正常。从技术上讲，不排除减去拟合解后存在小的负值，因此对于图像数据，需要额外的步骤来处理这一点。

这里是代码

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy.optimize import least_squares

def peak( x,y, a, x0, y0, s):
    """
    Just a symmetric peak for testing
    """
    return a * np.exp( -( (x - x0 )**2 + ( y - y0 )**2 ) / 2 / s**2 )

def second_order( xx, yy, aa, bb, cc, dd, ee, ff ):
    """
    Assuming that the base can be approximated by a second order equation
    generalization to higher orders should be straight forward
    """
    out = aa * xx**2 + 2 * bb * xx * yy + cc * yy**2 + dd * xx + ee * yy + ff
    return out


def residual_function( params, xa, ya, za, extracost, slope ):
    """
    cost function. Calculates difference from zero-plane
    with ultra high cost for negative values.
    previous solutions to similar problems have shown that sometimes 
    the optimization process has to be iterated with increasing 
    parameter slope ( and maybe extracost )
    """
    aa, bb, cc, dd, ee, ff = params
    ###subtract such that values become as small as possible
    ###
    diffarray = za - second_order( xa, ya, aa, bb, cc, dd, ee, ff )
    diffarray = diffarray.flatten( )
    ### BUT high costs for negative values
    cost = np.fromiter( (  -extracost * ( np.tanh( slope * x ) - 1 ) / 2.0 for x in diffarray ), np.float )
    return np.abs( cost ) + np.abs( diffarray )


### some test data
xl = np.linspace( -3, 5, 50 )
yl = np.linspace( -1, 7, 60 )

XX, YY = np.meshgrid( xl, yl )
 

VV = second_order( XX, YY, 0.1, 0.2, 0.08, 0.28, 1.9, 1.3 ) 
VV = VV + peak( XX, YY, 65, 1, 2, 0.3 )
# ~VV = VV + peak( XX, YY, 55, 3, 4, 0.5 )
# ~VV = VV + peak( XX, YY, 55, -1, 0, 0.4 )
# ~VV = VV + peak( XX, YY, 55, -3, 6, 0.7 )

### the optimization
result = least_squares(residual_function,  x0=( 0.0, 0.0, 0.00, 0.0, 0.0, 0 ), args=( XX, YY, VV, 1e4, 50 ) )
print result
print result.x
subtractme = second_order( XX, YY, *(result.x) ) 
nobase = VV - subtractme

### plotting
fig = plt.figure()
ax = fig.add_subplot( 1, 2, 1, projection='3d' )
ax.plot_surface( XX, YY, VV)
bx = fig.add_subplot( 1, 2, 2, projection='3d' )
bx.plot_surface( XX, YY, nobase)
plt.show()

提供

<< [0.092135 0.18974991 0.06144773 0.37054049 2.05096262 0.88314363]

和