0001 function [dY,dZ] = make_orthogonal_vec(dX)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 [N, d] = size(dX);
0013 if d ~= 3, error('data dimension error'); end
0014
0015 dY = zeros(N, d);
0016 dZ = zeros(N, d);
0017
0018
0019 xx = sum(dX.^2, 2);
0020 ix = find( xx > 0);
0021
0022 dx = dX(ix,:);
0023 xx = sqrt(xx(ix));
0024 N = size(dx,1);
0025
0026
0027 e1 = [ones(N,1), zeros(N,2)];
0028 e2 = [zeros(N,1), ones(N,1), zeros(N,1)];
0029
0030
0031 dX = vb_repmultiply(dx, 1./xx);
0032
0033
0034 dy1 = vb_cross2(dX,e1);
0035 dy2 = vb_cross2(dX,e2);
0036
0037 dyy1 = sqrt(sum(dy1.^2, 2));
0038 dyy2 = sqrt(sum(dy2.^2, 2));
0039
0040
0041 [tmp, jx] = max([dyy1, dyy2],[],2);
0042
0043 ix1 = find( jx == 1);
0044 ix2 = find( jx == 2);
0045
0046
0047 dy = zeros(N,3);
0048 dy(ix1,:) = dy1(ix1,:);
0049 dy(ix2,:) = dy2(ix2,:);
0050
0051
0052 dyy = sqrt(sum(dy.^2, 2));
0053 dy = vb_repmultiply(dy, 1./dyy);
0054
0055
0056
0057 dz = vb_cross2(dx,dy);
0058
0059
0060 dzz = sqrt(sum(dz.^2, 2));
0061 dz = vb_repmultiply(dz, 1./dzz);
0062
0063 dY(ix,:) = dy;
0064 dZ(ix,:) = dz;
0065
0066 return
0067
0068
0069 dyz = vb_cross2(dy,dz);
0070 XYZ = sum(dx.*dyz, 2);
0071
0072 ix = find(XYZ < 0);
0073
0074 dz(ix,:) = - dz(ix,:);
0075
0076 if ~isempty(ix)
0077 fprintf('# of reversed direction=%d\n',length(ix))
0078 end
0079
0080 DEBUG = 0;
0081
0082
0083 if DEBUG == 1
0084 dxerr = sum(abs(sqrt(sum(dX.^2, 2)) - 1))
0085 dyerr = sum(abs(sqrt(sum(dY.^2, 2)) - 1))
0086 dyerr = sum(abs(sqrt(sum(dZ.^2, 2)) - 1))
0087
0088 xyerr = sum(abs(sum(dX.*dY,2)))
0089 yzerr = sum(abs(sum(dY.*dZ,2)))
0090 zxerr = sum(abs(sum(dZ.*dX,2)))
0091 end