问题描述
,并希望删除其基础基线,使其看起来像:
图像总是不同的,通常具有一些峰,并且具有绝对的总偏移量,并且底面是倾斜/弯曲/非线性的。
我正在考虑使用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]
和