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

Question

我有一张类似的图片：

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

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

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

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

Answer 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');

Answer 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);

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

Answer 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]

和