function f= dd(r,Smax)
if (Smax <= 1) | (Smax/2 == round(Smax/2)) | (Smax ~= round(Smax))
error('SMAX must be an odd integer > 1.')
end
[M, N] = size(r);
f = r;
f(:) = 0;
alreadyProcessed = false(size(r));
for k = 3:2:Smax
zmin = ordfilt2(r, 1, ones(k, k), 'symmetric');
zmax = ordfilt2(r, k * k, ones(k, k), 'symmetric');
zmed = medfilt2(r, [k k], 'symmetric');
processUsingLevelB = (zmed > zmin) & (zmax > zmed) & ...
~alreadyProcessed;
zB = (r > zmin) & (zmax > r);
outputZxy = processUsingLevelB & zB;
outputZmed = processUsingLevelB & ~zB;
f(outputZxy) = r(outputZxy);
f(outputZmed) = zmed(outputZmed);
alreadyProcessed = alreadyProcessed | processUsingLevelB;
if all(alreadyProcessed(:))
break;
end
end
f(~alreadyProcessed) = zmed(~alreadyProcessed);