-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsmoother.m
98 lines (59 loc) · 1.36 KB
/
smoother.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
function [u,v,w]=smoother(u,v,w,U,V,W,L)
global noisy_input;
[x, y, z] = size(U);
flag = 0;
it = 0;
while flag ~= 2
if noisy_input ~= 1 && flag ~=1 && noisy_input == 1
[u,v,w] = inserter(u,v,w,U,V,W);
end
[omegax,omegay,omegaz]= omega3d_2th_order(u,v,w,x,y,z);
g1=cnv3d(x,y,z,omegax,omegay,omegaz);
Z = x*y*z;
L_E2 = speye(3*Z,3*Z)-.01*L ;
[ps1] = cgs(L_E2,g1',10e-4,100);
c = 0;
for m=1:z
for j=1:y
for i=1:x
c=c+1;
omegax(i,j,m)=ps1(c);
omegay(i,j,m)=ps1(c+x*y*z);
omegaz(i,j,m)=ps1(c+2*x*y*z);
end
end
end
[Oomegax,Oomegay,Oomegaz]= omega3d_2th_order(omegax,omegay,omegaz,x,y,z);
omx = Oomegax+u;
omy = Oomegay+v;
omz = Oomegaz+w;
[g1]=cnv3d(x,y,z,omx,omy,omz);
Z = x*y*z;
L_E2 = speye(3*Z,3*Z)-L ;
[ps1] = cgs(L_E2,g1',10e-4,100);
c = 0;
for m=1:z
for j=1:y
for i=1:x
c=c+1;
u(i,j,m)=ps1(c);
v(i,j,m)=ps1(c+x*y*z);
w(i,j,m)=ps1(c+2*x*y*z);
end
end
end
[u] = Energy(u,U);
[v] = Energy(v,V);
[w] = Energy(w,W);
[u,v,w] = inserter(u,v,w,U,V,W);
% [q,RMS1,divergence] = quality_checker(u,v,w,U,V,W)
if it == 0
it = it+1;
Mean_angle_deviation = 360;
else
[flag,Mean_angle_deviation,it]=termination_criterion(u,v,w,u_pr,v_pr,w_pr,it,Mean_angle_deviation,flag)
end
u_pr = u;
v_pr = v;
w_pr = w;
end