Skip to content

Commit 665d9d9

Browse files
committed
rough version of test data
1 parent fbd6c5c commit 665d9d9

10 files changed

Lines changed: 2069 additions & 0 deletions

_test_data/MFD_tests/Wrapping2.m

Lines changed: 168 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,168 @@
1+
% Trova le tangenti ai punti dati spostandoli in un riferimento col centro
2+
% della crf al centro degli assi. Quando non c'e' wrapping calcola la
3+
% distanza cmq tra i punti e mette come tangenti i punti iniziali. Infine
4+
% se i punti sono sulla crf non usa ottimizzazione ma identifica tangente
5+
% col punto.
6+
% lato = 1;% 1 dx o sotto---- -1 sx o sopra
7+
% problematico GA perche' sulla crf
8+
function [L,T1_glob,T2_glob,contact] = Wrapping2(centro_glob,R,P1,P2,lato,graphics)
9+
10+
% lavoro sempre con circonferenza al centro degli assi: l'algoritmo è più
11+
% preciso.
12+
v1 = P1-centro_glob;
13+
v2 = P2-centro_glob;
14+
15+
dist_P = v2-v1;
16+
if cross([1 0 0 ],[v1/norm(v1) 0])>0 == [0 0 1]
17+
% display('P1 sopra la circonferenza')
18+
else
19+
% display('P1 sotto la circonferenza')
20+
lato = lato*-1;
21+
end
22+
if cross([1 0 0 ],[v2/norm(v2) 0])<0 == [0 0 1]
23+
% display('P2 sotto la circonferenza')
24+
else
25+
% display('P2 sopra la circonferenza')
26+
end
27+
28+
if cross([1 0 0 ],[v1/norm(v1) 0])>0 == [0 0 1] & cross([1 0 0 ],[v2/norm(v2) 0])>0 == [0 0 1]
29+
% display('P1 e P2 sopra la circonferenza')
30+
lato = lato*-1;
31+
end
32+
%direzione: se P1 è sopra P2:
33+
%-----> lato = 1 wrapping a dx
34+
%-----> lato = -1 wrapping a sx
35+
% se P1 e sotto P2 è il contrario
36+
37+
%stima della distanza tra centro e vettore che unisce i punti (limit
38+
%condition)
39+
h = cross([-lato*v1,0],[dist_P,0]/norm([dist_P,0]));
40+
41+
if h(3)>R
42+
% display(' '),display('NO WRAPPING!!!'),display(' ')
43+
L = norm([dist_P,0]);
44+
contact = 0;
45+
T1_glob = P1;
46+
T2_glob = P2;
47+
if strcmp(graphics,'on') == 1
48+
%grafica normalizzata
49+
% plot(v1(1), v1(2),'o'),plot(v2(1), v2(2),'o')
50+
% plot([v1(1),v2(1)], [v1(2),v2(2)],'g')
51+
%grafica globale
52+
t = 0:(2*pi)/300:2*pi;
53+
CRF = [R*cos(t'), R*sin(t')];
54+
plot(centro_glob(:,1)+CRF(:,1), centro_glob(:,2)+CRF(:,2),'k'),hold on
55+
plot([P1(1),P2(1)], [P1(2),P2(2)],'k')
56+
plot(P1(1), P1(2),'o','MarkerEdgeColor','k','MarkerFaceColor','w');
57+
plot(P2(1), P2(2),'o','MarkerEdgeColor','k','MarkerFaceColor','w')
58+
59+
end
60+
return
61+
else
62+
% display(' '),display('WRAP!!!'),display(' ')
63+
contact = 1;
64+
end
65+
66+
%% algorithm
67+
% pos punto tangenza + vettore tangente = pos_punto est
68+
%% tangente nel verso delle theta positive se lato = 1
69+
% l'algoritmo chiude OT1P1 in pratica
70+
dist1 = sqrt(norm(v1)^2-R^2);
71+
72+
fun_min1= @(x) R*[cos(x) sin(x)]+lato*dist1*[-sin(x) cos(x)]-v1;
73+
options1=optimset('TolFun',10e-13,'TolX',10e-11);
74+
75+
x0 = pi/4;
76+
[sol1,fval,exitflag] = fsolve(fun_min1,x0,options1);
77+
% retry with new initial point
78+
if exitflag == -2
79+
x0 = sol1;
80+
[sol1,fval,exitflag] = fsolve(fun_min1,x0,options1);
81+
end
82+
% if no solution then it doesn't work
83+
if exitflag == -2
84+
display('La tangente superiore non converge!!!')
85+
waitforbuttonpress
86+
end
87+
T1 = [R*cos(sol1), R*sin(sol1)];
88+
if(dist1<10e-6)
89+
T1 = v1;
90+
% display('PUNTO SULLA CRF');
91+
% waitforbuttonpress
92+
end
93+
%% tangente nel verso opposto al precedente
94+
% l'algoritmo chiude OT1P2 in pratica
95+
%tangente nel verso delle theta negative
96+
dist2 = sqrt(norm(v2)^2-R^2);
97+
fun_min2= @(x) R*[cos(x) sin(x)]-lato*dist2*[-sin(x) cos(x)]-v2;
98+
x0 = -pi/4;
99+
options2=optimset('TolFun',10e-11,'TolX',10e-11);
100+
[sol2,fval2,exitflag2] = fsolve(fun_min2,x0,options2);
101+
% retry with new initial point
102+
if exitflag2 == -2
103+
x0 = sol+20/180*pi;
104+
[sol2,fval2,exitflag2] = fsolve(fun_min2,x0,options2);
105+
end
106+
% if no solution then it doesn't work
107+
if exitflag2 == -2
108+
display('La tangente inferiore non converge!!!')
109+
waitforbuttonpress
110+
end
111+
T2 = [R*cos(sol2), R*sin(sol2)];
112+
if(dist2<10e-6)
113+
T2 = v2;
114+
% display('PUNTO SULLA CRF');
115+
% waitforbuttonpress
116+
end
117+
%% length mi fido
118+
% arco (archi minori di 180)
119+
L_arch = R*asin(norm(cross([T1/R,0], [T2/R,0])));
120+
121+
if L_arch<0
122+
display('ARCO DI WRAPPING MINORE DI ZERO: FORSE WRAPPING NON CORRETTO?')
123+
waitforbuttonpress
124+
end
125+
%dist lineare
126+
Lin_vect1 = (v1-T1); L1 = hypot(Lin_vect1(:,1),Lin_vect1(:,2));
127+
Lin_vect2 = (v2-T2); L2 = hypot(Lin_vect2(:,1),Lin_vect2(:,2));
128+
%controllato GA e' affidabile nel wrapping
129+
% L1
130+
% L2
131+
% L_arch
132+
% waitforbuttonpress
133+
134+
L = L1+L2+L_arch;
135+
136+
%% points in global coordinates
137+
T1_glob = T1+centro_glob;
138+
T2_glob = T2+centro_glob;
139+
% local graphics of the wrapping in normalized coordinates
140+
%% grafica punti base
141+
142+
if strcmp(graphics,'on') == 1
143+
%% grafica normalizzata (no casi di no wrapping)
144+
% t = 0:(2*pi)/300:2*pi;
145+
% CRF = [R*cos(t'), R*sin(t')];
146+
% line(CRF(:,1), CRF(:,2)),hold on
147+
% % punti esterni
148+
% plot(v1(1), v1(2),'o'),plot(v2(1), v2(2),'o')
149+
% axis equal
150+
% % tangenti
151+
% plot([v1(1),T1(1)], [v1(2),T1(2)],'r')
152+
% plot([v2(1),T2(1)], [v2(2),T2(2)])
153+
%% grafica globale (no casi di no wrapping)
154+
155+
t = 0:(2*pi)/300:2*pi;
156+
CRF = [R*cos(t'), R*sin(t')];
157+
plot(centro_glob(:,1)+CRF(:,1), centro_glob(:,2)+CRF(:,2),'k'),hold on
158+
% punti esterni
159+
160+
161+
% tangenti
162+
plot([P1(1),T1_glob(1)], [P1(2),T1_glob(2)],'k');
163+
plot([P2(1),T2_glob(1)], [P2(2),T2_glob(2)],'k');
164+
plot(P1(1), P1(2),'o','MarkerEdgeColor','k','MarkerFaceColor','w');
165+
plot(P2(1), P2(2),'o','MarkerEdgeColor','k','MarkerFaceColor','w');
166+
axis equal
167+
end
168+
clc
Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
clear;clc;
2+
RotAngle = -90:2:90;
3+
4+
GLOBAL_ANAT = importdata('GLOBAL_ANATOM_MuscleForceDirection_vectors.sto');
5+
LOCAL_ANAT = importdata('LOCAL_ANATOM_MuscleForceDirection_vectors.sto');
6+
GLOBAL_EFFECT = importdata('GLOBAL_EFFECT_MuscleForceDirection_vectors.sto');
7+
LOCAL_EFFECT = importdata('LOCAL_EFFECT_MuscleForceDirection_vectors.sto');
8+
9+
%% Muscle 1: MuscleStraight
10+
%P1(0,0,0) in ground to P2(0.3, 0 ,0) in body1
11+
MUS1.v1_Global_Anat_Matlab = [cosd(RotAngle)', sind(RotAngle)',zeros(length(RotAngle),1)];
12+
% in the sto file, v1 = data(:,2:4) and v2 = -v1 = data(:,5:7)
13+
MUS1.v1_Global_Anat_OpenSim =GLOBAL_ANAT.data(:,2:4);
14+
MUS1.DiffGlobalAnat = MUS1.v1_Global_Anat_Matlab-MUS1.v1_Global_Anat_OpenSim;
15+
16+
%% Muscle 2: MuscleViaPoint
17+
% Effective: the effective attachments of muscle 2 are the same as muscle 1
18+
v1_Global_Effect_Matlab = MUS1.v1_Global_Anat_Matlab;
19+
% in the sto file, v1 = data(:,8:10) and v2 = -v1 = data(:,11:13)
20+
v1_Global_Effect_OpenSim = GLOBAL_EFFECT.data(:,8:10);
21+
DiffGlobalEffect = v1_Global_Effect_Matlab-v1_Global_Effect_OpenSim;
22+
23+
%% Muscle 3: MuscleWrap
24+
%P1(0.30000, 0.0 ,0.0) in ground to P2(-0.30000, 0.0 ,0.0)
25+
P1 = [0.30000*ones(length(RotAngle),1), 0.0*ones(length(RotAngle),1) ,0.0*ones(length(RotAngle),1)];
26+
P2 = -0.3*[cosd(RotAngle)', sind(RotAngle)',zeros(length(RotAngle),1)];
27+
28+
% [L,T1_glob,T2_glob,contact] = Wrapping2([0 0 0]',0.05,P1,P2,1,'on')
29+
30+
count = 0;
31+
for k =1:length(P1)
32+
[L(k,:),T1_glob(k,:),T2_glob(k,:),contact(k,:)] = Wrapping2([0 0 ],0.05,P1(k,1:2),P2(k,1:2),1,'on');
33+
v1(k,:) = (T1_glob(k,:)-P1(k,1:2))/norm(T1_glob(k,:)-P1(k,1:2));
34+
end
35+
36+
v1_GlobalAnat = [v1, zeros(length(P1),1)];
37+
v1_GlobalAnat_OpenSim = GLOBAL_ANAT.data(:,14:16);
38+
v1_GlobalAnat-v1_GlobalAnat_OpenSim
39+
Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
Results
2+
version=1
3+
nRows=91
4+
nColumns=2
5+
inDegrees=yes
6+
endheader
7+
time AngleRotation
8+
0 -90
9+
0.01 -88
10+
0.02 -86
11+
0.03 -84
12+
0.04 -82
13+
0.05 -80
14+
0.06 -78
15+
0.07 -76
16+
0.08 -74
17+
0.09 -72
18+
0.1 -70
19+
0.11 -68
20+
0.12 -66
21+
0.13 -64
22+
0.14 -62
23+
0.15 -60
24+
0.16 -58
25+
0.17 -56
26+
0.18 -54
27+
0.19 -52
28+
0.2 -50
29+
0.21 -48
30+
0.22 -46
31+
0.23 -44
32+
0.24 -42
33+
0.25 -40
34+
0.26 -38
35+
0.27 -36
36+
0.28 -34
37+
0.29 -32
38+
0.3 -30
39+
0.31 -28
40+
0.32 -26
41+
0.33 -24
42+
0.34 -22
43+
0.35 -20
44+
0.36 -18
45+
0.37 -16
46+
0.38 -14
47+
0.39 -12
48+
0.4 -10
49+
0.41 -8
50+
0.42 -6
51+
0.43 -4
52+
0.44 -2
53+
0.45 0
54+
0.46 2
55+
0.47 4
56+
0.48 6
57+
0.49 8
58+
0.5 10
59+
0.51 12
60+
0.52 14
61+
0.53 16
62+
0.54 18
63+
0.55 20
64+
0.56 22
65+
0.57 24
66+
0.58 26
67+
0.59 28
68+
0.6 30
69+
0.61 32
70+
0.62 34
71+
0.63 36
72+
0.64 38
73+
0.65 40
74+
0.66 42
75+
0.67 44
76+
0.68 46
77+
0.69 48
78+
0.7 50
79+
0.71 52
80+
0.72 54
81+
0.73 56
82+
0.74 58
83+
0.75 60
84+
0.76 62
85+
0.77 64
86+
0.78 66
87+
0.79 68
88+
0.8 70
89+
0.81 72
90+
0.82 74
91+
0.83 76
92+
0.84 78
93+
0.85 80
94+
0.86 82
95+
0.87 84
96+
0.88 86
97+
0.89 88
98+
0.9 90

0 commit comments

Comments
 (0)