Appendices - Feature Extraction and Image Processing - page 335

direct and Gaussian averaging, the salt and pepper contributes to the final result. To run it
again, with a different selection of noise points, just select the function noisy_p
:=addcondiments() and run it again. Each time, you'll get a different pattern of
noise. Check the filtering still works.
There is of course a median operator in Mathcad, but we thought we'd show you how
it worked. Median filters can be implemented in an adaptive manner, or using non-square
templates (e.g. cross or line, usually for computational reasons). We can get a Mathcad
median by:

their_med(pic,winsize):= newpic ← zero(pic)

winsize
2

half floor

for x ∈ half..cols(pic)-half-1
for y ∈ half..rows(pic)-half-1
newpic y,x ← median(submatrix
(pic,y-half,y+half,x-half,x+half))

newpic

This gives you a median operator for an arbitrary template size.
Finally, the last statistic is the mode . This is the peak of the probability distribution (the
value most likely to occur). One way to estimate its value is to use the truncated median
filter . It operates by taking the median of the distribution resulting from truncation of the
distribution within a window at a point beyond the mean.
Let's have a picture of an artery to work on:

noisy_p:=READBMP(artery)
Now, here's the code:

trun_med(p,wsze):= newpic ← zero(p)

wsze
2

ha floor

for x ∈ ha..cols(p)-ha-1
for y ∈ ha..rows(p)-ha-1
win ← submatrix(p,y-ha,y+ha,x-ha,x+ha)
med ← median(win)
ave ← mean(win)
upper ←

2·med-min(win)

lower ← 2·med-max(win)

cc ← 0
for i ∈ 0..wsze-1
for j ∈ 0..wsze-1
trun cc ← win j,i if(win j,i <upper)·(med<ave)
trun cc ← win j,i if(win j,i >lower)·(med>ave)
cc ← cc+1
newpic y,x ← median(trun)

newpic