在Matlab中移除图像偏移二维基线

如何解决在Matlab中移除图像偏移二维基线

我有一张类似的图片：

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

图像总是不同的，通常具有一些峰，并且具有绝对的总偏移量，并且底面是倾斜/弯曲/非线性的。

我正在考虑使用1D基准线拟合和减法技术处理常见信号频谱，并创建2D基准线图像，然后在数值上彼此相减。但是我不能完全以2D方式解决它。

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

解决方法

在我看来，我们可以采用某种高通滤波器来解决您的问题。一种方法是使用模糊滤波器（某种低通滤波器），然后从原始图像中减去模糊部分（称为“不清晰蒙版”）。因此，对于低通滤波，可以使用带高斯的卷积运算器或仅使用普通盒滤波器。另外，您也可以使用中值滤波器，这就是我在这里所做的：

%% 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,z0_);title('recovered baseline');
subplot(223);mesh(x,z0_-z0);title('error');

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

,

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

1。预处理

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

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矩阵在峰值处包含一个值，在其他区域包含零。到目前为止，我们可以通过在主图像中乘以该矩阵来删除基线。但是最好先软化发现区域的边缘，以使峰的边缘不会突然降至零。

FC = imfilter(double(P),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,'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),'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,2);
    surf(F,'none'); hold on;
    contour3(F,2)
    h3 = subplot(1,3);
    surf(T,'none');
    
    hlink = linkprop([h1 h2 h3],{'CameraPosition','CameraUpVector','xlim','ylim','zlim','clim'});
    set(h1,[0 max(I(:))])
    set(h1,[0 size(I,1)])
    set(h1,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),20),'replicate');
[X,Y] = meshgrid(linspace(0,1,size(I,2)),linspace(0,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);

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

,

我在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,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,ff = params
    ###subtract such that values become as small as possible
    ###
    diffarray = za - second_order( xa,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,65,2,0.3 )
# ~VV = VV + peak( XX,55,4,0.5 )
# ~VV = VV + peak( XX,-1,0.4 )
# ~VV = VV + peak( XX,-3,6,0.7 )

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

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

提供

<< [0.092135 0.18974991 0.06144773 0.37054049 2.05096262 0.88314363]

和