Retinal-Vessel-Segmentation/MatchFilterAndGaussDerKernel.m at master · farkoo/Retinal-Vessel-Segmentation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
% MF and FDOG kernel construction
% Inputs:
%       sigma   = scale value
%       YLength = length of neighborhood along y-axis
%       theta   = orientation
%       type    = MF or FDOG
% output:
%       kernel  = Gaussian kernel
function [kernel] = MatchFilterAndGaussDerKernel(sigma, YLength,theta, type)

    if type == 0
        widthOfTheKernel = ceil(sqrt((6*ceil(sigma)+1)^2 + YLength^2));
        if mod(widthOfTheKernel,2) == 0
            widthOfTheKernel = widthOfTheKernel + 1;
        end
        halfLength = (widthOfTheKernel - 1) / 2;
        row = 1;

        for y = halfLength:-1:-halfLength
            col = 1;
            for x = -halfLength:halfLength
                xPrime = x * cos(theta) + y * sin(theta);
                yPrime = y * cos(theta) - x * sin(theta);
                if abs(xPrime) > 3.5*ceil(sigma)
                    matchFilterKernel(row,col) = 0;
                elseif abs(yPrime) > (YLength-1) / 2
                    matchFilterKernel(row,col) = 0;
                else
                    matchFilterKernel(row,col) = -exp(-.5*(xPrime/sigma)^2)/(sqrt(2*pi)*sigma);
                end
                col = col + 1;
            end
            row = row + 1;
        end
        mean=sum(sum(matchFilterKernel))/sum(sum(matchFilterKernel<0));
        for i = 1:widthOfTheKernel
            for j =1:widthOfTheKernel
                if matchFilterKernel(i,j) < 0
                    matchFilterKernel(i,j)=matchFilterKernel(i,j)- mean;
                end
            end
        end
        kernel = matchFilterKernel;
    else
        widthOfTheKernel = ceil(sqrt( (6*ceil(sigma)+1)^2 + YLength^2));
        if mod(widthOfTheKernel,2) == 0
            widthOfTheKernel = widthOfTheKernel + 1;
        end
        halfLength = (widthOfTheKernel - 1) / 2;
        row = 1;
        for y = halfLength:-1:-halfLength
            col = 1;
            for x = -halfLength:halfLength
                xPrime = x * cos(theta) + y * sin(theta);
                yPrime = y * cos(theta) - x * sin(theta);
                if abs(xPrime) > 3.5*ceil(sigma)
                    GaussDerivativeKernel(row,col) = 0;
                elseif abs(yPrime) > (YLength-1) / 2
                    GaussDerivativeKernel(row,col) = 0;
                else
                    GaussDerivativeKernel(row,col)= -exp(-.5*(xPrime/sigma)^2)*xPrime/(sqrt(2*pi)*sigma^3);
                end
                col = col + 1;
            end
            row = row + 1;
        end
        kernel = GaussDerivativeKernel;
    end

end