function [ gaussian ] = gaussian( x, sigma, order, normalize ) if isempty(normalize) normalize = false; end gaussian_base = exp(-(x.*x)/(2*sigma^2)); % we will do the same as ITK. this is just done to be able to % compare gaussian results with those coming from ITK. if (normalize) gaussian_base = sigma*gaussian_base; end if order == 0 gaussian = 1/(sigma*sqrt(2*pi)) .* gaussian_base; elseif order == 1 gaussian = -x./(sigma^2*sqrt(2*pi)) .* gaussian_base; elseif order == 2 gaussian = - gaussian_base ./ (sigma^3*sqrt(2*pi)) + x.*x.*gaussian_base / (sigma^5*sqrt(2*pi)); end end