edge4d(A)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function S = edge4d(A)
%
% 4-D edge Sobel edge detector
%
% INPUT:
%   A       - 4-d double/single tensor
%
% OUTPUT:
%   S       - 4-d double/single tensor
%               Sobel edge magnitude values at every location (x,y,z,t)
%
% Class support of input:
%     float: double, single
% 
% 
% Copyright 2011-     by Rene Fung
%
% Permission is granted for anyone to copy, use, or modify this
% software and accompanying documents for any uncommercial
% purposes, provided this copyright notice is retained, and note is
% made of any changes that have been made. This software and
% documents are distributed without any warranty, express or
% implied
%
% Last Modified: 09/30/2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function S = edge4d(A)

A = double(A);

% 1D Sobel operators
h = [1 2 1];
hd = [1 0 -1];

% 4D Sobel kernels in x,y,z,t directions
for x = 1 : 3
    for y = 1 : 3
        for z = 1 : 3
            for t = 1 : 3
                h4x(x,y,z,t) = hd(x) * h(y) * h(z) * h(t);
                h4y(x,y,z,t) = h(x) * hd(y) * h(z) * h(t);
                h4z(x,y,z,t) = h(x) * h(y) * hd(z) * h(t);
                h4t(x,y,z,t) = h(x) * h(y) * h(z) * hd(t);
            end
        end
    end
end

% Convolute the n-d data with Sobel kernels
Sx = imfilter(A,h4x);
Sy = imfilter(A,h4y);
Sz = imfilter(A,h4z);
St = imfilter(A,h4t);

% Normalize gradient magnitude from Sobel operator
S = Sx.^2 + Sy.^2 + Sz.^2 + St.^2;
S = S ./ max(S(:));