make_orthogonal_vec_old

 make orthogonal vector to dX in 3D space 
  [dY,dZ] = make_orthogonal_vec(dX)
 dY & dZ is orthogonal to dX and consist 2D-axis of orthogonal plane
  dX , dY , dZ: [N x 3]

 2014-11-10 Masa-aki Sato

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [dY,dZ] = make_orthogonal_vec(dX)
0002 % make orthogonal vector to dX in 3D space
0003 %  [dY,dZ] = make_orthogonal_vec(dX)
0004 % dY & dZ is orthogonal to dX and consist 2D-axis of orthogonal plane
0005 %  dX , dY , dZ: [N x 3]
0006 %
0007 % 2014-11-10 Masa-aki Sato
0008 %
0009 % Copyright (C) 2011, ATR All Rights Reserved.
0010 % License : New BSD License(see VBMEG_LICENSE.txt)
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 %% Zero vector check
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 % reference vector
0027 e1 = [ones(N,1), zeros(N,2)];
0028 e2 = [zeros(N,1), ones(N,1), zeros(N,1)];
0029 
0030 % normalization
0031 dX  = vb_repmultiply(dX, 1./xx);
0032 
0033 % make orthogonal vector to dx by cross product
0034 dy1 = vb_cross2(dX,e1);
0035 dy2 = vb_cross2(dX,e2);
0036 % norm of dy1 & dy2
0037 dyy1 = sqrt(sum(dy1.^2, 2));
0038 dyy2 = sqrt(sum(dy2.^2, 2));
0039 
0040 % select vector with larger norm
0041 [tmp, jx] = max([dyy1, dyy2],[],2);
0042 
0043 ix1 = find( jx == 1);
0044 ix2 = find( jx == 2);
0045 
0046 % dy : orthogonal vector to dx
0047 dy = zeros(N,3);
0048 dy(ix1,:) = dy1(ix1,:);
0049 dy(ix2,:) = dy2(ix2,:);
0050 
0051 % normalization
0052 dyy = sqrt(sum(dy.^2, 2));
0053 dY  = vb_repmultiply(dy, 1./dyy);
0054 
0055 % dY & dZ is orthogonal to dX and consist 2D-axis of orthogonal plane
0056 % dx,dy,dz make right-hand coordinate
0057 dz = vb_cross2(dX,dy);
0058 
0059 % normalization
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 %--------------- END -----------------
0068 % check right-hand coordinate
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 % check
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