79 #define GENERAL_ERROR 1 87 #define EPSILON 1.0e-10 89 #define EPSILON_T 1.0e-5 90 #define EPSILON_U 1.0e-5 91 #define EPSILON_V 1.0e-5 99 #define B_0(T) (1.0-2.0*(T)+(T)*(T)) 100 #define B_1(T) (2.0*(1.0-(T))*(T)) 101 #define B_2(T) ((T)*(T)) 103 #define dB_0(T) (-2.0+2.0*(T)) 104 #define dB_1(T) (2.0-4.0*(T)) 105 #define dB_2(T) (2.0*(T)) 107 #define swap(VAL1, VAL2, TMP) {(TMP) = (VAL1); (VAL1) = (VAL2); (VAL2) = (TMP);} 109 #define copy_vec(DEST, SRC) {(DEST).x = (SRC).x; (DEST).y = (SRC).y; (DEST).z = (SRC).z;} 110 #define size_vec(VEC) ((VEC).x * (VEC).x + (VEC).y * (VEC).y + (VEC).z * (VEC).z) 112 #define dist_point(PNT1, PNT2) (((PNT1).x - (PNT2).x) * ((PNT1).x - (PNT2).x) + \ 113 ((PNT1).y - (PNT2).y) * ((PNT1).y - (PNT2).y) + \ 114 ((PNT1).z - (PNT2).z) * ((PNT1).z - (PNT2).z)) 115 #define aver_point(RES, PNT1, PNT2) {(RES).x = ((PNT1).x + (PNT2).x) / 2.0; \ 116 (RES).y = ((PNT1).y + (PNT2).y) / 2.0; \ 117 (RES).z = ((PNT1).z + (PNT2).z) / 2.0;} 119 #define mul_vec(VEC, VALUE) {(VEC).x *= (VALUE); (VEC).y *= (VALUE); (VEC).z *= (VALUE);} 121 #define div_vec(VEC, VALUE) {(VEC).x /= (VALUE); (VEC).y /= (VALUE); (VEC).z /= (VALUE);} 123 #define add_vec(RES, VEC1, VEC2) {(RES).x = (VEC1).x + (VEC2).x; \ 124 (RES).y = (VEC1).y + (VEC2).y; \ 125 (RES).z = (VEC1).z + (VEC2).z;} 127 #define sub_vec(RES, VEC1, VEC2) {(RES).x = (VEC1).x - (VEC2).x; \ 128 (RES).y = (VEC1).y - (VEC2).y; \ 129 (RES).z = (VEC1).z - (VEC2).z;} 132 #define dot_product(VEC1, VEC2) ((VEC1).x * (VEC2).x + (VEC1).y * (VEC2).y + (VEC1).z * (VEC2).z) 134 #define cross_product(RES, VEC1, VEC2) {(RES).x = (VEC1).y * (VEC2).z - (VEC1).z * (VEC2).y; \ 135 (RES).y = (VEC1).z * (VEC2).x - (VEC1).x * (VEC2).z; \ 136 (RES).z = (VEC1).x * (VEC2).y - (VEC1).y * (VEC2).x;} 139 #define array(MAT,ROWS,COLS,ROW,COL) MAT[(ROW)*(COLS)+(COL)] 143 #define expand_ellipsoid(ELLIPSOID) {(ELLIPSOID).max += epsilon; \ 144 (ELLIPSOID).min += epsilon; \ 145 (ELLIPSOID).mid += epsilon;} 146 #define shrink_ellipsoid(ELLIPSOID) {(ELLIPSOID).max -= epsilon; \ 147 (ELLIPSOID).min -= epsilon; \ 148 (ELLIPSOID).mid -= epsilon;} 150 #define expand_ellipse(ELLIPSE) {(ELLIPSE).max += epsilon; \ 151 (ELLIPSE).min += epsilon;} 152 #define shrink_ellipse(ELLIPSE) {(ELLIPSE).max -= epsilon; \ 153 (ELLIPSE).min -= epsilon;} 373 static void error_message (
int exit_code,
const char *format, ...);
878 if (a == 0.0)
return true;
879 if (-zero<a && a<zero)
return true;
934 for (
int i=0; i<3; i++) {
956 point.
x = L.
max; point.
y = 0.0; point.
z = 0.0;
960 point.
x = L.
max; point.
y = L.
min; point.
z = 0.0;
966 point.
x = L.
max; point.
y = 0.0; point.
z = L.
mid;
969 point.
x = 0.0; point.
y = L.
min; point.
z = 0.0;
972 point.
x = 0.0; point.
y = L.
min; point.
z = L.
mid;
975 point.
x = 0.0; point.
y = 0.0; point.
z = L.
mid;
10167 double max, max1, max2;
10168 double dist, dist1, dist2;
10169 double u1, uu1, uuu1, u2, uu2, uuu2, delta_u =
EPSILON_U;
10170 double v1, vv1, vvv1, v2, vv2, vvv2, delta_v =
EPSILON_V;
10171 int oo1, ooo1, oo2, ooo2;
10174 if(ellipsoid1 ->
id == ellipsoid2 ->
id)
return false;
10176 max1 = ellipsoid1 -> max +
epsilon;
10177 max2 = ellipsoid2 -> max +
epsilon;
10181 dist =
dist_point(ellipsoid1 -> center, ellipsoid2 -> center);
10182 if (dist > max * max)
return false;
10192 dist1 =
dist_point(ellipsoid1 -> center, pnt2);
10193 if(dist1 > max1 * max1)
return false;
10200 dist2 =
dist_point(ellipsoid2 -> center, pnt1);
10201 if(dist2 > max2 * max2)
return false;
10212 if(fabs(u1 - uu1) < delta_u && fabs(v1 - vv1) < delta_v)
break;
10221 if(fabs(u2 - uu2) < delta_u && fabs(v2 - vv2) < delta_v)
break;
10229 if(oo1 != ooo1 || oo2 != ooo2)
return true;
10231 u1 = (uu1 + uuu1) / 2.0;
10232 v1 = (vv1 + vvv1) / 2.0;
10233 u2 = (uu2 + uuu2) / 2.0;
10234 v2 = (vv2 + vvv2) / 2.0;
10507 double x_rate, y_rate, z_rate, value;
10511 copy_vec(trans.
z, ellipsoid -> middle);
10515 x_rate = pnt.
x / (ellipsoid -> max +
epsilon);
10516 y_rate = pnt.
y / (ellipsoid -> min +
epsilon);
10517 z_rate = pnt.
z / (ellipsoid -> mid +
epsilon);
10519 value = x_rate * x_rate + y_rate * y_rate + z_rate * z_rate - 1.0;
10520 if (value <= 0.0)
return true;
10948 point_rec grad_u, grad_v, vec, point_spec, pnt_spec, normal, p;
10950 double xx, yy, zz, dist, norm, rate, rate1 = 0.25, rate2 = 1.5, delta_u =
EPSILON_U, delta_v =
EPSILON_V;
10951 double u = 0.5, v = 0.5, du, dv, duu, dvv, ddu, ddv, duv, dvu, d_u = 0.0, d_v = 0.0, max_du = 0.25, max_dv = 0.25;
10952 double norm_u, norm_v, det_x, det_y, det_z, abs_det_x, abs_det_y, abs_det_z, delta_x, delta_y, delta_z, dd;
10953 int sign = 0, sign_x = 0, sign_y = 0, sign_z = 0, iter = 0, octant = 0;
10957 copy_vec(trans.
z, ellipsoid -> middle);
10966 if(xx == 0.0 && yy == 0.0 && zz == 0.0)
10969 if(ellipsoid -> max != ellipsoid -> min && ellipsoid -> max != ellipsoid -> mid){
10970 if(yy == 0.0 && zz == 0.0){
10971 rate = ellipsoid -> min / ellipsoid -> max;
10972 if(fabs(xx) < ellipsoid -> max * (1.0 - rate * rate))
10975 if(yy == 0.0 && xx == 0.0){
10976 rate = ellipsoid -> min / ellipsoid -> mid;
10977 if(fabs(zz) < ellipsoid -> mid * (1.0 - rate * rate))
10986 if(xx != 0.0)sign_x = (xx > 0.0) ? 1 : -1;
10987 if(yy != 0.0)sign_y = (yy > 0.0) ? 1 : -1;
10988 if(zz != 0.0)sign_z = (zz > 0.0) ? 1 : -1;
11019 if(sign_y == 0 && sign_z == 0)sign = sign_x;
11020 if(sign_z == 0 && sign_x == 0)sign = sign_y;
11021 if(sign_x == 0 && sign_y == 0)sign = sign_z;
11041 if(norm_v <
EPSILON && u < delta_u / 2.0){
11042 if(d_u < delta_u)
break;
11044 d_u -= (delta_u / 2.0 - u);
11053 abs_det_x = fabs(det_x = normal.
x);
11054 abs_det_y = fabs(det_y = normal.
y);
11055 abs_det_z = fabs(det_z = normal.
z);
11057 norm = det_x * det_x + det_y * det_y + det_z * det_z;
11059 sub_vec(vec, point_spec, pnt_spec);
11062 delta_x = point_spec.
x - dist * normal.
x - pnt_spec.
x;
11063 delta_y = point_spec.
y - dist * normal.
y - pnt_spec.
y;
11064 delta_z = point_spec.
z - dist * normal.
z - pnt_spec.
z;
11066 if(abs_det_x > abs_det_y && abs_det_x > abs_det_z){
11067 du = (+ grad_v.
z * delta_y - grad_v.
y * delta_z) / det_x;
11068 dv = (- grad_u.
z * delta_y + grad_u.
y * delta_z) / det_x;
11071 if(abs_det_y > abs_det_z){
11072 du = (+ grad_v.
x * delta_z - grad_v.
z * delta_x) / det_y;
11073 dv = (- grad_u.
x * delta_z + grad_u.
z * delta_x) / det_y;
11076 du = (+ grad_v.
y * delta_x - grad_v.
x * delta_y) / det_z;
11077 dv = (- grad_u.
y * delta_x + grad_u.
x * delta_y) / det_z;
11081 if(fabs(du) < delta_u && fabs(dv) < delta_v)
break;
11086 if(ddu > max_du || ddv > max_dv){
11089 dd = (ddu > ddv) ? ddu : ddv;
11095 dd = du * d_u + dv * d_v;
11097 ddu = du + rate1 * d_u;
11098 ddv = dv + rate1 * d_v;
11099 if(d_u * ddu + d_v * ddv < 0.0){
11100 dd = -dd / rate1 / (d_u * d_u + d_v * d_v);
11106 ddu = du - rate2 * d_u;
11107 ddv = dv - rate2 * d_v;
11108 if(d_u * ddu + d_v * ddv > 0.0){
11109 dd = dd / rate2 / (d_u * d_u + d_v * d_v);
11116 duu = dvu = d_u = du;
11117 dvv = duv = d_v = dv;
11121 duu = (0.0 - u) / 2.0;
11122 if(du < -delta_u)duv = duu / du * dv;
11127 duu = (1.0 - u) / 2.0;
11128 if(du > delta_v)duv = duu / du * dv;
11134 dvv = (0.0 - v) / 2.0;
11135 if(dv < -delta_v)dvu = dvv / dv * du;
11140 dvv = (1.0 - v) / 2.0;
11141 if(dv > delta_v)dvu = dvv / dv * du;
11150 u += (d_u = (fabs(duu) < fabs(dvu)) ? duu : dvu);
11151 v += (d_v = (fabs(dvv) < fabs(duv)) ? dvv : duv);
11153 if(u < 0.0)u = 0.0;
11154 if(v < 0.0)v = 0.0;
11155 if(u > 1.0)u = 1.0;
11156 if(v > 1.0)v = 1.0;
11202 copy_vec(trans.
z, ellipsoid -> middle);
11205 if(octant == 2 || octant == 3 || octant == 6 || octant == 7)p.
x *= -1.0;
11206 if(octant == 4 || octant == 3 || octant == 8 || octant == 7)p.
y *= -1.0;
11207 if(octant == 5 || octant == 6 || octant == 7 || octant == 8)p.
z *= -1.0;
11244 double weight = 0.707106781;
11245 double B0u, B1u, B2u, B0v, B1v, B2v, Bx, By, Bz;
11246 double B00, B10, B20, B01, B11, B21, B02, B12, B22;
11250 B1u =
B_1(u) * weight;
11254 B1v =
B_1(v) * weight;
11269 sum_B = B00 + B10 + B20 + B01 + B11 + B21 + B02 + B12 + B22;
11271 Bx = B00 * ellipsoid -> pnt00.x + B10 * ellipsoid -> pnt10.x + B20 * ellipsoid -> pnt20.x
11272 + B01 * ellipsoid -> pnt01.x + B11 * ellipsoid -> pnt11.x + B21 * ellipsoid -> pnt21.x
11273 + B02 * ellipsoid -> pnt02.x + B12 * ellipsoid -> pnt12.x + B22 * ellipsoid -> pnt22.x;
11274 By = B00 * ellipsoid -> pnt00.y + B10 * ellipsoid -> pnt10.y + B20 * ellipsoid -> pnt20.y
11275 + B01 * ellipsoid -> pnt01.y + B11 * ellipsoid -> pnt11.y + B21 * ellipsoid -> pnt21.y
11276 + B02 * ellipsoid -> pnt02.y + B12 * ellipsoid -> pnt12.y + B22 * ellipsoid -> pnt22.y;
11277 Bz = B00 * ellipsoid -> pnt00.z + B10 * ellipsoid -> pnt10.z + B20 * ellipsoid -> pnt20.z
11278 + B01 * ellipsoid -> pnt01.z + B11 * ellipsoid -> pnt11.z + B21 * ellipsoid -> pnt21.z
11279 + B02 * ellipsoid -> pnt02.z + B12 * ellipsoid -> pnt12.z + B22 * ellipsoid -> pnt22.z;
11281 pnt -> x = Bx / sum_B;
11282 pnt -> y = By / sum_B;
11283 pnt -> z = Bz / sum_B;
11290 double weight = 0.707106781;
11291 double B0u, B1u, B2u, B0v, B1v, B2v, dB0u, dB1u, dB2u;
11292 double B00, B10, B20, B01, B11, B21, B02, B12, B22;
11293 double dB00, dB10, dB20, dB01, dB11, dB21, dB02, dB12, dB22;
11294 double sum_B, sum_dB, sum_B2, Bx, By, Bz, dBx, dBy, dBz;
11297 B1u =
B_1(u) * weight;
11301 B1v =
B_1(v) * weight;
11305 dB1u =
dB_1(u) * weight;
11332 sum_B = B00 + B10 + B20 + B01 + B11 + B21 + B02 + B12 + B22;
11333 sum_dB = dB00 + dB10 + dB20 + dB01 + dB11 + dB21 + dB02 + dB12 + dB22;
11335 sum_B2 = sum_B * sum_B;
11337 Bx = B00 * ellipsoid -> pnt00.x + B10 * ellipsoid -> pnt10.x + B20 * ellipsoid -> pnt20.x
11338 + B01 * ellipsoid -> pnt01.x + B11 * ellipsoid -> pnt11.x + B21 * ellipsoid -> pnt21.x
11339 + B02 * ellipsoid -> pnt02.x + B12 * ellipsoid -> pnt12.x + B22 * ellipsoid -> pnt22.x;
11340 By = B00 * ellipsoid -> pnt00.y + B10 * ellipsoid -> pnt10.y + B20 * ellipsoid -> pnt20.y
11341 + B01 * ellipsoid -> pnt01.y + B11 * ellipsoid -> pnt11.y + B21 * ellipsoid -> pnt21.y
11342 + B02 * ellipsoid -> pnt02.y + B12 * ellipsoid -> pnt12.y + B22 * ellipsoid -> pnt22.y;
11343 Bz = B00 * ellipsoid -> pnt00.z + B10 * ellipsoid -> pnt10.z + B20 * ellipsoid -> pnt20.z
11344 + B01 * ellipsoid -> pnt01.z + B11 * ellipsoid -> pnt11.z + B21 * ellipsoid -> pnt21.z
11345 + B02 * ellipsoid -> pnt02.z + B12 * ellipsoid -> pnt12.z + B22 * ellipsoid -> pnt22.z;
11347 dBx = dB00 * ellipsoid -> pnt00.x + dB10 * ellipsoid -> pnt10.x + dB20 * ellipsoid -> pnt20.x
11348 + dB01 * ellipsoid -> pnt01.x + dB11 * ellipsoid -> pnt11.x + dB21 * ellipsoid -> pnt21.x
11349 + dB02 * ellipsoid -> pnt02.x + dB12 * ellipsoid -> pnt12.x + dB22 * ellipsoid -> pnt22.x;
11350 dBy = dB00 * ellipsoid -> pnt00.y + dB10 * ellipsoid -> pnt10.y + dB20 * ellipsoid -> pnt20.y
11351 + dB01 * ellipsoid -> pnt01.y + dB11 * ellipsoid -> pnt11.y + dB21 * ellipsoid -> pnt21.y
11352 + dB02 * ellipsoid -> pnt02.y + dB12 * ellipsoid -> pnt12.y + dB22 * ellipsoid -> pnt22.y;
11353 dBz = dB00 * ellipsoid -> pnt00.z + dB10 * ellipsoid -> pnt10.z + dB20 * ellipsoid -> pnt20.z
11354 + dB01 * ellipsoid -> pnt01.z + dB11 * ellipsoid -> pnt11.z + dB21 * ellipsoid -> pnt21.z
11355 + dB02 * ellipsoid -> pnt02.z + dB12 * ellipsoid -> pnt12.z + dB22 * ellipsoid -> pnt22.z;
11357 grad -> x = (dBx * sum_B - Bx * sum_dB) / sum_B2;
11358 grad -> y = (dBy * sum_B - By * sum_dB) / sum_B2;
11359 grad -> z = (dBz * sum_B - Bz * sum_dB) / sum_B2;
11368 double weight = 0.707106781;
11369 double B0u, B1u, B2u, B0v, B1v, B2v, dB0v, dB1v, dB2v;
11370 double B00, B10, B20, B01, B11, B21, B02, B12, B22;
11371 double dB00, dB10, dB20, dB01, dB11, dB21, dB02, dB12, dB22;
11372 double sum_B, sum_dB, sum_B2, Bx, By, Bz, dBx, dBy, dBz;
11375 B1u =
B_1(u) * weight;
11379 B1v =
B_1(v) * weight;
11383 dB1v =
dB_1(v) * weight;
11410 sum_B = B00 + B10 + B20 + B01 + B11 + B21 + B02 + B12 + B22;
11411 sum_dB = dB00 + dB10 + dB20 + dB01 + dB11 + dB21 + dB02 + dB12 + dB22;
11413 sum_B2 = sum_B * sum_B;
11415 Bx = B00 * ellipsoid -> pnt00.x + B10 * ellipsoid -> pnt10.x + B20 * ellipsoid -> pnt20.x
11416 + B01 * ellipsoid -> pnt01.x + B11 * ellipsoid -> pnt11.x + B21 * ellipsoid -> pnt21.x
11417 + B02 * ellipsoid -> pnt02.x + B12 * ellipsoid -> pnt12.x + B22 * ellipsoid -> pnt22.x;
11418 By = B00 * ellipsoid -> pnt00.y + B10 * ellipsoid -> pnt10.y + B20 * ellipsoid -> pnt20.y
11419 + B01 * ellipsoid -> pnt01.y + B11 * ellipsoid -> pnt11.y + B21 * ellipsoid -> pnt21.y
11420 + B02 * ellipsoid -> pnt02.y + B12 * ellipsoid -> pnt12.y + B22 * ellipsoid -> pnt22.y;
11421 Bz = B00 * ellipsoid -> pnt00.z + B10 * ellipsoid -> pnt10.z + B20 * ellipsoid -> pnt20.z
11422 + B01 * ellipsoid -> pnt01.z + B11 * ellipsoid -> pnt11.z + B21 * ellipsoid -> pnt21.z
11423 + B02 * ellipsoid -> pnt02.z + B12 * ellipsoid -> pnt12.z + B22 * ellipsoid -> pnt22.z;
11425 dBx = dB00 * ellipsoid -> pnt00.x + dB10 * ellipsoid -> pnt10.x + dB20 * ellipsoid -> pnt20.x
11426 + dB01 * ellipsoid -> pnt01.x + dB11 * ellipsoid -> pnt11.x + dB21 * ellipsoid -> pnt21.x
11427 + dB02 * ellipsoid -> pnt02.x + dB12 * ellipsoid -> pnt12.x + dB22 * ellipsoid -> pnt22.x;
11428 dBy = dB00 * ellipsoid -> pnt00.y + dB10 * ellipsoid -> pnt10.y + dB20 * ellipsoid -> pnt20.y
11429 + dB01 * ellipsoid -> pnt01.y + dB11 * ellipsoid -> pnt11.y + dB21 * ellipsoid -> pnt21.y
11430 + dB02 * ellipsoid -> pnt02.y + dB12 * ellipsoid -> pnt12.y + dB22 * ellipsoid -> pnt22.y;
11431 dBz = dB00 * ellipsoid -> pnt00.z + dB10 * ellipsoid -> pnt10.z + dB20 * ellipsoid -> pnt20.z
11432 + dB01 * ellipsoid -> pnt01.z + dB11 * ellipsoid -> pnt11.z + dB21 * ellipsoid -> pnt21.z
11433 + dB02 * ellipsoid -> pnt02.z + dB12 * ellipsoid -> pnt12.z + dB22 * ellipsoid -> pnt22.z;
11435 grad -> x = (dBx * sum_B - Bx * sum_dB) / sum_B2;
11436 grad -> y = (dBy * sum_B - By * sum_dB) / sum_B2;
11437 grad -> z = (dBz * sum_B - Bz * sum_dB) / sum_B2;
12333 delta.
x = global -> x - origin -> x;
12334 delta.
y = global -> y - origin -> y;
12335 delta.
z = global -> z - origin -> z;
12344 global -> x =
dot_product(*local, trans -> x) + origin -> x;
12345 global -> y =
dot_product(*local, trans -> y) + origin -> y;
12346 global -> z =
dot_product(*local, trans -> z) + origin -> z;
12353 swap(trans -> x.
y, trans -> y.
x, tmp);
12354 swap(trans -> x.
z, trans -> z.
x, tmp);
12355 swap(trans -> y.
z, trans -> z.
y, tmp);
12402 va_start(args, format);
12403 vsprintf(buffer, format, args);
12407 fprintf(stderr,
"\nError: %s\n\n", buffer);
12422 fprintf(stderr,
"Warning: %s\n", buffer);
12424 fprintf(stderr,
"%s\n", buffer);
void check_ellipsoid_rec_consistency(meso3d::ellipsoid_rec &L)
ADDED for muMech needs.
#define dist_point(PNT1, PNT2)
#define dot_product(VEC1, VEC2)
int check_ellipsoid_ellipsoid_overlap(ellipsoid_rec *ellipsoid1, ellipsoid_rec *ellipsoid2)
static void error_message(int exit_code, const char *format,...)
void transfom_ellipsoid_rec(meso3d::ellipsoid_rec &L)
ADDED for muMech needs.
static void get_ellipsoid_point(ellipsoid_rec *ellipsoid, double u, double v, point_rec *pnt)
static bool check_point_inside_ellipsoid(point_rec *point, ellipsoid_rec *ellipsoid)
(ellipsoid greater by epsilon)
static void invert_transformation(trans_matrix *trans)
static void transform_from_local_to_global(point_rec *global, point_rec *local, point_rec *origin, trans_matrix *trans)
point_rec middle
Semiaxes. (must be righthanded Cartesian coordinate system)
bool isZero(double zero, double a)
static void get_ellipsoid_point_octant(ellipsoid_rec *ellipsoid, int octant, double u, double v, point_rec *pnt)
static void transform_from_global_to_local(point_rec *global, point_rec *local, point_rec *origin, trans_matrix *trans)
static int project_point_to_ellipsoid(point_rec *point, ellipsoid_rec *ellipsoid, point_rec *pnt, double *u_par, double *v_par)
(exact)
static double get_ellipsoid_v_gradient(ellipsoid_rec *ellipsoid, double u, double v, point_rec *v_grad)
#define div_vec(VEC, VALUE)
#define swap(VAL1, VAL2, TMP)
MESO3d - library of functions for geometry analysis of ellipsoids and preprocesor for T3d...
#define copy_vec(DEST, SRC)
#define sub_vec(RES, VEC1, VEC2)
#define cross_product(RES, VEC1, VEC2)
static double get_ellipsoid_u_gradient(ellipsoid_rec *ellipsoid, double u, double v, point_rec *u_grad)