Двумерный согласованный фильтр

Question

Я хочу реализовать двумерный согласованный фильтр для извлечения кровеносных сосудов в соответствии со статьей «Обнаружение кровеносных сосудов на изображениях сетчатки с использованием двумерных согласованных фильтров» Чаудхури и др., IEEE Trans.по медицинской визуализации, 1989 г. ( есть PDF на веб-сайте автора ).

Краткое описание: поперечное сечение кровеносного сосуда имеет гауссовское распределение, и поэтому я хочу использовать гауссовское соответствиефильтр для увеличения SNR.Такое ядро ​​может быть математически выражено как:

K(x,y) = -exp(-x^2/2*sigma^2)  for |x|<3*sigma,  |y|<L/2

L здесь - длина судна с фиксированной ориентацией.Экспериментально sigma=1.5 и L = 7.

Мой код MATLAB для этой части:

s = 1.5; %sigma
t = -3*s:3*s;
theta=0:15:165; %different rotations
%one dimensional kernel
x = 1/sqrt(6*s)*exp(-t.^2/(2*s.^2));

L=7;
%two dimensional gaussian kernel
x2 = repmat(x,L,1);

Рассмотрим реакцию этого фильтра для пикселя, принадлежащего фоновой сетчатке.Предполагая, что фон имеет постоянную интенсивность с нулевым средним аддитивным гауссовским белым шумом, ожидаемое значение выходного сигнала фильтра в идеале должно быть равно нулю.Поэтому ядро ​​свертки модифицируется путем вычитания среднего значения s(t) из самой функции.Среднее значение ядра определяется следующим образом: m = Sum(K(x,y))/(number of points).

Таким образом, сверточная маска, используемая в этом алгоритме, определяется следующим образом: K(x, y) = K(x,y) - m.

Мой код MATLAB:

m = sum(x2(:))/(size(x2,1)*size(x2,2));
x2 = x2-m; 

Судно может быть ориентировано под любым углом 0<theta<180, и согласованный отклик фильтра максимален, когда он выровнен по theta+- 90 (распределение поперечного сечения по Гауссу, а не само судно).

Таким образом, нам нужно повернуть согласованный фильтр 12 раз с шагом 15 градусов.

Мой код MATLAB прилагается, но я не получаю желаемого результата.Любая помощь приветствуется.

%apply rotated matched filter on image
r = {};
for k = 1:12
    x3=imrotate(x2,theta(k),'crop');%figure;imagesc(x3);colormap gray;   
    r{k}=conv2(img,x3);
end
w=[];h = zeros(584,565);
for i = 1:565
    for j = 1:584
        for k = 1:32
            w= [w ,r{k}(j,i)];
        end
        h(j,i)=max(abs(w));
        w = [];
    end
end
%show result
figure('Name','after matched filter');imagesc(h);colormap gray

Для вращения я использовал imrotate, что мне кажется более разумным, но в статье все иначе: предположим, p=[x,y] будет дискретной точкой в ​​ядре.Для вычисления коэффициентов в повернутом ядре у нас есть [u,v] = p*Rotation_Matrix.

Rotation_Matrix=[cos(theta),sin(theta);-sin(theta),cos(theta)]

И ядро:

K(x,y) = -exp(-u^2/2*s^2)

Но новое ядро ​​больше не имеет гауссовой формы.Использование imrotate сохраняет гауссову форму.Так в чем же преимущество использования матрицы вращения?

Входное изображение:

Выход:

Согласованная фильтрация помогает увеличить SNR, но также усиливается фоновый шум.Правильно ли я использую imrotate для вращения ядра?Моя главная проблема с матрицей вращения, почему и каков правильный код для ее реализации.

Answer 1

Причина построения фильтра из его аналитического выражения для каждого поворота, а не с использованием imrotate, заключается в том, что экстент фильтра не является круглым, и, следовательно, вращение приводит к появлению «новых» значений пикселей и выталкивает некоторые другие пиксели изядро.Кроме того, вращение ядра, построенного так, как здесь (плавный переход вдоль одного направления, шаг по краю вдоль другого измерения), требует разных методов интерполяции вдоль каждого измерения, что imrotate не может сделать.Полученное повернутое ядро ​​всегда будет неправильным.

Обе эти проблемы можно легко увидеть при отображении созданного вами ядра вместе с двумя повернутыми версиями:

На этом экране появляются дополнительные проблемы: ядро ​​не центрируется на пикселе, что приводит к смещению выходного сигнала на половину пикселя.

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

Вращенные ядра могутбыть построено следующим образом:

m = max(ceil(3*s),(L-1)/2);
[x,y] = meshgrid(-m:m,-m:m); % non-rotated coordinate system, contains (0,0)
t = pi/6;                    % angle in radian
u = cos(t)*x - sin(t)*y;     % rotated coordinate system
v = sin(t)*x + cos(t)*y;     % rotated coordinate system
N = (abs(u) <= 3*s) & (abs(v) <= L/2);   % domain
k = exp(-u.^2/(2*s.^2));     % kernel
k = k - mean(k(N));
k(~N) = 0;                   % set kernel outside of domain to 0

Это результат для трех поворотов, использованных в примере выше (серый цвет по краям ядра соответствует значению 0, черные пиксели имеют отрицательное значение):

Другая проблема заключается в том, что вы используете conv2 со стандартной формой вывода 'full', вы шоздесь следует использовать 'same', чтобы выходные данные фильтра совпадали с входными данными.

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

img = im2double(rgb2gray(img));

s = 1.5; %sigma
L = 7;
theta = 0:15:165; %different rotations

out = zeros(size(img));

m = max(ceil(3*s),(L-1)/2);
[x,y] = meshgrid(-m:m,-m:m); % non-rotated coordinate system, contains (0,0)
for t = theta
   t = t / 180 * pi;        % angle in radian
   u = cos(t)*x - sin(t)*y; % rotated coordinate system
   v = sin(t)*x + cos(t)*y; % rotated coordinate system
   N = (abs(u) <= 3*s) & (abs(v) <= L/2); % domain
   k = exp(-u.^2/(2*s.^2)); % kernel
   k = k - mean(k(N));
   k(~N) = 0;               % set kernel outside of domain to 0

   res = conv2(img,k,'same');
   out = max(out,res);
end

out = out/max(out(:)); % force output to be in [0,1] interval that MATLAB likes
imwrite(out,'so_result.png')

Я получаю следующий вывод: