function SpatialMedian

% read an image
I = imread('cameraman.bmp');
J = imnoise(I,'speckle',0.02);

% call spatial median filter with 3-by-3 neighborhood for both data and
% regularization edges
U = DoSpatialMedian(J,0.5,1,1,50);

K = medfilt2(J,[3 3]);

figure(1),
subplot(2,2,1),imagesc(I),colormap(gray),axis image;
title('Cameraman');
subplot(2,2,2),imagesc(J),axis image;
title('Speckle noise added');
subplot(2,2,3),imagesc(K),axis image;
title('3-by-3 Median filter');
subplot(2,2,4),imagesc(U),axis image;
title('3-by-3 Spatial median filter');

function U = DoSpatialMedian(I,lambda,edge_conn,data_conn,tau)
% 3x3 neighborhood is the minimum allowed neighborhood here for both
% edge_conn and data_conn
% 1 = 3-by-3, 2 = 5-by-5, etc.
% tau: edge strength parameter


I = double(I);
minI = min(I(:));
I(:) = I - minI; % make I non-negative
[h,w] = size(I);

[X,Y] = meshgrid(1:w,1:h);
N = Y + h*(X-1);


nData=((2*data_conn+1)^2)*h*w;
nind=0;
for i=0:edge_conn,
    for j=0:edge_conn,
        if i==0||j==0,
            nind=nind+(h-i)*(w-j);
        else
            nind=nind+2*(h-i)*(w-j);
        end
    end
end
nind = nind-h*w;
tail = zeros(nind+nData,1);
head = zeros(nind+nData,1);
ind=0;
for i=0:edge_conn,
    for j=0:edge_conn,
        if i==0||j==0,
            if i==0 && j==0, % set data edges here
                tail(1:nData) = h*w+1;
                for k=1:(2*data_conn+1)^2,
                    head(ind+1:ind+h*w) = N(:);
                    ind = ind+h*w;
                end
            else % regularization edges
                tail(ind+1:ind+(h-i)*(w-j)) = N(1:end-i,1:end-j);
                head(ind+1:ind+(h-i)*(w-j)) = N(1+i:end,1+j:end);
                ind = ind + (h-i)*(w-j);
            end
        else % regularization edges
            tail(ind+1:ind+(h-i)*(w-j)) = N(1:end-i,1:end-j);
            head(ind+1:ind+(h-i)*(w-j)) = N(1+i:end,1+j:end);
            ind=ind+(h-i)*(w-j);
            tail(ind+1:ind+(h-i)*(w-j)) = N(1:end-i,1+j:end);
            head(ind+1:ind+(h-i)*(w-j)) = N(1+i:end,1:end-j);
            ind=ind+(h-i)*(w-j);
        end
    end
end

ub = lambda*ones(nData+nind,1);
for n=nData+1:nData+nind,
    ub(n) = exp(-((I(head(n))-I(tail(n)))/tau)^2);
end
lb = -ub;
supply = zeros(h*w+1,1);
obj = zeros(nData+nind,1);
ind=0;
for i=-data_conn:data_conn,
    for j=-data_conn:data_conn,
        obj(ind+1:ind+h*w)=-circshift(I,[i j]);
        ind = ind + h*w;
    end
end
tic;
[a, b, v] = cplexnet(obj, ub, lb, tail, head, supply, [], [], [], []);
toc;
U = reshape(v(1:end-1),h,w)-v(end);
U(:) = U + minI; % adjust for minI