1 #include "euler.h"
2
3 /* Euler transformation functions */
4
5 PG_FUNCTION_INFO_V1(spheretrans_in);
6 PG_FUNCTION_INFO_V1(spheretrans_from_float8);
7 PG_FUNCTION_INFO_V1(spheretrans_from_float8_and_type);
8 PG_FUNCTION_INFO_V1(spheretrans_equal);
9 PG_FUNCTION_INFO_V1(spheretrans_not_equal);
10 PG_FUNCTION_INFO_V1(spheretrans_phi);
11 PG_FUNCTION_INFO_V1(spheretrans_theta);
12 PG_FUNCTION_INFO_V1(spheretrans_psi);
13 PG_FUNCTION_INFO_V1(spheretrans_type);
14 PG_FUNCTION_INFO_V1(spheretrans_zxz);
15 PG_FUNCTION_INFO_V1(spheretrans);
16 PG_FUNCTION_INFO_V1(spheretrans_invert);
17 PG_FUNCTION_INFO_V1(spheretrans_trans);
18 PG_FUNCTION_INFO_V1(spheretrans_trans_inv);
19 PG_FUNCTION_INFO_V1(spheretrans_point);
20 PG_FUNCTION_INFO_V1(spheretrans_point_inverse);
21
22 /*
23 * Checks and modifies the Euler transformation.
24 */
25 static void
spheretrans_check(SEuler * e)26 spheretrans_check(SEuler *e)
27 {
28 SPoint sp[3];
29
30 sp[0].lat = sp[1].lat = sp[2].lat = 0.0;
31 sp[0].lng = e->phi;
32 sp[1].lng = e->theta;
33 sp[2].lng = e->psi;
34 spoint_check(&sp[0]);
35 spoint_check(&sp[1]);
36 spoint_check(&sp[2]);
37 e->phi = sp[0].lng;
38 e->theta = sp[1].lng;
39 e->psi = sp[2].lng;
40 }
41
42 Datum
spheretrans_in(PG_FUNCTION_ARGS)43 spheretrans_in(PG_FUNCTION_ARGS)
44 {
45 SEuler *se = (SEuler *) palloc(sizeof(SEuler));
46 char *c = PG_GETARG_CSTRING(0);
47 unsigned char etype[3];
48 int i;
49
50 void sphere_yyparse(void);
51
52 init_buffer(c);
53 sphere_yyparse();
54
55 if (get_euler(&se->phi, &se->theta, &se->psi, etype))
56 {
57
58 for (i = 0; i < 3; i++)
59 {
60 switch (i)
61 {
62 case 0:
63 se->phi_a = etype[i];
64 break;
65 case 1:
66 se->theta_a = etype[i];
67 break;
68 case 2:
69 se->psi_a = etype[i];
70 break;
71 }
72 }
73 spheretrans_check(se);
74 }
75 else
76 {
77 reset_buffer();
78 pfree(se);
79 se = NULL;
80 elog(ERROR, "spheretrans_in: parse error");
81 }
82 reset_buffer();
83 PG_RETURN_POINTER(se);
84 }
85
86 Datum
spheretrans_from_float8(PG_FUNCTION_ARGS)87 spheretrans_from_float8(PG_FUNCTION_ARGS)
88 {
89 SEuler *se = (SEuler *) palloc(sizeof(SEuler));
90
91 se->phi = PG_GETARG_FLOAT8(0);
92 se->theta = PG_GETARG_FLOAT8(1);
93 se->psi = PG_GETARG_FLOAT8(2);
94 seuler_set_zxz(se);
95 spheretrans_check(se);
96 PG_RETURN_POINTER(se);
97 }
98
99 Datum
spheretrans_from_float8_and_type(PG_FUNCTION_ARGS)100 spheretrans_from_float8_and_type(PG_FUNCTION_ARGS)
101 {
102 SEuler *se;
103 Datum d[3];
104 int i;
105 char *c = PG_GETARG_CSTRING(3);
106 unsigned char t = 0;
107
108 d[0] = PG_GETARG_DATUM(0);
109 d[1] = PG_GETARG_DATUM(1);
110 d[2] = PG_GETARG_DATUM(2);
111 se = (SEuler *) DatumGetPointer(
112 DirectFunctionCall3(spheretrans_from_float8,
113 d[0], d[1], d[2]));
114
115 for (i = 0; i < 3; i++)
116 {
117 switch (c[i])
118 {
119 case 'x':
120 case 'X':
121 t = EULER_AXIS_X;
122 break;
123 case 'y':
124 case 'Y':
125 t = EULER_AXIS_Y;
126 break;
127 case 'z':
128 case 'Z':
129 t = EULER_AXIS_Z;
130 break;
131 default:
132 t = 0;
133 }
134
135 if (t == 0)
136 {
137 pfree(se);
138 elog(ERROR, "invalid axis format");
139 }
140 switch (i)
141 {
142 case 0:
143 se->phi_a = t;
144 break;
145 case 1:
146 se->theta_a = t;
147 break;
148 case 2:
149 se->psi_a = t;
150 break;
151 }
152 }
153 PG_RETURN_POINTER(se);
154
155 }
156
157 void
seuler_set_zxz(SEuler * se)158 seuler_set_zxz(SEuler *se)
159 {
160 se->phi_a = EULER_AXIS_Z;
161 se->theta_a = EULER_AXIS_X;
162 se->psi_a = EULER_AXIS_Z;
163 }
164
165 bool
strans_eq(const SEuler * e1,const SEuler * e2)166 strans_eq(const SEuler *e1, const SEuler *e2)
167 {
168 SPoint in[2], p[4];
169
170 in[0].lng = 0.0;
171 in[0].lat = 0.0;
172 in[1].lng = PIH;
173 in[1].lat = 0.0;
174
175 euler_spoint_trans(&p[0], &in[0], e1);
176 euler_spoint_trans(&p[1], &in[1], e1);
177 euler_spoint_trans(&p[2], &in[0], e2);
178 euler_spoint_trans(&p[3], &in[1], e2);
179
180 return spoint_eq(&p[0], &p[2]) && spoint_eq(&p[1], &p[3]);
181 }
182
183 Datum
spheretrans_equal(PG_FUNCTION_ARGS)184 spheretrans_equal(PG_FUNCTION_ARGS)
185 {
186 SEuler *e1 = (SEuler *) PG_GETARG_POINTER(0);
187 SEuler *e2 = (SEuler *) PG_GETARG_POINTER(1);
188
189 PG_RETURN_BOOL(strans_eq(e1, e2));
190 }
191
192 Datum
spheretrans_not_equal(PG_FUNCTION_ARGS)193 spheretrans_not_equal(PG_FUNCTION_ARGS)
194 {
195 SEuler *e1 = (SEuler *) PG_GETARG_POINTER(0);
196 SEuler *e2 = (SEuler *) PG_GETARG_POINTER(1);
197
198 PG_RETURN_BOOL(!strans_eq(e1, e2));
199 }
200
201 Datum
spheretrans_phi(PG_FUNCTION_ARGS)202 spheretrans_phi(PG_FUNCTION_ARGS)
203 {
204 SEuler *se = (SEuler *) PG_GETARG_POINTER(0);
205
206 PG_RETURN_POINTER(&se->phi);
207 }
208
209 Datum
spheretrans_theta(PG_FUNCTION_ARGS)210 spheretrans_theta(PG_FUNCTION_ARGS)
211 {
212 SEuler *se = (SEuler *) PG_GETARG_POINTER(0);
213
214 PG_RETURN_POINTER(&se->theta);
215 }
216
217 Datum
spheretrans_psi(PG_FUNCTION_ARGS)218 spheretrans_psi(PG_FUNCTION_ARGS)
219 {
220 SEuler *se = (SEuler *) PG_GETARG_POINTER(0);
221
222 PG_RETURN_POINTER(&se->psi);
223 }
224
225 Datum
spheretrans_type(PG_FUNCTION_ARGS)226 spheretrans_type(PG_FUNCTION_ARGS)
227 {
228 SEuler *se = (SEuler *) PG_GETARG_POINTER(0);
229 BpChar *result = (BpChar *) palloc(3 + VARHDRSZ);
230 char ret[4];
231 int i;
232 unsigned char t = 0;
233
234 for (i = 0; i < 3; i++)
235 {
236 switch (i)
237 {
238 case 0:
239 t = se->phi_a;
240 break;
241 case 1:
242 t = se->theta_a;
243 break;
244 case 2:
245 t = se->psi_a;
246 break;
247 }
248 switch (t)
249 {
250 case EULER_AXIS_X:
251 ret[i] = 'X';
252 break;
253 case EULER_AXIS_Y:
254 ret[i] = 'Y';
255 break;
256 case EULER_AXIS_Z:
257 ret[i] = 'Z';
258 break;
259 }
260 }
261 ret[3] = '\0';
262 SET_VARSIZE(result, 3 + VARHDRSZ);
263 memcpy((void *) VARDATA(result), (void *) &ret[0], 3);
264
265 PG_RETURN_BPCHAR_P(result);
266 }
267
268 void
spheretrans_inv(SEuler * se)269 spheretrans_inv(SEuler *se)
270 {
271 float8 lng[3];
272 const unsigned char c = se->phi_a;
273
274 lng[2] = -se->phi;
275 lng[1] = -se->theta;
276 lng[0] = -se->psi;
277 se->phi = lng[0];
278 se->theta = lng[1];
279 se->psi = lng[2];
280 se->phi_a = se->psi_a;
281 se->psi_a = c;
282 }
283
284 void
spheretrans_inverse(SEuler * se_out,const SEuler * se_in)285 spheretrans_inverse(SEuler *se_out, const SEuler *se_in)
286 {
287 memcpy((void *) se_out, (void *) se_in, sizeof(SEuler));
288 spheretrans_inv(se_out);
289 }
290
291 void
strans_zxz(SEuler * ret,const SEuler * se)292 strans_zxz(SEuler *ret, const SEuler *se)
293 {
294 if (se->phi_a == EULER_AXIS_Z &&
295 se->theta_a == EULER_AXIS_X &&
296 se->psi_a == EULER_AXIS_Z)
297 {
298 memcpy((void *) ret, (void *) se, sizeof(SEuler));
299 }
300 else
301 {
302 SEuler tmp;
303
304 tmp.psi = 0.0;
305 tmp.theta = 0.0;
306 tmp.phi = 0.0;
307 seuler_set_zxz(&tmp);
308 seuler_trans_zxz(ret, se, &tmp);
309 }
310 }
311
312 Datum
spheretrans_zxz(PG_FUNCTION_ARGS)313 spheretrans_zxz(PG_FUNCTION_ARGS)
314 {
315 SEuler *si = (SEuler *) PG_GETARG_POINTER(0);
316 SEuler *ret = (SEuler *) palloc(sizeof(SEuler));
317
318 strans_zxz(ret, si);
319 PG_RETURN_POINTER(ret);
320 }
321
322 Datum
spheretrans(PG_FUNCTION_ARGS)323 spheretrans(PG_FUNCTION_ARGS)
324 {
325 Datum d = PG_GETARG_DATUM(0);
326
327 PG_RETURN_DATUM(d);
328 }
329
330 Datum
spheretrans_invert(PG_FUNCTION_ARGS)331 spheretrans_invert(PG_FUNCTION_ARGS)
332 {
333 SEuler *se = (SEuler *) PG_GETARG_POINTER(0);
334 SEuler *ret = (SEuler *) palloc(sizeof(SEuler));
335
336 spheretrans_inverse(ret, se);
337 PG_RETURN_POINTER(ret);
338 }
339
340 void
seuler_trans_zxz(SEuler * out,const SEuler * in,const SEuler * se)341 seuler_trans_zxz(SEuler *out, const SEuler *in, const SEuler *se)
342 {
343 SPoint sp[4];
344
345 sp[0].lng = 0.0;
346 sp[0].lat = 0.0;
347 sp[1].lng = PIH;
348 sp[1].lat = 0.0;
349 euler_spoint_trans(&sp[2], &sp[0], in);
350 euler_spoint_trans(&sp[3], &sp[1], in);
351 euler_spoint_trans(&sp[0], &sp[2], se);
352 euler_spoint_trans(&sp[1], &sp[3], se);
353 spherevector_to_euler(out, &sp[0], &sp[1]);
354 }
355
356 Datum
spheretrans_trans(PG_FUNCTION_ARGS)357 spheretrans_trans(PG_FUNCTION_ARGS)
358 {
359 SEuler *se1 = (SEuler *) PG_GETARG_POINTER(0);
360 SEuler *se2 = (SEuler *) PG_GETARG_POINTER(1);
361 SEuler *out = (SEuler *) palloc(sizeof(SEuler));
362
363 seuler_trans_zxz(out, se1, se2);
364 PG_RETURN_POINTER(out);
365 }
366
367 Datum
spheretrans_trans_inv(PG_FUNCTION_ARGS)368 spheretrans_trans_inv(PG_FUNCTION_ARGS)
369 {
370 SEuler *se1 = (SEuler *) PG_GETARG_POINTER(0);
371 SEuler *se2 = (SEuler *) PG_GETARG_POINTER(1);
372 SEuler *out = (SEuler *) palloc(sizeof(SEuler));
373 SEuler tmp;
374
375 spheretrans_inverse(&tmp, se2);
376 seuler_trans_zxz(out, se1, &tmp);
377 spheretrans_check(out);
378 PG_RETURN_POINTER(out);
379 }
380
381 void
euler_spoint_trans(SPoint * out,const SPoint * in,const SEuler * se)382 euler_spoint_trans(SPoint *out, const SPoint *in, const SEuler *se)
383 {
384 Vector3D v, o;
385
386 spoint_vector3d(&v, in);
387 euler_vector_trans(&o, &v, se);
388 vector3d_spoint(out, &o);
389 }
390
391 Datum
spheretrans_point(PG_FUNCTION_ARGS)392 spheretrans_point(PG_FUNCTION_ARGS)
393 {
394 SPoint *sp = (SPoint *) PG_GETARG_POINTER(0);
395 SEuler *se = (SEuler *) PG_GETARG_POINTER(1);
396 SPoint *out = (SPoint *) palloc(sizeof(SPoint));
397
398 euler_spoint_trans(out, sp, se);
399 PG_RETURN_POINTER(out);
400 }
401
402 Datum
spheretrans_point_inverse(PG_FUNCTION_ARGS)403 spheretrans_point_inverse(PG_FUNCTION_ARGS)
404 {
405 Datum sp = PG_GETARG_DATUM(0);
406 SEuler *se = (SEuler *) PG_GETARG_POINTER(1);
407 SEuler tmp;
408 Datum ret;
409
410 spheretrans_inverse(&tmp, se);
411 ret = DirectFunctionCall2(spheretrans_point,
412 sp, PointerGetDatum(&tmp));
413 PG_RETURN_DATUM(ret);
414 }
415
416 /*
417 * Transforms a spherical vector from 'spb' to 'spe' into an inverse Euler
418 * transformation. Returns true if the transformation was successful.
419 */
420 static bool
spherevector_to_euler_inv(SEuler * se,const SPoint * spb,const SPoint * spe)421 spherevector_to_euler_inv(SEuler *se, const SPoint *spb, const SPoint *spe)
422 {
423 if (spoint_eq(spb, spe))
424 {
425 return false;
426 }
427 else
428 {
429 Vector3D vbeg, vend, vtmp;
430 SPoint spt[2];
431 SEuler set;
432
433 spoint_vector3d(&vbeg, spb);
434 spoint_vector3d(&vend, spe);
435 vector3d_cross(&vtmp, &vbeg, &vend);
436 vector3d_spoint(&spt[0], &vtmp);
437 set.phi = -spt[0].lng - PIH;
438 set.theta = spt[0].lat - PIH;
439 set.psi = 0.0;
440 seuler_set_zxz(&set);
441 euler_spoint_trans(&spt[1], spb, &set);
442 set.psi = -spt[1].lng;
443 memcpy((void *) se, (void *) &set, sizeof(SEuler));
444 }
445
446 return true;
447 }
448
449 bool
spherevector_to_euler(SEuler * se,const SPoint * spb,const SPoint * spe)450 spherevector_to_euler(SEuler *se, const SPoint *spb, const SPoint *spe)
451 {
452 bool ret;
453
454 ret = spherevector_to_euler_inv(se, spb, spe);
455 if (ret)
456 {
457 spheretrans_inv(se);
458 }
459 return ret;
460 }
461
462 void
euler_vector_trans(Vector3D * out,const Vector3D * in,const SEuler * se)463 euler_vector_trans(Vector3D *out, const Vector3D *in, const SEuler *se)
464 {
465 int i;
466 unsigned char t;
467 const double *a;
468 double u[3], vr[3], sa, ca;
469
470 t = 0;
471 a = NULL;
472 u[0] = in->x;
473 u[1] = in->y;
474 u[2] = in->z;
475
476 for (i = 0; i < 3; i++)
477 {
478
479 switch (i)
480 {
481 case 0:
482 a = &se->phi;
483 t = se->phi_a;
484 break;
485 case 1:
486 a = &se->theta;
487 t = se->theta_a;
488 break;
489 case 2:
490 a = &se->psi;
491 t = se->psi_a;
492 break;
493 }
494
495 if (FPzero(*a))
496 {
497 continue;
498 }
499
500 sa = sin(*a);
501 ca = cos(*a);
502
503 switch (t)
504 {
505
506 case EULER_AXIS_X:
507 vr[0] = u[0];
508 vr[1] = ca * u[1] - sa * u[2];
509 vr[2] = sa * u[1] + ca * u[2];
510 break;
511 case EULER_AXIS_Y:
512 vr[0] = ca * u[0] + sa * u[2];
513 vr[1] = u[1];
514 vr[2] = ca * u[2] - sa * u[0];
515 break;
516 case EULER_AXIS_Z:
517 vr[0] = ca * u[0] - sa * u[1];
518 vr[1] = sa * u[0] + ca * u[1];
519 vr[2] = u[2];
520 break;
521
522 }
523 memcpy((void *) &u[0], (void *) &vr[0], sizeof(u));
524
525 }
526 out->x = u[0];
527 out->y = u[1];
528 out->z = u[2];
529 }
530