51 |
|
insert_vert(vert, miga[i]->rbfv[1]); |
52 |
|
} |
53 |
|
/* should be just 3 vertices */ |
54 |
< |
if ((vert[3] == NULL) | (vert[4] != NULL)) |
54 |
> |
if ((vert[2] == NULL) | (vert[3] != NULL)) |
55 |
|
return(0); |
56 |
|
/* identify edge 0 */ |
57 |
|
for (i = 3; i--; ) |
85 |
|
return(1); |
86 |
|
} |
87 |
|
|
88 |
< |
/* Determine if we are close enough to the given edge */ |
88 |
> |
/* Determine if we are close enough to an edge */ |
89 |
|
static int |
90 |
|
on_edge(const MIGRATION *ej, const FVECT ivec) |
91 |
|
{ |
92 |
< |
double cos_a = DOT(ej->rbfv[0]->invec, ivec); |
93 |
< |
double cos_b = DOT(ej->rbfv[1]->invec, ivec); |
94 |
< |
double cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec); |
95 |
< |
double cos_aplusb = cos_a*cos_b - |
96 |
< |
sqrt((1.-cos_a*cos_a)*(1.-cos_b*cos_b)); |
92 |
> |
double cos_a, cos_b, cos_c, cos_aplusb; |
93 |
> |
/* use triangle inequality */ |
94 |
> |
cos_a = DOT(ej->rbfv[0]->invec, ivec); |
95 |
> |
if (cos_a <= 0) |
96 |
> |
return(0); |
97 |
|
|
98 |
< |
return(cos_aplusb - cos_c < .01); |
98 |
> |
cos_b = DOT(ej->rbfv[1]->invec, ivec); |
99 |
> |
if (cos_b <= 0) |
100 |
> |
return(0); |
101 |
> |
|
102 |
> |
cos_aplusb = cos_a*cos_b - sqrt((1.-cos_a*cos_a)*(1.-cos_b*cos_b)); |
103 |
> |
if (cos_aplusb <= 0) |
104 |
> |
return(0); |
105 |
> |
|
106 |
> |
cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec); |
107 |
> |
|
108 |
> |
return(cos_c - cos_aplusb < .001); |
109 |
|
} |
110 |
|
|
111 |
|
/* Determine if we are inside the given triangle */ |
126 |
|
return(sgn2 == sgn3); |
127 |
|
} |
128 |
|
|
129 |
+ |
/* Test and set for edge */ |
130 |
+ |
static int |
131 |
+ |
check_edge(unsigned char *emap, int nedges, const MIGRATION *mig, int mark) |
132 |
+ |
{ |
133 |
+ |
int ejndx, bit2check; |
134 |
+ |
|
135 |
+ |
if (mig->rbfv[0]->ord > mig->rbfv[1]->ord) |
136 |
+ |
ejndx = mig->rbfv[1]->ord + (nedges-1)*mig->rbfv[0]->ord; |
137 |
+ |
else |
138 |
+ |
ejndx = mig->rbfv[0]->ord + (nedges-1)*mig->rbfv[1]->ord; |
139 |
+ |
|
140 |
+ |
bit2check = 1<<(ejndx&07); |
141 |
+ |
|
142 |
+ |
if (emap[ejndx>>3] & bit2check) |
143 |
+ |
return(0); |
144 |
+ |
if (mark) |
145 |
+ |
emap[ejndx>>3] |= bit2check; |
146 |
+ |
return(1); |
147 |
+ |
} |
148 |
+ |
|
149 |
|
/* Compute intersection with the given position over remaining mesh */ |
150 |
|
static int |
151 |
|
in_mesh(MIGRATION *miga[3], unsigned char *emap, int nedges, |
153 |
|
{ |
154 |
|
MIGRATION *ej1, *ej2; |
155 |
|
RBFNODE *tv; |
126 |
– |
int ejndx; |
156 |
|
/* check visitation record */ |
157 |
< |
if (mig->rbfv[0]->ord > mig->rbfv[1]->ord) |
129 |
< |
ejndx = mig->rbfv[1]->ord + (nedges-1)*mig->rbfv[0]->ord; |
130 |
< |
else |
131 |
< |
ejndx = mig->rbfv[0]->ord + (nedges-1)*mig->rbfv[1]->ord; |
132 |
< |
if (emap[ejndx>>3] & 1<<(ejndx&07)) /* tested already? */ |
157 |
> |
if (!check_edge(emap, nedges, mig, 1)) |
158 |
|
return(0); |
134 |
– |
emap[ejndx>>3] |= 1<<(ejndx&07); /* else mark & test it */ |
159 |
|
if (on_edge(mig, ivec)) { |
160 |
|
miga[0] = mig; /* close enough to edge */ |
161 |
|
return(1); |
163 |
|
/* do triangles either side */ |
164 |
|
for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL; |
165 |
|
ej1 = nextedge(mig->rbfv[0],ej1)) { |
166 |
< |
if (ej1 == mig) |
167 |
< |
continue; |
168 |
< |
tv = opp_rbf(mig->rbfv[0],ej1); |
169 |
< |
for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2)) |
170 |
< |
if (opp_rbf(tv,ej2) == mig->rbfv[1]) { |
171 |
< |
if (in_mesh(miga, emap, nedges, ivec, ej1)) |
172 |
< |
return(1); |
173 |
< |
if (in_mesh(miga, emap, nedges, ivec, ej2)) |
174 |
< |
return(1); |
175 |
< |
if (in_tri(mig->rbfv[0], mig->rbfv[1], |
176 |
< |
tv, ivec)) { |
177 |
< |
miga[0] = mig; |
178 |
< |
miga[1] = ej1; |
179 |
< |
miga[2] = ej2; |
180 |
< |
return(1); |
181 |
< |
} |
166 |
> |
if (ej1 == mig) |
167 |
> |
continue; |
168 |
> |
tv = opp_rbf(mig->rbfv[0],ej1); |
169 |
> |
for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2)) |
170 |
> |
if (opp_rbf(tv,ej2) == mig->rbfv[1]) { |
171 |
> |
int do_ej1 = check_edge(emap, nedges, ej1, 0); |
172 |
> |
int do_ej2 = check_edge(emap, nedges, ej2, 0); |
173 |
> |
if (do_ej1 && in_mesh(miga, emap, nedges, ivec, ej1)) |
174 |
> |
return(1); |
175 |
> |
if (do_ej2 && in_mesh(miga, emap, nedges, ivec, ej2)) |
176 |
> |
return(1); |
177 |
> |
/* check just once */ |
178 |
> |
if (do_ej1 & do_ej2 && in_tri(mig->rbfv[0], |
179 |
> |
mig->rbfv[1], tv, ivec)) { |
180 |
> |
miga[0] = mig; |
181 |
> |
miga[1] = ej1; |
182 |
> |
miga[2] = ej2; |
183 |
> |
return(1); |
184 |
|
} |
185 |
+ |
} |
186 |
|
} |
187 |
< |
return(0); |
187 |
> |
return(0); /* not near this edge */ |
188 |
|
} |
189 |
|
|
190 |
|
/* Find edge(s) for interpolating the given vector, applying symmetry */ |
359 |
|
for (j = 0; j < mtx_ncols(miga[0]); j++) |
360 |
|
for (k = (mtx_coef(miga[0],i,j) > FTINY) * |
361 |
|
mtx_ncols(miga[2]); k--; ) |
362 |
< |
n += (mtx_coef(miga[2],i,k) > FTINY && |
362 |
> |
n += (mtx_coef(miga[2],i,k) > FTINY || |
363 |
|
mtx_coef(miga[1],j,k) > FTINY); |
364 |
|
#ifdef DEBUG |
365 |
|
fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n", |
402 |
|
double rad2k; |
403 |
|
FVECT vout; |
404 |
|
int pos[2]; |
405 |
< |
if ((mb <= FTINY) | (mc <= FTINY)) |
405 |
> |
if ((mb <= FTINY) & (mc <= FTINY)) |
406 |
|
continue; |
407 |
|
rbf2k = &miga[2]->rbfv[1]->rbfa[k]; |
408 |
|
rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact); |