2017-07-19 24 views
1

私は2つの平面の間の角度を計算する必要があります。 1つの平面は手平面であり、2つ目は前腕平面である。私はその飛行機の法線を計算し、MATLABで式atan2(norm(cross(var.n1,var.n2)),dot(var.n1,var.n2));を使用しました。正と負のピークを特徴とする手首の屈曲/伸展角度を見たいが、この式では正のピークしか得られない。Matlab atan2 2つの3Dベクトルの間の角度

%% Script to compute the angles of wrist flexion/extension and adduction/abduction based on Vicon data 
% REF: Cheryl et al., March 2008 
% Order of the markers: 1.WRR 2.WRU 3.FAU 4.FAR 5.CMC2 6.CMC5 7.MCP5 8.MCP2 

clc 
close all 
clear all 

%% Initialization 
% dir_kinematic = input('Path of the folder containing the kinematic files: ','s'); 
dir_kinematic = 'C:\Users\Utente\Desktop\TESI\Vicon\test.mat'; 
cd(dir_kinematic); 

fileList = getAllFiles(dir_kinematic,0); % Get names of all kinematic files 
f=1; 
%% Conversion to angles 
for f = 1:length(fileList) 
    if ~isempty(strfind(fileList{f},'mat')) % Take only mat files 


     % 0. Loading 
     load(fileList{f}); 

     % 1. Filtering 
     frameRate = kinematic.framerate; 
     n = 9; 
     Wn = 2/(frameRate/2); 
     ftype = 'low'; 
     [b,a] = butter(n,Wn,ftype); 
     kinematic.x = filtfilt(b,a,kinematic.x); 
     kinematic.y = filtfilt(b,a,kinematic.y); 
     kinematic.z = filtfilt(b,a,kinematic.z); 


     % 2. Create vectors 
     var.n=length(kinematic.x); 

     % Forearm plane 
     var.FAU_WRU=[kinematic.x(:,2)-kinematic.x(:,3),kinematic.y(:,2)-kinematic.y(:,3),kinematic.z(:,2)-kinematic.z(:,3)]; 
     var.WRR_WRU=[kinematic.x(:,2)-kinematic.x(:,1),kinematic.y(:,2)-kinematic.y(:,1),kinematic.z(:,2)-kinematic.z(:,1)]; 

     % Hand plane 
     var.CMC5_MCP5=[kinematic.x(:,7)-kinematic.x(:,6),kinematic.y(:,7)-kinematic.y(:,6),kinematic.z(:,7)-kinematic.z(:,6)]; 
     var.MCP2_MCP5=[kinematic.x(:,7)-kinematic.x(:,8),kinematic.y(:,7)-kinematic.y(:,8),kinematic.z(:,7)-kinematic.z(:,8)]; 

     % Transpose 
     var.FAU_WRU = var.FAU_WRU'; 
     var.WRR_WRU = var.WRR_WRU'; 
     var.CMC5_MCP5 = var.CMC5_MCP5'; 
     var.MCP2_MCP5 = var.MCP2_MCP5'; 

     % 3. Calculate angle of wrist flexion/extension 
     % Cross vector function for all time => create normal vector plane 
     var.forearm_n=[]; 
     var.hand_n=[]; 
     var.theta_rad=[]; 

     for i = 1:var.n % Loop through experiment 

      % vector x and y of the forearm plane 
      var.v1=var.FAU_WRU(:,i); % take x,y,z of the vector for every time 
      var.v2=var.WRR_WRU(:,i); 

      % vector x and y of the hand plane 
      var.v3=var.CMC5_MCP5(:,i); 
      var.v4=var.MCP2_MCP5(:,i); 

      var.forearm_n= [var.forearm_n, cross(var.v1,var.v2)]; 
      var.hand_n=[var.hand_n, cross(var.v3,var.v4)]; 

     end 

     % Calculate angle 
     for i = 1:var.n 

      var.n1=(var.forearm_n(:,i)); 
      var.n2=var.hand_n(:,i); 

      var.scalar_product(i) = dot(var.n1,var.n2); 

      %Equation (2) of the paper 
      var.theta_rad=[var.theta_rad, atan2(norm(cross(var.n1,var.n2)),dot(var.n1,var.n2))]; % result in radian 

      angle.flex_deflex_wrist{f}=(var.theta_rad*180)/pi; 

     end 

     % 4. Calculate angle of wrist adduction/abduction 

     % Projection vector onto plane 

     var.MCP2_MCP5_forearmproj=[]; 
     var.WRR_WRU_forearmproj=[]; 
     var.rad_ul_angle_rad=[]; 

     for i=1:var.n 

      %take x,y,z of the vector for each time 
      var.v1=var.MCP2_MCP5(:,i); 
      var.v2=var.WRR_WRU(:,i); 

      % vector x and y of the forearm plane 
      var.vfx=var.FAU_WRU(:,i); % take x,y,z of the vector for every time 
      var.vfy=var.WRR_WRU(:,i); 

      %projection of vector MCP2_MCP5 and WRR_WRU onto forearm plane 
      var.squNorm1=(norm(var.vfx)*norm(var.vfx)); 
      var.squNorm2=(norm(var.vfy)*norm(var.vfy)); 

      var.MCP2_MCP5_forearmproj=[var.MCP2_MCP5_forearmproj,((((var.v1')*var.vfx)*var.vfx/var.squNorm1)+(((var.v1')*var.vfy)*var.vfy/var.squNorm2))]; 
      var.WRR_WRU_forearmproj=[var.WRR_WRU_forearmproj,((var.vfx*((var.v2')*var.vfx/var.squNorm1))+(var.vfy*((var.v2')*var.vfy/var.squNorm2)))]; 
     end 

     % Calculate angle 

     for i=1:var.n 

      var.n1=var.MCP2_MCP5_forearmproj(:,i)'; 
      var.n2=var.WRR_WRU_forearmproj(:,i); 
      var.product=var.n1*var.n2; 

      var.rad_ul_angle_rad=[var.rad_ul_angle_rad, atan2(norm(cross(var.n1,var.n2)),dot(var.n1,var.n2))];% en rad 
      angle.rad_ul_wrist{f}=(var.rad_ul_angle_rad*180)/pi; 
     end 

    end 
end 

私のアングルが常にポジティブである理由を知りたいですか?私は正と負のピークを見る必要があります... 助けてくれてありがとう!

答えて

1

2つの平面の間の角度は法線のそれと同じですので、私は2つのベクトルの間の角度をに制限します。

我々は、二つのベクトル(B)間のクロス積は、それらの両方に垂直な別のベクターであることを知っています。

私たちは角度Tと+ Tを区別する問題に取り組んでいます。ここでTはある角度です。 ABは同じである一方であるため

これは(、AB | | AXB) にatan2:あなたが使用される式を使用して

、これら二つの角度が原因で使用される基本的な式自体に、同じ結果を与えるだろうどちらの場合も、axbはとちょうどの2つのベクトルの法線の符号が異なります。このベクトルのノルムを計算すると、その符号に関する情報が失われるため、関数は常に正の値を返します。

あなたはそれがあなたの角度が正または負であるかどうかを判断するために、X bの符号を追跡する必要が

について何ができますか。

注:携帯電話から返信しているため、フォーマットやコードを追加できないため、まもなくその回答が更新されます。

関連する問題