1
2
3 /*
4 A* -------------------------------------------------------------------
5 B* This file contains source code for the PyMOL computer program
6 C* Copyright (c) Schrodinger, LLC.
7 D* -------------------------------------------------------------------
8 E* It is unlawful to modify or remove this copyright notice.
9 F* -------------------------------------------------------------------
10 G* Please see the accompanying LICENSE file for further information.
11 H* -------------------------------------------------------------------
12 I* Additional authors of this source file include:
13 -*
14 -*
15 -*
16 Z* -------------------------------------------------------------------
17 */
18 #ifndef _H_Vector
19 #define _H_Vector
20
21 #include"os_predef.h"
22 #include"os_gl.h"
23 #include<math.h>
24
25 /* NOTE THIS VERSION USES RADIANS BY DEFAULT! */
26
27
28 /* NOTE: Matrices are assumed to be row-major (C-like not
29 * OpenGL-like) unless explictly labeled as per the following
30 * conventions:
31 *
32 * row-major: 33f, 33d, 44f, 44d, R33f, R33d, R44f, R44d
33 * column-major: C33f, C33d, C44f, C44d
34 */
35
36 #define cPI 3.14159265358979323846 /* pi */
37
38 #define GET_BIT(val, bit) (((val) >> (bit)) & 1)
39 #define SET_BIT(val, bit) (val) |= (1 << (bit))
40 #define SET_BIT_OFF(val, bit) (val) &= ~(1 << (bit))
41 #define SET_BIT_TO(val, bit, v) {if(v) SET_BIT(val, bit); else SET_BIT_OFF(val, bit);}
42
43 short countBits(unsigned long bits);
44 short countBitsInt(int bits);
45
46 typedef float Vector3f[3]; /* for local vars only - use float* for parameters */
47 typedef float Vector4f[4];
48 typedef int Vector3i[3];
49
50 typedef float Matrix33f[3][3];
51 typedef double Matrix33d[3][3];
52 typedef float Matrix53f[5][3];
53 typedef double Matrix53d[5][3];
54
55 unsigned int optimizer_workaround1u(unsigned int value);
56
57 float get_random0to1f(void);
58
59 float deg_to_rad(float angle);
60 float rad_to_deg(float angle);
61
62 void normalize23f(const float *v1, float *v2);
63 void normalize3d(double *v1);
64 void normalize2f(float *v1);
65 void normalize4f(float *v1);
66
67 void clamp3f(float *v1);
68 void get_divergent3f(const float *src, float *dst);
69 void get_random3f(float *x);
70 void scatter3f(float *v, float weight);
71 void wiggle3f(float *v, const float *p, const float *s);
72 void extrapolate3f(const float *v1, const float *unit, float *result);
73
74 void mix3f(const float *v1, const float *v2, float fxn, float *v3);
75 void mix3d(const double *v1, const double *v2, double fxn, double *v3);
76
77 void get_system3f(float *x, float *y, float *z); /* make random system */
78 void get_system1f3f(float *x, float *y, float *z); /* make system in direction of x */
79 void get_system2f3f(float *x, float *y, float *z); /* make system in direction of x, perp to x,y */
80
81 double dot_product3d(const double *v1, const double *v2);
82 void remove_component3d(const double *v1, const double *unit, double *result);
83 void cross_product3d(const double *v1, const double *v2, double *cross);
84 void scale3d(const double *v1, const double v0, double *v2);
85 void add3d(const double *v1, const double *v0, double *v2);
86
87 double distance_line2point3f(const float *base, const float *normal, const float *point,
88 float *alongNormalSq);
89 double distance_halfline2point3f(const float *base, const float *normal, const float *point,
90 float *alongNormalSq);
91
92 int equal3f(const float *v1, const float *v2);
93
94 int pymol_roundf(float f);
95
96 float get_angle3f(const float *v1, const float *v2);
97 float get_dihedral3f(const float *v0, const float *v1, const float *v2, const float *v3);
98 double length3d(const double *v1);
99
100 void min3f(const float *v1, const float *v2, float *v3);
101 void max3f(const float *v1, const float *v2, float *v3);
102
103 void dump3i(const int *v, const char *prefix);
104 void dump2f(const float *v, const char *prefix);
105 void dump3f(const float *v, const char *prefix);
106 void dump3d(const double *v, const char *prefix);
107 void dump4f(const float *v, const char *prefix);
108 void dump33f(const float *m, const char *prefix);
109 void dump33d(const double *m, const char *prefix);
110 void dump44f(const float *m, const char *prefix);
111 void dump44d(const double *m, const char *prefix);
112
113 void copy44f(const float *src, float *dst);
114 void copy44d(const double *src, double *dst);
115
116 void identity33f(float *m1);
117 void identity33d(double *m);
118 void identity44f(float *m1);
119 void identity44d(double *m1);
120
121 bool is_identityf(int n, const float *m, float threshold=1.0E-6F);
122 bool is_allclosef(int nrow,
123 const float *m1, int ncol1,
124 const float *m2, int ncol2, float threshold=1.0E-6F);
125 bool is_diagonalf(int nrow,
126 const float *m, int ncol=0, float threshold=1.0E-6F);
127 double determinant33f(const float *m, int ncol=3);
128
129 void glOrtho44f(float *m1,
130 GLfloat left, GLfloat right,
131 GLfloat bottom, GLfloat top,
132 GLfloat nearVal, GLfloat farVal);
133 void glFrustum44f(float *m1,
134 GLfloat left, GLfloat right,
135 GLfloat bottom, GLfloat top,
136 GLfloat nearVal, GLfloat farVal);
137
138 void copy44f44f(const float *src, float *dst);
139 void copy44d44f(const double *src, float *dst);
140 void copy44f44d(const float *src, double *dst);
141
142 void copy44d33f(const double *src, float *dst);
143 void copy44f33f(const float *src, float *dst);
144 void copy33f44d(const float *src, double *dst);
145 void copy33f44f(const float *src, float *dst);
146 void copy3d3f(const double *v1, float *v2);
147 void copy3f3d(const float *v1, double *v2);
148
149
150 /* in the following matrix multiplies and transformations:
151 the last two matrices can be the same matrix! */
152
153 void transpose33f33f(const float *m1, float *m2);
154 void transpose33d33d(const double *m1, double *m2);
155 void transpose44f44f(const float *m1, float *m2);
156 void transpose44d44d(const double *m1, double *m2);
157
158 void transform33f3f(const float *m1, const float *m2, float *m3);
159 void transform33Tf3f(const float *m1, const float *m2, float *m3); /* uses transpose */
160
161 void transform44f3f(const float *m1, const float *m2, float *m3);
162 void transform44f4f(const float *m1, const float *m2, float *m3);
163
164 void transform44d3f(const double *m1, const float *m2, float *m3);
165 void transform44d3d(const double *m1, const double *m2, double *m3);
166 void inverse_transformC44f3f(const float *m1, const float *m2, float *m3);
167 void inverse_transform44f3f(const float *m1, const float *m2, float *m3);
168 void inverse_transform44d3f(const double *m1, const float *m2, float *m3);
169 void inverse_transform44d3d(const double *m1, const double *m2, double *m3);
170 void transform44f3fas33f3f(const float *m1, const float *m2, float *m3);
171 void transform44d3fas33d3f(const double *m1, const float *m2, float *m3);
172
173 void multiply33f33f(const float *m1, const float *m2, float *m3);
174 void multiply33d33d(const double *m1, const double *m2, double *m3);
175
176
177 /* as matrix types */
178
179 void matrix_transform33f3f(const Matrix33f m1, const float *v1, float *v2);
180 void matrix_inverse_transform33f3f(const Matrix33f m1, const float *v1, float *v2);
181
182 void rotation_to_matrix33f(const float *axis, float angle, Matrix33f mat);
183 void matrix_multiply33f33f(Matrix33f m1, Matrix33f m2, Matrix33f m3);
184 void matrix_multiply33d33d(Matrix33d m1, Matrix33d m2, Matrix33d m3);
185
186
187 /* A 4x4 TTT matrix is really a 3x3 rotation matrix with two translation vectors:
188 (1) a pre-translation stored in forth row, first three columns.
189 (2) and a post-translation stored in forth column, first three rows.
190 There are certain cases where this representation is more convenient.
191 */
192 void combineTTT44f44f(const float *m1, const float *m2, float *m3);
193 void transformTTT44f3f(const float *m1, const float *m2, float *m3);
194 void transform_normalTTT44f3f(const float *m1, const float *m2, float *m3);
195 void initializeTTT44f(float *m);
196
197 void multiply44d44d44d(const double *left, const double *right, double *product);
198 void left_multiply44d44d(const double *left, double *right);
199 void right_multiply44d44d(double *left, const double *right);
200
201 void multiply44f44f44f(const float *left, const float *right, float *product);
202 void left_multiply44f44f(const float *left, float *right);
203 void right_multiply44f44f(float *left, const float *right);
204
205 void reorient44d(double *matrix);
206
207 void recondition33d(double *matrix);
208 void recondition44d(double *matrix);
209
210
211 /* invert a 4x4 homogenous that contains just rotation & tranlation
212 (e.g. no scaling & fourth row is 0,0,0,1) */
213 void invert_special44d44d(const double *original, double *inv);
214 void invert_special44f44f(const float *original, float *inv);
215
216 void invert_rotation_only44d44d(const double *original, double *inv);
217
218 void convertTTTfR44d(const float *ttt, double *homo);
219 void convertTTTfR44f(const float *ttt, float *homo);
220 void convertR44dTTTf(const double *homo, float *ttt);
221 void convert44d44f(const double *dbl, float *flt);
222 void convert44f44d(const float *flt, double *dbl);
223
224 void get_rotation_about3f3fTTTf(float angle, const float *dir, const float *origin, float *ttt);
225
226
227 /* end revised matrix routines */
228
229
230 /*------------------------------------------------------------------------*/
231
232
233 /* OLD MATRIX STUFF below NEEDS REWORKING */
234
235 void rotation_matrix3f(float angle, float x, float y, float z, float *m);
236
237 typedef float *oMatrix5f[5]; /* PHASE THESE OUT! - THEY CAUSE PROBLEMS! */
238
239 typedef float *oMatrix3f[3];
240
241 typedef float *oMatrix3d[3];
242
243
244 /*void matcopy ( oMatrix5f to, oMatrix5f from );
245 void mattran ( oMatrix5f nm, oMatrix5f om, int axis, float dist );
246 void matrot ( oMatrix5f nm, oMatrix5f om, int axis, float angle );*/
247
248 void matrix_to_rotation(Matrix53f rot, float *axis, float *angle);
249 void rotation_to_matrix(Matrix53f rot, const float *axis, float angle);
250
251 void transform3d3f(const oMatrix3d m1, const float *v1, float *v2);
252 void transform33d3f(const Matrix33d m1, const float *v1, float *v2);
253 void transform5f3f(const oMatrix5f m, const float *v1, float *v2);
254
255 void mult4f(const float *vsrc, float val, float *vdest);
256 void mult3f(const float *vsrc, float val, float *vdest);
257 float max3(float val1, float val2, float val3);
258 float ave3(float val1, float val2, float val3);
259 float ave2(float val1, float val2);
260 void white4f(float *rgba, float value);
261 void add4f(const float *v1, const float *v2, float *sum);
262
263 int countchrs(const char *str, char ch);
264
265 float smooth(float x, float power);
266
267 void subdivide(int n, float *x, float *y);
268
269 //-------------------------------------------------------------------------
270 // Small inline functions
271 // (many of these were macros up to PyMOL 1.7.6)
272 //-------------------------------------------------------------------------
273
274 static const float _0f_inline = 0.0F;
275 static const double _0d_inline = 0.0;
276 static const float _1f_inline = 1.0F;
277 static const double _1d_inline = 1.0;
278 static const float R_SMALL_inline = 0.000000001F;
279 static const double R_SMALLd_inline = 0.000000001;
280
281 #define normalize3f inline_normalize3f
282 #define sqrt1f inline_sqrt1f
283 #define sqrt1d inline_sqrt1d
284 #define diff3f inline_diff3f
285 #define diffsq3f inline_diffsq3f
286 #define within3f inline_within3f
287 #define within3fsq inline_within3fsq
288 #define within3fret inline_within3fret
289 #define remove_component3f inline_remove_component3f
290 #define project3f inline_project3f
291
set3f(float * v1,float x,float y,float z)292 inline void set3f(float * v1, float x, float y, float z) {
293 v1[0] = x;
294 v1[1] = y;
295 v1[2] = z;
296 }
297
298 template <typename T>
zero3(T * v1)299 inline void zero3(T * v1) {
300 v1[0] = 0;
301 v1[1] = 0;
302 v1[2] = 0;
303 }
304
305 #define zero3f zero3
306 #define zero3i zero3
307
ones3f(float * v1)308 inline void ones3f(float * v1) {
309 v1[0] = 1.0F;
310 v1[1] = 1.0F;
311 v1[2] = 1.0F;
312 }
313
314 /**
315 * Takes the dot product of two three-term vectors
316 * @param v1 an array of three constant floats
317 * @param v2 an array of three constant floats
318 * @return the dot product
319 */
dot_product3f(const float * v1,const float * v2)320 inline float dot_product3f(const float * v1, const float * v2) {
321 return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
322 }
323
324 /**
325 * Sets the values of the second vector to the additive inverses of each
326 * respective value of the first vector, leaving the first vector unchanged.
327 * @param v1 an array of three constant floats
328 * @param v2 an array of three constant floats
329 */
invert3f3f(const float * v1,float * v2)330 inline void invert3f3f(const float * v1, float * v2) {
331 v2[0] = -v1[0];
332 v2[1] = -v1[1];
333 v2[2] = -v1[2];
334 }
335
invert3f(float * v)336 inline void invert3f(float * v) {
337 invert3f3f(v, v);
338 }
339
scale3f(const float * v1,float v0,float * v2)340 inline void scale3f(const float * v1, float v0, float * v2) {
341 v2[0] = v1[0] * v0;
342 v2[1] = v1[1] * v0;
343 v2[2] = v1[2] * v0;
344 }
345
346 template <typename S, typename D>
copyN(const S * src,D * dst,int N)347 void copyN(const S * src, D * dst, int N) {
348 for (int i = 0; i < N; ++i)
349 dst[i] = src[i];
350 }
351
352 template <typename S, typename D>
copy2(const S * src,D * dst)353 void copy2(const S * src, D * dst) {
354 dst[0] = src[0];
355 dst[1] = src[1];
356 }
357
358 template <typename S, typename D>
copy3(const S * src,D * dst)359 void copy3(const S * src, D * dst) {
360 dst[0] = src[0];
361 dst[1] = src[1];
362 dst[2] = src[2];
363 }
364
365 template <typename S, typename D>
copy4(const S * src,D * dst)366 void copy4(const S * src, D * dst) {
367 dst[0] = src[0];
368 dst[1] = src[1];
369 dst[2] = src[2];
370 dst[3] = src[3];
371 }
372
373 #define copy2f copy2<float, float>
374 #define copy3f copy3
375 #define copy3d copy3<double, double>
376 #define copy4f copy4
377
add3f(const float * v1,const float * v2,float * v3)378 inline void add3f(const float * v1, const float * v2, float * v3) {
379 v3[0] = v1[0] + v2[0];
380 v3[1] = v1[1] + v2[1];
381 v3[2] = v1[2] + v2[2];
382 }
383
subtract3f(const float * v1,const float * v2,float * v3)384 inline void subtract3f(const float * v1, const float * v2, float * v3) {
385 v3[0] = v1[0] - v2[0];
386 v3[1] = v1[1] - v2[1];
387 v3[2] = v1[2] - v2[2];
388 }
389
lengthsq3f(const float * v1)390 inline float lengthsq3f(const float * v1) {
391 return (v1[0] * v1[0]) + (v1[1] * v1[1]) + (v1[2] * v1[2]);
392 }
393
average3f(const float * v1,const float * v2,float * avg)394 inline void average3f(const float * v1, const float * v2, float * avg) {
395 (avg)[0] = ((v1)[0]+(v2)[0])/2;
396 (avg)[1] = ((v1)[1]+(v2)[1])/2;
397 (avg)[2] = ((v1)[2]+(v2)[2])/2;
398 }
399
cross_product3f(const float * v1,const float * v2,float * cross)400 inline void cross_product3f(const float * v1, const float * v2, float * cross) {
401 cross[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
402 cross[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
403 cross[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
404 }
405
inline_sqrt1f(float f)406 inline double inline_sqrt1f(float f)
407 { /* no good as a macro because f is used twice */
408 if(f > _0f_inline)
409 return (sqrt(f));
410 else
411 return (_0d_inline);
412 }
413
inline_sqrt1d(double f)414 inline double inline_sqrt1d(double f)
415 { /* no good as a macro because f is used twice */
416 if(f > _0d_inline)
417 return (sqrt(f));
418 else
419 return (_0d_inline);
420 }
421
length3f(const float * v1)422 inline float length3f(const float * v1) {
423 return sqrt1f(lengthsq3f(v1));
424 }
425
length2f(const float * v1)426 inline float length2f(const float * v1) {
427 return sqrt1f((v1[0] * v1[0]) + (v1[1] * v1[1]));
428 }
429
inline_normalize3f(float * v1)430 inline void inline_normalize3f(float *v1)
431 {
432 double vlen = length3f(v1);
433 if(vlen > R_SMALLd_inline) {
434 float inV = (float) (_1d_inline / vlen);
435 v1[0] *= inV;
436 v1[1] *= inV;
437 v1[2] *= inV;
438 } else {
439 v1[0] = v1[1] = v1[2] = _0f_inline;
440 }
441 }
442
inline_diff3f(const float * v1,const float * v2)443 inline double inline_diff3f(const float *v1, const float *v2)
444 {
445 float dx, dy, dz;
446 dx = (v1[0] - v2[0]);
447 dy = (v1[1] - v2[1]);
448 dz = (v1[2] - v2[2]);
449 return (sqrt1d(dx * dx + dy * dy + dz * dz));
450 }
451
inline_diffsq3f(const float * v1,const float * v2)452 inline float inline_diffsq3f(const float *v1, const float *v2)
453 {
454 float dx, dy, dz;
455 dx = (v1[0] - v2[0]);
456 dy = (v1[1] - v2[1]);
457 dx = dx * dx;
458 dz = (v1[2] - v2[2]);
459 dy = dy * dy;
460 return (dz * dz + (dx + dy));
461 }
462
inline_within3f(const float * v1,const float * v2,float dist)463 inline int inline_within3f(const float *v1, const float *v2, float dist)
464 {
465 float dx, dy, dz, dist2;
466 dx = (float) fabs(v1[0] - v2[0]);
467 dy = (float) fabs(v1[1] - v2[1]);
468 if(dx > dist)
469 return (0);
470 dz = (float) fabs(v1[2] - v2[2]);
471 dx = dx * dx;
472 if(dy > dist)
473 return (0);
474 dy = dy * dy;
475 dist2 = dist * dist;
476 if(dz > dist)
477 return (0);
478 return (((dx + dy) + dz * dz) <= dist2);
479 }
480
inline_within3fsq(const float * v1,const float * v2,float dist,float dist2)481 inline int inline_within3fsq(const float *v1, const float *v2, float dist, float dist2)
482 {
483 /* manually optimized to take advantage of parallel execution units */
484 float dx, dy, dz;
485 dx = v1[0] - v2[0];
486 dy = v1[1] - v2[1];
487 dx = (float) fabs(dx);
488 dy = (float) fabs(dy);
489 if(dx > dist)
490 return (0);
491 dz = v1[2] - v2[2];
492 dx = dx * dx;
493 if(dy > dist)
494 return (0);
495 dz = (float) fabs(dz);
496 dy = dy * dy;
497 if(dz > dist)
498 return (0);
499 dx = dx + dy;
500 dz = dz * dz;
501 if(dx > dist2)
502 return (0);
503 return ((dx + dz) <= (dist2));
504 }
505
inline_within3fret(const float * v1,const float * v2,float cutoff,const float cutoff2,float * diff,float * dist)506 inline int inline_within3fret(const float *v1, const float *v2, float cutoff,
507 const float cutoff2, float *diff, float *dist)
508 {
509 float dx, dy, dz, dist2;
510 dx = (float) fabs((diff[0] = v1[0] - v2[0]));
511 dy = (float) fabs((diff[1] = v1[1] - v2[1]));
512 if(dx > cutoff)
513 return 0;
514 dz = (float) fabs((diff[2] = v1[2] - v2[2]));
515 dx = dx * dx;
516 if(dy > cutoff)
517 return 0;
518 dy = dy * dy;
519 if(dz > cutoff)
520 return 0;
521 if((dist2 = ((dx + dy) + dz * dz)) > cutoff2)
522 return 0;
523 *dist = (float) sqrt1f(dist2);
524 return 1;
525 }
526
inline_remove_component3f(const float * v1,const float * unit,float * result)527 inline void inline_remove_component3f(const float *v1, const float *unit, float *result)
528 {
529 float dot;
530
531 dot = v1[0] * unit[0] + v1[1] * unit[1] + v1[2] * unit[2];
532 result[0] = v1[0] - unit[0] * dot;
533 result[1] = v1[1] - unit[1] * dot;
534 result[2] = v1[2] - unit[2] * dot;
535 }
536
inline_project3f(const float * v1,const float * v2,float * proj)537 inline float inline_project3f(const float *v1, const float *v2, float *proj)
538 {
539 float dot;
540
541 dot = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
542 proj[0] = v2[0] * dot;
543 proj[1] = v2[1] * dot;
544 proj[2] = v2[2] * dot;
545
546 return (dot);
547 }
548
549 #endif
550