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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190 | #include "GreenCoordinates3D.hpp"
#include <iostream>
#include <math.h>
#include <map>
namespace MeshKit {
GreenCoordinates3D::GreenCoordinates3D(MKCore* core, vector< vector<double> > iNodes)
{
mk_core = core;
InteriorNodes.resize(iNodes.size());
for (unsigned int i = 0; i < iNodes.size(); i++){
InteriorNodes[i].resize(3);
for (int j = 0; j < 3; j++)
InteriorNodes[i][j] = iNodes[i][j];
}
}
//Setup the cages for surface mesh
void GreenCoordinates3D::SetupCages(vector<vector<double> > cageNodes, vector< vector<int> > cageFaces, vector<vector<double> > normal)<--- The function 'SetupCages' is never used.
{
int num_CageNodes = cageNodes.size(), num_CageFaces = cageFaces.size();
CageNodes.resize(num_CageNodes);
CageFaces.resize(num_CageFaces);
for (int i = 0; i < num_CageNodes; i++){
int size_j = cageNodes[i].size();
CageNodes.resize(size_j);
for (int j = 0; j < size_j; j++)
CageNodes[i][j] = cageNodes[i][j];
}
for (int i = 0; i < num_CageFaces; i++){
int size_j = CageFaces[i].size();
CageFaces.resize(size_j);
for (int j = 0; j < size_j; j++)
CageFaces[i][j] = cageFaces[i][j];
}
Normals.resize(normal.size());
for (unsigned int i = 0; i < normal.size(); i++){
unsigned int size_j = normal[i].size();
Normals[i].resize(size_j);
for (unsigned int j = 0; j < size_j; j++)
Normals[i][j] = normal[i][j];
}
}
//calculate the green coordinates
void GreenCoordinates3D::Execute()
{
//Initialization
double error = 1.0e-5;
weight_t.resize(InteriorNodes.size());
weight_v.resize(InteriorNodes.size());
for (unsigned int i = 0; i < weight_t.size(); i++){
weight_t[i].resize(CageFaces.size());
for (unsigned int j = 0; j < CageFaces.size(); j++)
weight_t[i][j] = 0.0;
}
for (unsigned int i = 0; i < weight_v.size(); i++){
weight_v[i].resize(CageNodes.size());
for (unsigned int j = 0; j < CageNodes.size(); j++)
weight_v[i][j] = 0.0;
}
//Main Loop--Calculate the weight for cage nodes and cage face normal
//From the appendix of the paper: Green Coordinates, Yaron Lipman, David Levin, Daniel Cohen-Or
for (unsigned int i = 0; i < InteriorNodes.size(); i++){
for (unsigned int j = 0; j < CageFaces.size(); j++){
vector<vector<double> > node(3, vector<double>(3));
for (int k = 0; k < 3; k++){
node[0][k] = CageNodes[CageFaces[j][0]][k] - InteriorNodes[i][k];
node[1][k] = CageNodes[CageFaces[j][1]][k] - InteriorNodes[i][k];
node[2][k] = CageNodes[CageFaces[j][2]][k] - InteriorNodes[i][k];
}
double mag = node[0][0]*Normals[j][0] + node[1][1]*Normals[j][1] + node[2][2]*Normals[j][2];
vector<double> p(3);
for (int k = 0; k < 3; k++)
p[k] = mag*Normals[j][k];
vector<double> s(3),I(3),II(3);
vector<vector<double> > q(3, vector<double>(3)), N(3, vector<double>(3));
for (int k = 0; k < 3; k++){
vector<double> ivector(3), jvector(3);
for (int l = 0; l < 3; l++){
ivector[l] = node[k][l] - p[l];
jvector[l] = node[(k+1)%3][l]- p[l];
}
vector<double> crossproduct(3);
crossproduct[0] = ivector[1]*jvector[2] - ivector[2]*jvector[1];
crossproduct[1] = -1.0*(ivector[0]*jvector[2]-ivector[2]*jvector[0]);
crossproduct[2] = ivector[0]*jvector[1] - ivector[1]*jvector[0];
double dotproduct = crossproduct[0]*Normals[j][0]+crossproduct[1]*Normals[j][1]+crossproduct[2]*Normals[j][2];
s[k] = dotproduct>0? 1.0:-1.0;
vector<double> tmp(3);
tmp[0] = 0.0; tmp[1] = 0.0; tmp[2] = 0.0;
I[k] = GCTriInt(p, node[k], node[(k+1)%3], tmp);
II[k] = GCTriInt(tmp, node[k], node[(k+1)%3], tmp);
q[k][0] = node[(k+1)%3][1]*node[k][2]-node[(k+1)%3][2]*node[k][1];
q[k][1] = -1.0*(node[(k+1)%3][0]*node[k][2]-node[(k+1)%3][2]*node[k][0]);
q[k][2] = node[(k+1)%3][0]*node[k][1]-node[(k+1)%3][1]*node[k][0];
mag = sqrt(pow(q[k][0],2) + pow(q[k][1],2) + pow(q[k][2], 2));
N[k][0] = q[k][0]/mag;
N[k][1] = q[k][1]/mag;
N[k][2] = q[k][2]/mag;
}
double I_scalar = -1.0*fabs(s[0]*I[0]+s[1]*I[1]+s[2]*I[2]);
weight_t[i][j] = -1.0*I_scalar;
vector<double> w(3);
w[0] = I_scalar*Normals[j][0] + N[0][0]*I_scalar*I[0] + N[1][0]*I_scalar*I[1] + N[2][0]*I_scalar*I[2];
w[1] = I_scalar*Normals[j][1] + N[0][1]*I_scalar*I[0] + N[1][1]*I_scalar*I[1] + N[2][1]*I_scalar*I[2];
w[2] = I_scalar*Normals[j][2] + N[0][2]*I_scalar*I[0] + N[1][2]*I_scalar*I[1] + N[2][2]*I_scalar*I[2];
double w_mag = sqrt(pow(w[0],2)+pow(w[1],2)+pow(w[2],2));
if (w_mag > error)//this guarantees the local deformation property
for (int k = 0; k < 3; k++)
weight_v[i][CageFaces[j][k]] += (N[(k+1)%3][0]*w[0]+N[(k+1)%3][1]*w[1]+N[(k+1)%3][2]*w[2])/(N[(k+1)%3][0]*node[k][0]+N[(k+1)%3][1]*node[k][1]+N[(k+1)%3][2]*node[k][2]);
}
}
}
double GreenCoordinates3D::GCTriInt(vector<double> p, vector<double> v1, vector<double> v2, vector<double> iNodes)
{
//calculate the angle
double dotproduct_21p1 = (v2[0]-v1[0])*(p[0]-v1[0]) + (v2[1]-v1[1])*(p[1]-v1[1]) + (v2[2]-v1[2])*(p[2]-v1[2]);
double dotproduct_1p2p = (v1[0]-p[0])*(v2[0]-p[0]) + (v1[1]-p[1])*(v2[1]-p[1]) + (v1[2]-p[2])*(v2[2]-p[2]);
double mag_21 = sqrt((v2[0]-v1[0])*(v2[0]-v1[0])+(v2[1]-v1[1])*(v2[1]-v1[1])+(v2[2]-v1[2])*(v2[2]-v1[2]));
double mag_p1 = sqrt((p[0]-v1[0])*(p[0]-v1[0])+(p[1]-v1[1])*(p[1]-v1[1])+(p[2]-v1[2])*(p[2]-v1[2]));
double mag_p2 = sqrt((p[0]-v2[0])*(p[0]-v2[0])+(p[1]-v2[1])*(p[1]-v2[1])+(p[2]-v2[2])*(p[2]-v2[2]));
double alpha = acos(dotproduct_21p1/(mag_21*mag_p1));
double belta = acos(dotproduct_1p2p/(mag_p1*mag_p2));
double lambda = pow(mag_p1,2)*pow(sin(alpha), 2);
double c = pow(p[0]-iNodes[0],2) + pow(p[1]-iNodes[1],2) + pow(p[2]-iNodes[2],3);
double pi = atan(1.0)*4.0;
double S = sin(pi-alpha), C = cos(pi-alpha);
double theta_pi_alpha = -1.0*(S >0? 1.0:-1.0)*0.5*(2.0*sqrt(C)*atan(sqrt(c)*C/sqrt(lambda+S*S*c)) + sqrt(lambda)*log(2*sqrt(lambda)*pow(S,2)*(1-2*c*C/(c*(1+C)+lambda+sqrt(pow(lambda,2)+lambda*c*pow(S,2))))/pow(1-C,2)));
S = sin(pi-alpha-belta);
C = cos(pi-alpha-belta);
double theta_pi_alpha_belta = -1.0*(S >0? 1.0:-1.0)*0.5*(2.0*sqrt(C)*atan(sqrt(c)*C/sqrt(lambda+S*S*c)) + sqrt(lambda)*log(2*sqrt(lambda)*pow(S,2)*(1-2*c*C/(c*(1+C)+lambda+sqrt(pow(lambda,2)+lambda*c*pow(S,2))))/pow(1-C,2)));
return -1.0*fabs(theta_pi_alpha-theta_pi_alpha_belta-sqrt(c)*belta)/(4*pi);
}
//input new cage and calculate deformed interior points based on the existing weights for cage vertices and cage face normals
void GreenCoordinates3D::GetDeformedVertices(vector< vector<double> > nodes, vector< vector<double> > norm, vector<vector<double> > &ReturnNodes)<--- The function 'GetDeformedVertices' is never used.
{
ReturnNodes.resize(InteriorNodes.size());
//loop over the interior nodes
for (unsigned int i = 0; i < InteriorNodes.size(); i++){
ReturnNodes[i].resize(3);
for (int j = 0; j < 3; j++)
ReturnNodes[i][j] = 0.0;
//consider the effect of cage nodes
for (unsigned int j = 0; j < nodes.size(); j++)
for (int k = 0; k < 3; k++)
ReturnNodes[i][k] += weight_v[i][j]*nodes[j][k];
//consider the effect of cage face normals
for (unsigned int j = 0; j < norm.size(); j++)
for (int k = 0; k < 3; k++)
ReturnNodes[i][k] += weight_t[i][j]*norm[j][k];
}
}
GreenCoordinates3D::~GreenCoordinates3D()
{
cout << "It is over now in calculating the green coordinates for " << endl;
}
}
|