Description of vectorize

0001 function [I,D] = vectorize_inverse(M);
0002 % vectorize_inverse: fast computation of 3x3 and 4x4 inverses
0003 % [Inv,Det] = vectorize_inverse(Mat);
0004 % Mat is a 3x3xN or 4x4xN matrix
0005 % Inv is inverse where Inv(:,:,k) = inv(Mat(:,:,k))
0006 % Det is vector of determinants
0007  
0008 % (C) 2022 A Adler. License GPL v2 or v3
0009 
0010 if ischar(M) && strcmp(M,'UNIT_TEST'); do_unit_test; return; end
0011 
0012 szM = size(M);
0013 if szM(1) ~= szM(2)
0014     error('vectorize_inverse: matrix must be square');
0015 end
0016 switch szM(1)
0017    case 3; [I,D] = vectorize_3x3inv(M);
0018    case 4; [I,D] = vectorize_4x4inv(M);
0019    otherwise;
0020       error('Only 3x3 and 4x4 matrices possible')
0021 end
0022 
0023 function [I,D] = vectorize_3x3inv(M)
0024    a= M(1,1,:); b= M(1,2,:); c= M(1,3,:);
0025    d= M(2,1,:); e= M(2,2,:); f= M(2,3,:);
0026    g= M(3,1,:); h= M(3,2,:); i= M(3,3,:);
0027 
0028    D= a.*(e.*i-f.*h) - b.*(d.*i-f.*g) + c.*(d.*h-e.*g);
0029    if any(abs(D) < eps) 
0030       warning('Determinant close to zero');
0031    end
0032 
0033    I= (1./D).* ...
0034       [e.*i - f.*h, c.*h - b.*i, b.*f - c.*e;
0035        f.*g - d.*i, a.*i - c.*g, c.*d - a.*f;
0036        d.*h - e.*g, b.*g - a.*h, a.*e - b.*d];
0037    
0038 
0039 function [I,D] = vectorize_4x4inv(M)
0040 % Adapted from the Mesa3D implementation of
0041 % gluInvertMatrix(const double m[16], double invOut[16])
0042 % Available in Mesa3D (License in MIT)
0043 %
0044 % positions of 1 removed in rev6202
0045 
0046 
0047 % Precalculate pieces to speed up
0048    M21 = M(2,1,:); M31 = M(3,1,:); M41 = M(4,1,:);
0049    M22 = M(2,2,:); M32 = M(3,2,:); M42 = M(4,2,:);
0050    M23 = M(2,3,:); M33 = M(3,3,:); M43 = M(4,3,:);
0051    M24 = M(2,4,:); M34 = M(3,4,:); M44 = M(4,4,:);
0052 
0053    I11= M22.*M33.*M44 - ...
0054              M22.*M43.*M34 - ...
0055              M23.*M32.*M44 + ...
0056              M23.*M42.*M34 + ...
0057              M24.*M32.*M43 - ...
0058              M24.*M42.*M33;
0059 
0060    I12=-M33.*M44 + ...
0061              M43.*M34 + ...
0062              M32.*M44 - ...
0063              M42.*M34 - ...
0064              M32.*M43 + ...
0065              M42.*M33;
0066 
0067    I13= M23.*M44 - ...
0068              M43.*M24 - ...
0069              M22.*M44 + ...
0070              M42.*M24 + ...
0071              M22.*M43 - ...
0072              M42.*M23;
0073 
0074    I14=-M23.*M34 + ...
0075              M33.*M24 + ...
0076              M22.*M34 - ...
0077              M32.*M24 - ...
0078              M22.*M33 + ...
0079              M32.*M23;
0080 
0081    I21=-M21.*M33.*M44 + ...
0082              M21.*M43.*M34 + ...
0083              M23.*M31.*M44 - ...
0084              M23.*M41.*M34 - ...
0085              M24.*M31.*M43 + ...
0086              M24.*M41.*M33;
0087 
0088    I22= M33.*M44 - ...
0089              M43.*M34 - ...
0090              M31.*M44 + ...
0091              M41.*M34 + ...
0092              M31.*M43 - ...
0093              M41.*M33;
0094 
0095    I23=-M23.*M44 + ...
0096              M43.*M24 + ...
0097              M21.*M44 - ...
0098              M41.*M24 - ...
0099              M21.*M43 + ...
0100              M41.*M23;
0101 
0102    I24= M23.*M34 - ...
0103              M33.*M24 - ...
0104              M21.*M34 + ...
0105              M31.*M24 + ...
0106              M21.*M33 - ...
0107              M31.*M23;
0108 
0109    I31= M21.*M32.*M44 - ...
0110              M21.*M42.*M34 - ...
0111              M22.*M31.*M44 + ...
0112              M22.*M41.*M34 + ...
0113              M24.*M31.*M42 - ...
0114              M24.*M41.*M32;
0115 
0116    I32=-M32.*M44 + ...
0117              M42.*M34 + ...
0118              M31.*M44 - ...
0119              M41.*M34 - ...
0120              M31.*M42 + ...
0121              M41.*M32;
0122 
0123    I33= M22.*M44 - ...
0124              M42.*M24 - ...
0125              M21.*M44 + ...
0126              M41.*M24 + ...
0127              M21.*M42 - ...
0128              M41.*M22;
0129 
0130    I34=-M22.*M34 + ...
0131              M32.*M24 + ...
0132              M21.*M34 - ...
0133              M31.*M24 - ...
0134              M21.*M32 + ...
0135              M31.*M22;
0136 
0137    I41=-M21.*M32.*M43 + ...
0138              M21.*M42.*M33 + ...
0139              M22.*M31.*M43 - ...
0140              M22.*M41.*M33 - ...
0141              M23.*M31.*M42 + ...
0142              M23.*M41.*M32;
0143 
0144    I42= M32.*M43 - ...
0145              M42.*M33 - ...
0146              M31.*M43 + ...
0147              M41.*M33 + ...
0148              M31.*M42 - ...
0149              M41.*M32;
0150 
0151    I43=-M22.*M43 + ...
0152              M42.*M23 + ...
0153              M21.*M43 - ...
0154              M41.*M23 - ...
0155              M21.*M42 + ...
0156              M41.*M22;
0157 
0158    I44= M22.*M33 - ...
0159              M32.*M23 - ...
0160              M21.*M33 + ...
0161              M31.*M23 + ...
0162              M21.*M32 - ...
0163              M31.*M22;
0164 
0165    D =      I11 + ...
0166        M21.*I12 + ...
0167        M31.*I13 + ...
0168        M41.*I14;
0169 
0170    if any(abs(D) < eps) 
0171       warning('Determinant close to zero');
0172    end
0173 
0174    I = bsxfun(@times,...
0175        [I11, I12, I13, I14;
0176         I22, I22, I23, I24;
0177         I32, I32, I33, I34;
0178         I42, I42, I43, I44], (1./D));
0179 
0180 
0181 
0182 
0183 
0184 
0185 function do_unit_test
0186    M=reshape((1:16).^(0.5),4,4); M(1,:) = 1;
0187    M(2:end,:,:) = M(2:end,:,:) - ...
0188              mean(M(2:end,:,:),2);
0189    [I,D] = vectorize_4x4inv(M);
0190    iM = inv(M);
0191    unit_test_cmp('vectorize_4x4inv', ...
0192        I(:,2:end) ,iM(:,2:end),1e-09);
0193    unit_test_cmp('vectorize_4x4inv',D,det(M),1e-12);
0194    M=reshape((1:9).^(0.5),3,3); M(1,:) = 1;
0195    [I,D] = vectorize_3x3inv(M);
0196    unit_test_cmp('vectorize_3x3inv',I,inv(M),1e-12);
0197    unit_test_cmp('vectorize_3x3inv',D,det(M),1e-12);
0198
vectorize_inverse

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
