algstdfilt.m

%This is what stdfilt actually runs after checking the input - Y.B.

function J = algstdfilt(I,h)
% Main algorithm used by stdfilt function. See stdfilt for more
% details

% No input validation is done in this function.


% Copyright 2013-2015 The MathWorks, Inc.


n = sum(h(:));

% If n = 1 then return default J (all zeros).
% Otherwise, calculate standard deviation. The formula for standard deviation
% can be rewritten in terms of the theoretical definition of
% convolution. However, in practise, use correlation in IMFILTER to avoid a
% flipped answer when NHOOD is asymmetric.
% conv1 = imfilter(I.^2,h,'symmetric') / (n-1); 
% conv2 = imfilter(I,h,'symmetric').^2 / (n*(n-1));
% std = sqrt(conv1-conv2).  
% These equations can be further optimized for speed.

n1 = n - 1;
if n ~= 1
  conv1 = imfilter(I.^2, h/n1 , 'symmetric');
  conv2 = imfilter(I, h/sqrt(n*n1), 'symmetric').^2;
  J = sqrt(max((conv1 - conv2),0));
else
  J = zeros(size(I));
end

error

$ nosetests -sx test_calculateSC.py:test_std_filt_m
F
======================================================================
FAIL: test_calculateSC.test_std_filt_m
----------------------------------------------------------------------
Traceback (most recent call last):
  File "C:\Users\mrbli_000\Desktop\Python\Anaconda\lib\site-packages\nose\case.py", line 198, in runTest
    self.test(*self.arg)
  File "C:\Users\mrbli_000\Desktop\f0cl\edenis-confidential\testrepo\edenis-v4msti\tests\calculateSC\test_calculateSC.py", line 70, in test_std_filt_m
    return std_filt_m(test_file,std_filter)
  File "<decorator-gen-2>", line 2, in std_filt_m
  File "C:\Users\mrbli_000\Desktop\f0cl\edenis-confidential\testrepo\edenis-v4msti\tests\calculateSC\rt_insert.py", line 8, in _wrapper
    return self.callback(wrapped, _args, _dargs, args, dargs)
  File "C:\Users\mrbli_000\Desktop\f0cl\edenis-confidential\testrepo\edenis-v4msti\tests\calculateSC\test_calculateSC.py", line 38, in callback
    assert np.allclose(mfunc_result, local_result), mfunc_name
AssertionError: stdfilt

----------------------------------------------------------------------
Ran 1 test in 0.158s

FAILED (failures=1)
stdfilt
octave:
[[ 1.67332005  2.          2.16794834  1.81659021  1.        ]
 [ 1.67332005  2.          2.16794834  1.81659021  1.        ]
 [ 1.67332005  2.          2.16794834  1.81659021  1.        ]
 [ 1.67332005  2.          2.16794834  1.81659021  1.        ]
 [ 1.67332005  2.          2.16794834  1.81659021  1.        ]]
python:
[[ 0.83666003  1.30384048  1.58113883  1.30384048  0.83666003]
 [ 0.83666003  1.30384048  1.58113883  1.30384048  0.83666003]
 [ 0.83666003  1.30384048  1.58113883  1.30384048  0.83666003]
 [ 0.83666003  1.30384048  1.58113883  1.30384048  0.83666003]
 [ 0.83666003  1.30384048  1.58113883  1.30384048  0.83666003]]
octave/python:
[[ 2.          1.53392998  1.37113092  1.39326109  1.19522861]
 [ 2.          1.53392998  1.37113092  1.39326109  1.19522861]
 [ 2.          1.53392998  1.37113092  1.39326109  1.19522861]
 [ 2.          1.53392998  1.37113092  1.39326109  1.19522861]
 [ 2.          1.53392998  1.37113092  1.39326109  1.19522861]]

std_filt_m.py

@mimic('stdfilt')
def std_filt_m(img, kernel):
    img = np.transpose(img)
    n = kernel.sum()
    n1 = n - 1
    c1 = filter2D(img**2, -1, kernel/n1, borderType = BORDER_REFLECT)
    c2 = filter2D(img, -1, kernel/np.sqrt(n*n1), borderType = BORDER_REFLECT)**2
    J = np.sqrt(np.maximum(c1-c2, 0))
    return np.transpose(J)

test_std_filt_m.py

def test_std_filt_m():
    window_size = 5
    std_filter = np.ones(window_size)
    test_file = (np.arange(25.)).reshape(5,5)
    return std_filt_m(test_file,std_filter)