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 |
---|