amedfilt2_fixpt.m

来自「This file contains the material presente」· M 代码 · 共 126 行

function J = amedfilt2_fixpt(I) %#eml
% 2-D Adaptive Median Filter
% This filter ignores edge effects and boundary conditions, as such, the
% output is a cropped version of the original image, where the amount
% cropped is equal to the maximum window size vertically and horizontally.

% Define smax as a constant
smax = 9;

% Initialize Output Image (J)
J = I;

% Calculate valid region limits for filter
[nrows ncols] = size(I);
ll = ceil(smax/2);
ul = floor(smax/2);

% We can remove this line from the loops
window_ind = -ul:ul;

% Loop over the entire image ignoring edge effects
for rows = ll:nrows-ul
    for cols = ll:ncols-ul
        
        region = I(rows+window_ind,cols+window_ind);
        centerpixel = region(ll,ll);

        for s = 3:2:smax
            
            % We can collapse the ROI calculations into a single function
            [rmin,rmax,rmed] = roi_stats(region,smax,s);

            % adapt region size
            if rmed > rmin && rmed < rmax
                if centerpixel <= rmin || centerpixel >= rmax
                    J(rows,cols) = rmed;
                end

                % stop adapting
                break;
            end
        end
    end
end



function [rmin,rmax,rmed] = roi_stats(region,smax,s)
F = fimath(...
'RoundMode', 'Nearest',...
'OverflowMode', 'wrap',...
'ProductMode', 'KeepLSB', 'ProductWordLength', 16,...
'SumMode', 'KeepLSB', 'SumWordLength', 16,...
'CastBeforeSum', true);

persistent histogram;
persistent pmin;
persistent pmax;

if isempty(histogram) || s==3
    histogram = zeros(256,1,'uint16');
    pmin = rescale(fi(255,0,8,0,F),8);
    pmax = fi(0,0,8,8,F);
end

% Limits for ROI
ll = ceil(smax/2)-floor(s/2);
ul = ceil(smax/2)+floor(s/2);

if s==3
    for i = ll:ul
        for j = ll:ul
            val = region(i,j);
            ind = uint16(fi(256*region(i,j),0,8,0));
            histogram(ind+1) = histogram(ind+1)+1;
            
            if val > pmax
                pmax = val;
            end

            if val < pmin
                pmin = val;
            end
        end
    end
else
    for i = ll:ul
        for j = ll:ul
            if  i==ll || i==ul || j==ll || j==ul
                val = region(i,j);
                ind = uint16(fi(256*region(i,j),0,8,0));
                histogram(ind+1) = histogram(ind+1)+1;
            
                if val > pmax
                    pmax = val;
                end

                if val < pmin
                    pmin = val;
                end
            end
        end
    end
end

rmin = pmin;
rmax = pmax;
rmed = rescale(fi(255,0,8,0,F),8);

cs = histogram(1);
if cs >= s*s/2
    rmed = fi(0,0,8,8,F);
end

for i = 2:256
    if rmed == 0
        break;
    else
        cs = cs+histogram(i);

        if cs >= s*s/2
            rmed = rescale(fi(i-1,0,8,0,F),8);
            break
        end
    end    
end