SeparateKernel.m

function [K1 KN ERR]=SeparateKernel(H)
% This function SEPARATEKERNEL will separate ( do decomposition ) any 
% 2D, 3D or nD kernel into 1D kernels. Ofcourse only a sub-set of Kernels
% are separable such as a Gaussian Kernel, but it will give least-squares
% sollutions for non-separatable kernels.
% 
% Separating a 3D or 4D image filter in to 1D filters will give an large
% speed-up in image filtering with for instance the function imfilter.
%
% [K1 KN ERR]=SeparateKernel(H);
%   
% inputs,
%   H : The 2D, 3D ..., ND kernel
%   
% outputs,
%   K1 : Cell array with the 1D kernels
%   KN : Approximation of the ND input kernel by the 1D kernels
%   ERR : The sum of absolute difference between approximation and input kernel
%
% 
% How the algorithm works:
% If we have a separable kernel like
% 
%  H = [1 2 1
%       2 4 2
%       3 6 3];
%
% We like to solve unknow 1D kernels,
%  a=[a(1) a(2) a(3)]
%  b=[b(1) b(2) b(3)]
%
% We know that,
%  H = a'*b
%
%       b(1)    b(2)    b(3)
%       --------------------
%  a(1)|h(1,1) h(1,2) h(1,3)
%  a(2)|h(2,1) h(2,2) h(2,3)
%  a(3)|h(3,1) h(3,2) h(3,3)
%
% Thus,
%  h(1,1) == a(1)*b(1)
%  h(2,1) == a(2)*b(1)
%  h(3,1) == a(3)*b(1)
%  h(4,1) == a(1)*b(2)
%  ...
%
% We want to solve this by using fast matrix (least squares) math,
%
%  c = M * d; 
%  
%  c a column vector with all kernel values H
%  d a column vector with the unknown 1D kernels 
%
% But matrices "add" values and we have something like  h(1,1) == a(1)*b(1);
% We solve this by taking the log at both sides 
% (We replace zeros by a small value. Whole lines/planes of zeros are
%  removed at forehand and re-added afterwards)
%
%  log( h(1,1) ) == log(a(1)) + log b(1))
%
% The matrix is something like this,
%
%      a1 a2 a3 b1 b2 b3    
% M = [1  0  0  1  0  0;  h11
%      0  1  0  1  0  0;  h21
%      0  0  1  1  0  0;  h31
%      1  0  0  0  1  0;  h21
%      0  1  0  0  1  0;  h22
%      0  0  1  0  1  0;  h23
%      1  0  0  0  0  1;  h31
%      0  1  0  0  0  1;  h32
%      0  0  1  0  0  1]; h33
%
% Least squares solution
%  d = exp(M\log(c))
%
% with the 1D kernels
%
%  [a(1);a(2);a(3);b(1);b(2);b(3)] = d
%
% The Problem of Negative Values!!!
%
% The log of a negative value is possible it gives a complex value, log(-1) = i*pi
% if we take the expontential it is back to the old value, exp(i*pi) = -1 
%
%  But if we use the solver with on of the 1D vectors we get something like, this :
%
%  input         result        abs(result)    angle(result) 
%   -1     -0.0026 + 0.0125i     0.0128         1.7744 
%    2      0.0117 + 0.0228i     0.0256         1.0958 
%   -3     -0.0078 + 0.0376i     0.0384         1.7744  
%    4      0.0234 + 0.0455i     0.0512         1.0958
%    5      0.0293 + 0.0569i     0.0640         1.0958
% 
% The absolute value is indeed correct (difference in scale is compensated
% by the order 1D vectors)
%
% As you can see the angle is correlated with the sign of the values. But I
% didn't found the correlation yet. For some matrices it is something like
%
%  sign=mod(angle(solution)*scale,pi) == pi/2;
%
% In the current algorithm, we just flip the 1D kernel values one by one.
% The sign change which gives the smallest error is permanently swapped. 
% Until swapping signs no longer decreases the error
%
% Examples,
%   a=permute(rand(5,1),[1 2 3 4])-0.5;
%   b=permute(rand(5,1),[2 1 3 4])-0.5;
%   c=permute(rand(5,1),[3 2 1 4])-0.5;
%   d=permute(rand(5,1),[4 2 3 1])-0.5;
%   H = repmat(a,[1 5 5 5]).*repmat(b,[5 1 5 5]).*repmat(c,[5 5 1 5]).*repmat(d,[5 5 5 1]);
%   [K,KN,err]=SeparateKernel(H);
%   disp(['Summed Absolute Error between Real and approximation by 1D filters : ' num2str(err)]);
%
%   a=permute(rand(3,1),[1 2 3])-0.5;
%   b=permute(rand(3,1),[2 1 3])-0.5;
%   c=permute(rand(3,1),[3 2 1])-0.5;
%   H = repmat(a,[1 3 3]).*repmat(b,[3 1 3 ]).*repmat(c,[3 3 1 ])
%   [K,KN,err]=SeparateKernel(H); err
%
%   a=permute(rand(4,1),[1 2 3])-0.5;
%   b=permute(rand(4,1),[2 1 3])-0.5;
%   H = repmat(a,[1 4]).*repmat(b,[4 1]);
%   [K,KN,err]=SeparateKernel(H); err
%
% Function is written by D.Kroon, uses "log idea" from A. J. Hendrikse, 
% University of Twente (July 2010)


% We first make some structure which contains information about
% the transformation from kernel to 1D kernel array, number of dimensions
% and other stuff
data=InitializeDataStruct(H);

% Make the matrix of c = M * d; 
M=makeMatrix(data);

% Solve c = M * d with least squares
warning('off','MATLAB:rankDeficientMatrix');
par=exp(M\log(abs(data.H(:))));

% Improve the values by solving the remaining difference
KN = Filter1DtoFilterND(par,data);
par2=exp(M\log(abs(KN(:)./data.H(:))));
par=par./par2;

% Change the sign of a 1D filtering value if it decrease the error
par = FilterCorrSign(par,data);

% Split the solution d in separate 1D kernels
K1 = ValueList2Filter1D(par,data);

% Re-add the removed zero rows/planes to the 1D vectors
K1=re_add_zero_rows(data,K1);

% Calculate the approximation of the ND kernel if using the 1D kernels
KN = Filter1DtoFilterND(par,data,K1);

% Calculate the absolute error
ERR =sum(abs(H(:)-KN(:)));

function par = FilterCorrSign(par,data)
Ert=zeros(1,length(par));
ERR=inf; t=0;
par=sign(rand(size(par))-0.5).*par;
while(t<ERR)
    % Calculate the approximation of the ND kernel if using the 1D kernels
    KN = Filter1DtoFilterND(par,data);
    % Calculate the absolute error
    ERR =sum(abs(data.H(:)-KN(:)));
    % Flip the sign of every 1D filter value, and look if the error
    % improves
    for i=1:length(par)
        par2=par; par2(i)=-par2(i);
        KN = Filter1DtoFilterND(par2,data);
        Ert(i) =sum(abs(data.H(:)-KN(:)));
    end
    % Flip the sign of the 1D filter value with the largest improvement
    [t,j]=min(Ert); if(t<ERR), par(j)=-par(j); end
end

function data=InitializeDataStruct(H)
data.sizeHreal=size(H);
data.nreal=ndims(H);
[H,preserve_zeros]=remove_zero_rows(H);
data.H=H;
data.n=ndims(H);
data.preserve_zeros=preserve_zeros;
data.H(H==0)=eps;
data.sizeH=size(data.H);
data.sep_parb=cumsum([1 data.sizeH(1:data.n-1)]);
data.sep_pare=cumsum(data.sizeH);
data.sep_parl=data.sep_pare-data.sep_parb+1;
data.par=(1:numel(H))+1;

function [H,preserve_zeros]=remove_zero_rows(H)
% Remove whole columns/rows/planes with zeros,
% because we know at forehand that they will give a kernel 1D value of 0
% and will otherwise increase the error in the end result.
preserve_zeros=zeros(numel(H),2); pz=0;
sizeH=size(H);
for i=1:ndims(H)
    H2D=reshape(H,size(H,1),[]);
    check_zero=~any(H2D,2);
    if(any(check_zero))
        zero_rows=find(check_zero);
        for j=1:length(zero_rows)
            pz=pz+1;
            preserve_zeros(pz,:)=[i zero_rows(j)];
            sizeH(1)=sizeH(1)-1;
        end
        H2D(check_zero,:)=[];
        H=reshape(H2D,sizeH);
    end
    H=shiftdim(H,1);
    sizeH=circshift(sizeH,[0 -1]);
    H=reshape(H,sizeH);
end
preserve_zeros=preserve_zeros(1:pz,:);

function K1=re_add_zero_rows(data,K1)
% Re-add the 1D kernel values responding to a whole column/row or plane
% of zeros
for i=1:size(data.preserve_zeros,1)
    di=data.preserve_zeros(i,1);
    pos=data.preserve_zeros(i,2);
    if(di>length(K1)), K1{di}=1; end
    val=K1{di};
    val=val(:);
    val=[val(1:pos-1);0;val(pos:end)];
    dim=ones(1,data.nreal); dim(di)=length(val);
    K1{di}=reshape(val,dim);
end

function M=makeMatrix(data)
 M = zeros(numel(data.H),sum(data.sizeH));
 K1 = (1:numel(data.H))';
 for i=1:data.n;
    p=data.par(data.sep_parb(i):data.sep_pare(i)); p=p(:);
    dim=ones(1,data.n); dim(i)=data.sep_parl(i);
    Ki=reshape(p(:),dim);
    dim=data.sizeH; dim(i)=1;
    K2=repmat(Ki,dim)-1;
    M(sub2ind(size(M),K1(:),K2(:)))=1;
 end
 
function Kt = Filter1DtoFilterND(par,data,K1)
if(nargin==2)
 Kt=ones(data.sizeH);
 for i=1:data.n
     p=par(data.sep_parb(i):data.sep_pare(i)); p=p(:);
     dim=ones(1,data.n); dim(i)=data.sep_parl(i);
     Ki=reshape(p(:),dim);
     dim=data.sizeH; dim(i)=1;
     Kt=Kt.*repmat(Ki,dim);
 end
else
  Kt=ones(data.sizeHreal);
  for i=1:data.n
    dim=data.sizeHreal; dim(i)=1;
    Kt=Kt.*repmat(K1{i},dim);
  end
end

function K = ValueList2Filter1D(par,data)
 K=cell(1,data.n);
 for i=1:data.n
     p=par(data.sep_parb(i):data.sep_pare(i)); p=p(:);
     dim=ones(1,data.n); dim(i)=data.sep_parl(i);
     K{i}=reshape(p(:),dim);
 end