Code covered by the BSD License  

Highlights from
Separate Kernel in 1D kernels

5.0

5.0 | 2 ratings Rate this file 30 Downloads (last 30 days) File Size: 3.99 KB File ID: #28218
image thumbnail

Separate Kernel in 1D kernels

by

 

16 Jul 2010 (Updated )

Decompose an arbitrary N dimensional filtering kernel into 1D kernels, for faster filtering

| Watch this File

File Information
Description

 This function SEPARATEKERNEL will separate ( do decomposition of ) any
  2D, 3D or nD kernel into 1D kernels. Of course only a sub-set of Kernels
  are separable such as a Gaussian Kernel, but it will give approximations for non-separable kernels.
  
  Separating a 3D or 5D image filter into 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 unknown 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 other 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

.
---------------------------------------------------------------------
Alternative SVD Method
---------------------------------------------------------------------
Cris Luengo made the same function but then based on the SVD of 2D kernels. His algorithm is faster with large kernels (because SVD is an build-in function), an deals more elegant with negative and zero values.
http://www.mathworks.com/matlabcentral/fileexchange/28238-kernel-decomposition

.
.
---------------------------------------------------------------------
Question
---------------------------------------------------------------------
If someone knows the mathematical / analytical solution for finding the signs of the 1D filter kernels, please let me know.

Acknowledgements

This file inspired Kernel Decomposition.

MATLAB release MATLAB 7.10 (R2010a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (2)
06 Sep 2010 dai zhengguo

Hi Dirk-Jan. Your code is excellent to be used. However, there is a small bug here, if you use the test code below:
for i = 1:20
v = [1; 2; 1] ;
h = [-1 0 1] ;
s = v*h ;
[K,KN,err] = SeparateKernel( s ) ;
err
end
Please check it.

12 Aug 2010 Cris Luengo

Hi Dirk-Jan. That's an elegant solution to the zero-value problem.

Updates
16 Jul 2010

Changed text on website

27 Jul 2010

Solved zero value problems

Contact us