Смещение изображения (базовая линия 2D) Удалить в Matlab

У меня есть изображение, похожее на это:

введите здесь описание изображения

и хотите удалить его базовую линию, чтобы она выглядела так:

введите здесь описание изображения

Изображение всегда разное, обычно имеет несколько пиков и имеет полное абсолютное смещение, а базовая поверхность наклонена/изогнута/нелинейна.

Я думал об использовании метода подгонки и вычитания базовой линии 1D для обычных спектров сигналов и создания изображения базовой линии 2D, а затем численного вычитания каждого из другого. Но в 2D никак не могу разобраться.

Это улучшенный вопрос, который я задавал раньше, но он должен быть более ясным.


person nolimits    schedule 08.10.2020    source источник
comment
Просто мозговой штурм, какие примерки/аппроксимации вы планируете использовать, когда мы вычисляем базовую линию 1D?   -  person MichaelTr7    schedule 10.10.2020
comment
Я предполагаю, что решение более тривиально, чем я думаю. Это может быть простая подгонка для всех линий x вдоль y. А может тогда и все y строк по x. чтобы поверхность подходила? Я не могу понять это в 2D.   -  person nolimits    schedule 11.10.2020
comment
Похоже, это можно решить, применив преобразование Фурье к изображению и удалив низкие частоты перед его инвертированием.   -  person Ander Biguri    schedule 12.10.2020


Ответы (3)


Мне кажется, что мы можем применить какой-то фильтр верхних частот, чтобы решить вашу проблему. Один из способов сделать это — использовать фильтр размытия (какой-нибудь низкочастотный фильтр) и вычесть размытую часть из оригинала (так называемое нерезкое маскирование). Таким образом, для фильтрации нижних частот вы можете использовать свертки с гауссовым фильтром или просто блочный фильтр. В качестве альтернативы вы также можете использовать медианный фильтр, что я и сделал здесь:

введите здесь описание изображения

%% 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');
person flawr    schedule 12.10.2020

В предыдущих ответах были предложены интересные математические методы удаления базовой линии. Но я предполагаю, что этот вопрос является продолжением вашего предыдущего вопросы, а под изображением вы подразумеваете, что ваши данные на самом деле являются изображением. Если это так, вы можете использовать методы обработки изображений, чтобы найти пики и сгладить области вокруг них.

введите здесь описание изображения

<сильный>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),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);

введите здесь описание изображения

введите здесь описание изображения

введите здесь описание изображения

Результаты для изображений, представленных в предыдущих вопросах:

введите здесь описание изображения

введите здесь описание изображения

person saastn    schedule 13.10.2020
comment
Большое спасибо саастн! Я попытался запустить ваш код, но меня смущает ошибка, которую я получаю из-за demoSteps. ›› removeBaseline(2020_10_5_12_57_7.png) Недостаточно входных аргументов. Ошибка в removeBaseline (строка 40), если demoSteps. ___ Как применить код к своему изображению? - person nolimits; 19.10.2020
comment
@nolimits, требуется два аргумента. Сначала оттенки серого, а затем логическое значение, указывающее, хотите ли вы, чтобы функция отображала этапы алгоритма. Добавил образец для этого в свой пост. - person saastn; 19.10.2020
comment
Спасибо! I2 — это LinkProp 1x1. Как я могу вернуть файл изображения? ›› imshow(I2) Ошибка при использовании imageDisplayValidateParams Ожидается, что входной номер 1, I, будет одного из следующих типов: числовой, логический - person nolimits; 19.10.2020
comment
@nolimits, извините, опять плохо. Образец был изменен. Но я очень не рекомендую вам использовать кусок кода, размещенный в Интернете, не зная, как он работает. Возвращается LinkProp, так что в последний рисунок, вы можете вращать все три графика вместе. - person saastn; 19.10.2020
comment
Спасибо большое! Я не так хорошо разбираюсь в кодировании, а ваш код выглядит таким сложным. Мне нужно время, чтобы обдумать это. Первый шаг состоял в том, чтобы заставить его работать, а затем посмотреть на каждый шаг, что он делает методом проб и ошибок. :) - person nolimits; 20.10.2020

У меня есть решение на Python, но, думаю, было бы несложно перенести его в MATLAB.

Он работает с кучей пиков. Я сделал несколько предположений, например, что есть несколько пиков. Работает с одним, но лучше, если пиков будет несколько. Пики могут перекрываться. Основным предположением, конечно же, является форма фона, но ее можно изменить, если существуют другие модели.
Основная идея состоит в том, чтобы вычесть функцию, но запрещающую отрицательные значения. Это делается с помощью функции дополнительных затрат, которую я оставляю дифференцируемой для минимизации. Как следствие, могут возникнуть проблемы со значениями, близкими к нулю. Такие случаи можно обработать, повторяя, насколько резкими будут дополнительные затраты. Можно начать с умеренного наклона, равного примерно единице, и повторить подбор с более крутым наклоном и начальными значениями из предыдущего подбора. Я сделал это на подобных проблемах, и это работает нормально. Технически не исключено, что после вычитания аппроксимирующего решения остаются небольшие отрицательные значения, поэтому для данных изображения потребуется дополнительный шаг, чтобы позаботиться об этом.

Вот код

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]

а также

исходные и исправленные данные

person mikuszefski    schedule 12.10.2020