{"id":98,"date":"2006-11-28T11:23:47","date_gmt":"2006-11-28T16:23:47","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=98"},"modified":"2019-10-22T16:35:53","modified_gmt":"2019-10-22T20:35:53","slug":"separable-convolution-part-2","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2006\/11\/28\/separable-convolution-part-2\/","title":{"rendered":"Separable convolution: Part 2"},"content":{"rendered":"
\r\n

Back in October<\/a> I introduced the concept of filter separability<\/i>. A two-dimensional filter s<\/i> is said to be separable if it can be written as the convolution of two one-dimensional filters v<\/i> and h<\/i>:\r\n <\/p>\r\n

<\/p>\r\n

I said then that \"next time\" I would explain how to determine whether a given filter is separable. Well, I guess I got side-tracked,\r\n but I'm back on topic now.\r\n <\/p>\r\n

This question gave me one of earliest opportunities at The MathWorks to wander down to company co-founder Cleve's<\/a> office and ask for advice. I asked, \"How can I determine if a matrix is an outer product of two vectors?\" Cleve was very\r\n helpful, as he always is, although I was a little embarrassed afterward that I hadn't figured it out myself. \"Go look at\r\n the rank<\/tt><\/a> function,\" Cleve told me.\r\n <\/p>\r\n

Of course. If a matrix is an outer product of two vectors, its rank is 1. Here are the key lines of code in rank<\/tt>:\r\n <\/p>

dbtype 15:19<\/span> rank<\/span><\/pre>
\r\n15    s = svd(A);\r\n16    if nargin==1\r\n17       tol = max(size(A)') * eps(max(s));\r\n18    end\r\n19    r = sum(s > tol);\r\n\r\n<\/pre>
So the test is this:  The rank of A<\/tt> is the number of nonzero singular values of A<\/tt>, with some numerical tolerance based on eps<\/tt> and the size of A<\/tt>.\r\n   <\/p>\r\n   
Let's try it with a few common filters.<\/p>\r\n   
An averaging filter should be obvious:<\/p>
averaging = ones(5,5) \/ 25;\r\nrank(averaging)<\/pre>
\r\nans =\r\n\r\n     1\r\n\r\n<\/pre>
The Sobel kernel:<\/p>
sobel = [-1 0 1; -2 0 2; -1 0 1];\r\nrank(sobel)<\/pre>
\r\nans =\r\n\r\n     1\r\n\r\n<\/pre>
The two-dimensional Gaussian is the only radially symmetric function that is also separable:<\/p>
gaussian = fspecial('gaussian'<\/span>);\r\nrank(gaussian)<\/pre>
\r\nans =\r\n\r\n     1\r\n\r\n<\/pre>
A disk is not separable:<\/p>
disk = fspecial('disk'<\/span>);\r\nrank(disk)<\/pre>
\r\nans =\r\n\r\n     5\r\n\r\n<\/pre>
So how can we determine the outer product vectors? The answer is to go back to the svd<\/tt> function.  Here's a snippet from the doc:  [U,S,V] = svd(X)<\/tt> produces a diagonal matrix S<\/tt> of the same dimension as X<\/tt>, with nonnegative diagonal elements in decreasing order, and unitary matrices U<\/tt> and V<\/tt> so that X = U*S*V'<\/tt>.\r\n   <\/p>\r\n   
A rank 1 matrix has only one nonzero singular value, so U*S*V'<\/tt> becomes U(:,1) * S(1,1) * V(:,1)'<\/tt>.  This is basically the outer product we were seeking.  Therefore, we want the first columns of U<\/tt> and V<\/tt>. (We have to remember also to use the nonzero singular value as a scale factor.)\r\n   <\/p>\r\n   
Let's try this with the Gaussian filter.<\/p>
[U,S,V] = svd(gaussian)<\/pre>
\r\nU =\r\n\r\n   -0.1329    0.9581   -0.2537\r\n   -0.9822   -0.1617   -0.0959\r\n   -0.1329    0.2364    0.9625\r\n\r\n\r\nS =\r\n\r\n    0.6420         0         0\r\n         0    0.0000         0\r\n         0         0    0.0000\r\n\r\n\r\nV =\r\n\r\n   -0.1329   -0.6945   -0.7071\r\n   -0.9822    0.1880    0.0000\r\n   -0.1329   -0.6945    0.7071\r\n\r\n<\/pre>
Now get the horizontal and vertical vectors from the first columns of U<\/tt> and V<\/tt>.\r\n   <\/p>
v = U(:,1) * sqrt(S(1,1))<\/pre>
\r\nv =\r\n\r\n   -0.1065\r\n   -0.7870\r\n   -0.1065\r\n\r\n<\/pre>
h = V(:,1)' * sqrt(S(1,1))<\/pre>
\r\nh =\r\n\r\n   -0.1065   -0.7870   -0.1065\r\n\r\n<\/pre>
I have chosen, somewhat arbitrarily to split the scale factor, S(1,1)<\/tt>, \"equally\" between v<\/tt> and h<\/tt>.\r\n   <\/p>\r\n   
Now check to make sure this works:<\/p>
gaussian - v*h<\/pre>
\r\nans =\r\n\r\n  1.0e-015 *\r\n\r\n   -0.0243   -0.1527   -0.0243\r\n   -0.0139   -0.1110   -0.0139\r\n   -0.0035         0   -0.0035\r\n\r\n<\/pre>
Except for normal floating-point roundoff differences, gaussian<\/tt> and v*h<\/tt> are equal.\r\n   <\/p>\r\n   
You can find code similar to this in the MATLAB function filter2<\/tt><\/a>, as well as in the Image Processing Toolbox function imfilter<\/tt>.\r\n   <\/p>