source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/qtfm/householder_matrix.m @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d

File size: 1.8 KB
Line 
1function [h] = householder_matrix(a, v)
2% h = householder_matrix(a, v) returns the Householder matrix that
3% will zero all elements of a except those corresponding to (any)
4% non-zero elements of v. If v has zero norm, an identity matrix is
5% returned. The matrix returned is either a left or right Householder
6% matrix dependent on whether the input parameters are column or row
7% vectors respectively.
8
9% Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan.
10% See the file : Copyright.m for further details.
11
12% Reference:
13%
14% Sangwine, S. J. and Le Bihan, N.,
15% Quaternion singular value decomposition based on bidiagonalization
16% to a real or complex matrix using quaternion Householder transformations,
17% Applied Mathematics and Computation, 182(1), 1 November 2006, 727-738,
18% DOI:10.1016/j.amc.2006.04.032.
19%
20% S. J. Sangwine and N. Le Bihan,
21% Quaternion Singular Value Decomposition based on Bidiagonalization
22% to a Real Matrix using Quaternion Householder Transformations,
23% arXiv:math.NA/0603251, 10 March 2006. Available at http://www.arxiv.org/.
24
25error(nargchk(2, 2, nargin)), error(nargoutchk(1, 1, nargout))
26
27if size(a) ~= size(v)
28    error('Input parameters must be vectors of the same size, both row, or both column')
29end
30
31[m, n] = size(a);
32
33col_vector = n == 1 & m == length(a);
34row_vector = m == 1 & n == length(a);
35
36if col_vector == 0 & row_vector == 0
37    error('Parameters must be either column or row vectors.');   
38end
39
40h = quaternion(eye(length(a)));
41
42n = norm(v); % If v has zero norm, we return an identity matrix.
43
44if n ~= 0
45    [u, zeta] = householder_vector(a, v ./ n);
46    if size(u) ~= size(a)
47        error('Result from function householder_vector is not of correct type (row, column)');
48    end
49    if col_vector
50        h = (1 ./ zeta) .* (h - u * u');
51    else
52        h = (h - u.' * conj(u)) .* (1 ./ zeta);
53    end
54end
Note: See TracBrowser for help on using the repository browser.