[37] | 1 | function C = conv2(w, x, y, z) |
---|
| 2 | % CONV2 Two dimensional convolution. |
---|
| 3 | % (Quaternion overloading of standard Matlab function.) |
---|
| 4 | |
---|
| 5 | % Copyright © 2006 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
| 6 | % See the file : Copyright.m for further details. |
---|
| 7 | |
---|
| 8 | % This function operates in a similar way to the standard Matlab function |
---|
| 9 | % apart from supporting both left and right coefficients. (Since quaternion |
---|
| 10 | % multiplication is not commutative, the general case requires both left |
---|
| 11 | % and right multiplication in the convolution product/summation.) The |
---|
| 12 | % Matlab function allows the first two parameters to be vectors - this is |
---|
| 13 | % not implemented as yet. Acceptable calling profiles are: |
---|
| 14 | % |
---|
| 15 | % C = conv2(A, B) - A is convolved on the left of B, that is A * B |
---|
| 16 | % C = conv2({L, R}, B) - The convolution is L * B * R. L and R must be of |
---|
| 17 | % the same size. |
---|
| 18 | % |
---|
| 19 | % An optional last parameter can specify 'shape' as for the standard Matlab |
---|
| 20 | % function. This is currently not implemented. |
---|
| 21 | |
---|
| 22 | error(nargchk(2, 4, nargin)), error(nargoutchk(0, 1, nargout)) |
---|
| 23 | |
---|
| 24 | if nargin == 4 |
---|
| 25 | error('conv2 does not yet support 4 parameters.') |
---|
| 26 | end |
---|
| 27 | |
---|
| 28 | % If two or three parameters are supplied, the first two can be either two |
---|
| 29 | % matrices, or a cell array and a matrix. The cell array must contain two |
---|
| 30 | % matrices of the same size, but they need not be quaternion valued, since |
---|
| 31 | % they could be real or complex. |
---|
| 32 | |
---|
| 33 | if iscell(w) |
---|
| 34 | % The first parameter must be a cell array of two elements, and the |
---|
| 35 | % elements must have the same size. |
---|
| 36 | |
---|
| 37 | if length(w) ~= 2 |
---|
| 38 | error('If a cell array, the first parameter must have two elements.') |
---|
| 39 | end |
---|
| 40 | |
---|
| 41 | L = w{1}; R = w{2}; |
---|
| 42 | |
---|
| 43 | if ~all(size(L) == size(R)) |
---|
| 44 | error('The elements of the cell array must have the same size.') |
---|
| 45 | end |
---|
| 46 | |
---|
| 47 | B = x; % The second parameter must be a matrix. |
---|
| 48 | % Should we do any checks on it here? |
---|
| 49 | else |
---|
| 50 | % The first parameter must be a matrix (or vector). It need not |
---|
| 51 | % be quaternion-valued, since it could be a left matrix of |
---|
| 52 | % coefficients. We treat the smaller matrix as a coefficient |
---|
| 53 | % matrix, but put it on the left or right according to its |
---|
| 54 | % position in the parameter list, and create a matrix of ones |
---|
| 55 | % of the same size for the coefficient on the opposite side. |
---|
| 56 | |
---|
| 57 | if numel(w) < numel(x) % Is this the best comparison that could be done? |
---|
| 58 | L = w; B = x; R = ones(size(w)); |
---|
| 59 | else |
---|
| 60 | L = ones(size(x)); B = w; R = x; |
---|
| 61 | end |
---|
| 62 | end |
---|
| 63 | |
---|
| 64 | if nargin == 3 |
---|
| 65 | |
---|
| 66 | % The third parameter must be a 'shape' parameter. Check for the |
---|
| 67 | % legal possibilities. |
---|
| 68 | |
---|
| 69 | if ~ischar(y) |
---|
| 70 | error('Third parameter must be a string.'); |
---|
| 71 | end |
---|
| 72 | |
---|
| 73 | if ~strcmp(y, 'full') && ~strcmp(y, 'same') && ~strcmp(y, 'valid') |
---|
| 74 | error(['The third parameter must be one of', ... |
---|
| 75 | '''full'', ''same'', or ''valid''']) |
---|
| 76 | else |
---|
| 77 | shape = y; |
---|
| 78 | end |
---|
| 79 | else |
---|
| 80 | shape = 'full'; % This wasn't specified, so use a default value. |
---|
| 81 | end |
---|
| 82 | |
---|
| 83 | % The variables L, B, R now contain the matrices to be convolved, in the |
---|
| 84 | % order L * B * R. |
---|
| 85 | |
---|
| 86 | [ma, na] = size(L); % This is the same as size(R), of course (we checked that above). |
---|
| 87 | [mb, nb] = size(B); |
---|
| 88 | mc = ma + mb - 1; |
---|
| 89 | nc = na + nb - 1; |
---|
| 90 | |
---|
| 91 | C = quaternion(zeros(mc, nc)); % Preconstruct the result array. |
---|
| 92 | |
---|
| 93 | for n1 = 1 : mc |
---|
| 94 | for n2 = 1 : nc |
---|
| 95 | |
---|
| 96 | % This code is vectorized: the two arrays are sliced to extract the |
---|
| 97 | % parts that must be multiplied element by element and summed to give |
---|
| 98 | % the current element of the output array. This is complicated by the |
---|
| 99 | % cases where one array overlaps the edge(s) of the other. |
---|
| 100 | |
---|
| 101 | k11 = max(1, n1 + 1 - mb); k12 = min(n1, ma); |
---|
| 102 | k21 = max(1, n2 + 1 - nb); k22 = min(n2, na); |
---|
| 103 | |
---|
| 104 | % We need to use indexing here, but since this is a class method, |
---|
| 105 | % normal array indexing doesn't work. See the comment in the file |
---|
| 106 | % svd.m where the same problem and solution occurs. The commented |
---|
| 107 | % statement below is what we would write if array indexing worked |
---|
| 108 | % here as per normal. |
---|
| 109 | % C(n1, n2) = ... |
---|
| 110 | % sum(L(k11:k12, k21:k22) * B(n1 + 1 - k1, n2 + 1 - k2) * ... |
---|
| 111 | % R(k11:k12, k21:k22)); |
---|
| 112 | |
---|
| 113 | s = substruct('()',{n1, n2}); |
---|
| 114 | t = substruct('()',{k11:k12, k21:k22}); |
---|
| 115 | C = subsasgn(C, s, ... |
---|
| 116 | sum(sum(subsref(L, t) .* ... |
---|
| 117 | flipud(fliplr(subsref(B, ... |
---|
| 118 | substruct('()', {n1 + 1 - k12 : n1 + 1 - k11, ... |
---|
| 119 | n2 + 1 - k22 : n2 + 1 - k21})))) .* ... |
---|
| 120 | subsref(R, t)) ... |
---|
| 121 | )); |
---|
| 122 | end |
---|
| 123 | end |
---|
| 124 | |
---|
| 125 | if ~strcmp(shape, 'full') |
---|
| 126 | error('The shape parameter is not implemented.') |
---|
| 127 | end |
---|