Openholo  v1.1
Open Source Digital Holographic Library
vec.h
Go to the documentation of this file.
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install, copy or use the software.
7 //
8 //
9 // License Agreement
10 // For Open Source Digital Holographic Library
11 //
12 // Openholo library is free software;
13 // you can redistribute it and/or modify it under the terms of the BSD 2-Clause license.
14 //
15 // Copyright (C) 2017-2024, Korea Electronics Technology Institute. All rights reserved.
16 // E-mail : contact.openholo@gmail.com
17 // Web : http://www.openholo.org
18 //
19 // Redistribution and use in source and binary forms, with or without modification,
20 // are permitted provided that the following conditions are met:
21 //
22 // 1. Redistribution's of source code must retain the above copyright notice,
23 // this list of conditions and the following disclaimer.
24 //
25 // 2. Redistribution's in binary form must reproduce the above copyright notice,
26 // this list of conditions and the following disclaimer in the documentation
27 // and/or other materials provided with the distribution.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the copyright holder or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 // This software contains opensource software released under GNU Generic Public License,
41 // NVDIA Software License Agreement, or CUDA supplement to Software License Agreement.
42 // Check whether software you use contains licensed software.
43 //
44 //M*/
45 
46 #ifndef __vec_h
47 #define __vec_h
48 // Description:
49 // Mathematical tools to handle n-dimensional vectors
50 //
51 // Author:
52 // Myung-Joon Kim
53 // Dae-Hyun Kim
54 
55 
56 #include "ivec.h"
57 #include "epsilon.h"
58 #include <math.h>
59 #include <stdio.h>
60 
61 namespace oph {
62 
66  struct __declspec(dllexport) vec2 {
67  Real v[2];
68  static const int n;
69 
70  inline vec2() { }
71  inline vec2(Real a) { v[1 - 1] = a; v[2 - 1] = a; }
72  inline vec2(Real v_1, Real v_2) { v[1 - 1] = v_1; v[2 - 1] = v_2; }
73  inline vec2(const ivec2& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; }
74  inline vec2(const vec2& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; }
75 
76 
77  inline vec2& operator=(const vec2& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; return *this; }
78  inline Real& operator[] (int i) { return v[i]; }
79  inline const Real& operator[] (int i) const { return v[i]; }
80  inline Real& operator() (int i) { return v[i % 2]; }
81  inline const Real& operator() (int i) const { return v[i % 2]; }
82 
83  bool unit();
84  Real length() const;
85 
86  inline bool is_zero() const { return (v[0] == 0.0 && v[1] == 0.0); }
87  inline bool is_tiny(Real tiny_tol = epsilon) const {
88  return (fabs(v[0]) <= tiny_tol && fabs(v[1]) <= tiny_tol);
89  }
90 
91  //
92  // returns 1: this and other vectors are parallel
93  // -1: this and other vectors are anti-parallel
94  // 0: this and other vectors are not parallel
95  // or at least one of the vectors is zero
96  int is_parallel(
97  const vec2&, // other vector
98  Real = angle_tolerance // optional angle tolerance (radians)
99  ) const;
100 
101  // returns true: this and other vectors are perpendicular
102  // false: this and other vectors are not perpendicular
103  // or at least one of the vectors is zero
104  bool is_perpendicular(
105  const vec2&, // other vector
106  Real = angle_tolerance // optional angle tolerance (radians)
107  ) const;
108 
109  //
110  // set this vector to be perpendicular to another vector
111  bool perpendicular( // Result is not unitized.
112  // returns false if input vector is zero
113  const vec2&
114  );
115 
116  //
117  // set this vector to be perpendicular to a line defined by 2 points
118  bool perpendicular(
119  const vec2&,
120  const vec2&
121  );
122  };
123 
124 
125 
126 
127 
128  //| binary op : componentwise
129 
130 
131  inline vec2 operator + (const vec2& a, const vec2& b)
132  {
133  vec2 c;
134  for (int i = 0; i < 2; ++i) { c[i] = a[i] + b[i]; }
135  return c;
136  }
137 
138  inline vec2 operator + (Real a, const vec2& b)
139  {
140  vec2 c;
141  for (int i = 0; i < 2; ++i) { c[i] = a + b[i]; }
142  return c;
143  }
144 
145  inline vec2 operator + (const vec2& a, Real b)
146  {
147  vec2 c;
148  for (int i = 0; i < 2; ++i) { c[i] = a[i] + b; }
149  return c;
150  }
151 
152 
153 
154  inline vec2 operator - (const vec2& a, const vec2& b)
155  {
156  vec2 c;
157  for (int i = 0; i < 2; ++i) { c[i] = a[i] - b[i]; }
158  return c;
159  }
160 
161  inline vec2 operator - (Real a, const vec2& b)
162  {
163  vec2 c;
164  for (int i = 0; i < 2; ++i) { c[i] = a - b[i]; }
165  return c;
166  }
167 
168  inline vec2 operator - (const vec2& a, Real b)
169  {
170  vec2 c;
171  for (int i = 0; i < 2; ++i) { c[i] = a[i] - b; }
172  return c;
173  }
174 
175 
176 
177  inline vec2 operator * (const vec2& a, const vec2& b)
178  {
179  vec2 c;
180  for (int i = 0; i < 2; ++i) { c[i] = a[i] * b[i]; }
181  return c;
182  }
183 
184  inline vec2 operator * (Real a, const vec2& b)
185  {
186  vec2 c;
187  for (int i = 0; i < 2; ++i) { c[i] = a * b[i]; }
188  return c;
189  }
190 
191  inline vec2 operator * (const vec2& a, Real b)
192  {
193  vec2 c;
194  for (int i = 0; i < 2; ++i) { c[i] = a[i] * b; }
195  return c;
196  }
197 
198 
199 
200  inline vec2 operator / (const vec2& a, const vec2& b)
201  {
202  vec2 c;
203  for (int i = 0; i < 2; ++i) { c[i] = a[i] / b[i]; }
204  return c;
205  }
206 
207  inline vec2 operator / (Real a, const vec2& b)
208  {
209  vec2 c;
210  for (int i = 0; i < 2; ++i) { c[i] = a / b[i]; }
211  return c;
212  }
213 
214  inline vec2 operator / (const vec2& a, Real b)
215  {
216  vec2 c;
217  for (int i = 0; i < 2; ++i) { c[i] = a[i] / b; }
218  return c;
219  }
220 
221 
222 
223  //| cumulative op : componentwise
224 
225 
226  inline vec2 operator += (vec2& a, const vec2& b)
227  {
228  return a = (a + b);
229  }
230 
231  inline vec2 operator += (vec2& a, Real b)
232  {
233  return a = (a + b);
234  }
235 
236 
237 
238  inline vec2 operator -= (vec2& a, const vec2& b)
239  {
240  return a = (a - b);
241  }
242 
243  inline vec2 operator -= (vec2& a, Real b)
244  {
245  return a = (a - b);
246  }
247 
248 
249 
250  inline vec2 operator *= (vec2& a, const vec2& b)
251  {
252  return a = (a * b);
253  }
254 
255  inline vec2 operator *= (vec2& a, Real b)
256  {
257  return a = (a * b);
258  }
259 
260 
261 
262  inline vec2 operator /= (vec2& a, const vec2& b)
263  {
264  return a = (a / b);
265  }
266 
267  inline vec2 operator /= (vec2& a, Real b)
268  {
269  return a = (a / b);
270  }
271 
272 
273 
274  //| logical op : componentwise
275 
276 
277  inline int operator == (const vec2& a, const vec2& b)
278  {
279  int c = 1;
280  for (int i = 0; i < 2; ++i) { c = c && a[i] == b[i]; }
281  return c;
282  }
283 
284  inline int operator == (Real a, const vec2& b)
285  {
286  int c = 1;
287  for (int i = 0; i < 2; ++i) { c = c && a == b[i]; }
288  return c;
289  }
290 
291  inline int operator == (const vec2& a, Real b)
292  {
293  int c = 1;
294  for (int i = 0; i < 2; ++i) { c = c && a[i] == b; }
295  return c;
296  }
297 
298 
299 
300  inline int operator < (const vec2& a, const vec2& b)
301  {
302  int c = 1;
303  for (int i = 0; i < 2; ++i) { c = c && a[i] < b[i]; }
304  return c;
305  }
306 
307  inline int operator < (Real a, const vec2& b)
308  {
309  int c = 1;
310  for (int i = 0; i < 2; ++i) { c = c && a < b[i]; }
311  return c;
312  }
313 
314  inline int operator < (const vec2& a, Real b)
315  {
316  int c = 1;
317  for (int i = 0; i < 2; ++i) { c = c && a[i] < b; }
318  return c;
319  }
320 
321 
322 
323  inline int operator <= (const vec2& a, const vec2& b)
324  {
325  int c = 1;
326  for (int i = 0; i < 2; ++i) { c = c && a[i] <= b[i]; }
327  return c;
328  }
329 
330  inline int operator <= (Real a, const vec2& b)
331  {
332  int c = 1;
333  for (int i = 0; i < 2; ++i) { c = c && a <= b[i]; }
334  return c;
335  }
336 
337  inline int operator <= (const vec2& a, Real b)
338  {
339  int c = 1;
340  for (int i = 0; i < 2; ++i) { c = c && a[i] <= b; }
341  return c;
342  }
343 
344 
345 
346  inline int operator > (const vec2& a, const vec2& b)
347  {
348  int c = 1;
349  for (int i = 0; i < 2; ++i) { c = c && a[i] > b[i]; }
350  return c;
351  }
352 
353  inline int operator > (Real a, const vec2& b)
354  {
355  int c = 1;
356  for (int i = 0; i < 2; ++i) { c = c && a > b[i]; }
357  return c;
358  }
359 
360  inline int operator > (const vec2& a, Real b)
361  {
362  int c = 1;
363  for (int i = 0; i < 2; ++i) { c = c && a[i] > b; }
364  return c;
365  }
366 
367 
368 
369  inline int operator >= (const vec2& a, const vec2& b)
370  {
371  int c = 1;
372  for (int i = 0; i < 2; ++i) { c = c && a[i] >= b[i]; }
373  return c;
374  }
375 
376  inline int operator >= (Real a, const vec2& b)
377  {
378  int c = 1;
379  for (int i = 0; i < 2; ++i) { c = c && a >= b[i]; }
380  return c;
381  }
382 
383  inline int operator >= (const vec2& a, Real b)
384  {
385  int c = 1;
386  for (int i = 0; i < 2; ++i) { c = c && a[i] >= b; }
387  return c;
388  }
389 
390 
391 
392  //| unary op : componentwise
393  inline vec2 operator - (const vec2& a)
394  {
395  vec2 c;
396  for (int i = 0; i < 2; ++i) { c[i] = -a[i]; }
397  return c;
398  }
399 
400  //| R^n -> R
401  inline Real sum(const vec2& a)
402  {
403  Real s = 0;
404 
405  s += a[1 - 1];
406  s += a[2 - 1];
407 
408  return s;
409  }
410 
411  inline Real inner(const vec2& a, const vec2& b)
412  {
413  vec2 tmp = a * b;
414  return sum(tmp);
415  }
416 
417  inline Real norm(const vec2& a)
418  {
419  return sqrt(inner(a, a));
420  }
421 
422  inline Real squaredNorm(const vec2& a) {
423  return inner(a, a);
424  }
425 
426  inline vec2 unit(const vec2& a)
427  {
428  Real n = norm(a);
429  if (n < epsilon)
430  return 0;
431  else
432  return a / n;
433  }
434 
435  inline Real angle(const vec2& a, const vec2& b)
436  {
437  Real ang = inner(unit(a), unit(b));
438  if (ang > 1 - epsilon)
439  return 0;
440  else if (ang < -1 + epsilon)
441  return M_PI;
442  else
443  return acos(ang);
444  }
445 
446  inline vec2 proj(const vec2& axis, const vec2& a)
447  {
448  vec2 u = unit(axis);
449  return inner(a, u) * u;
450  }
451 
452  inline vec2 absolute(const vec2& val)
453  {
454  return vec2(fabs(val[0]), fabs(val[1]));
455  }
456 
457  void store(FILE* fp, const vec2& v);
458  //int scan(FILE* fp, const vec2& v);
459 
460  int apx_equal(const vec2& a, const vec2& b);
461  int apx_equal(const vec2& a, const vec2& b, Real eps);
462 
466  struct __declspec(dllexport) vec3 {
467  Real v[3];
468  static const int n;
469 
470  inline vec3() { }
471  inline vec3(Real a) { v[1 - 1] = a; v[2 - 1] = a; v[3 - 1] = a; }
472  inline vec3(Real v_1, Real v_2, Real v_3) { v[1 - 1] = v_1; v[2 - 1] = v_2; v[3 - 1] = v_3; }
473  inline vec3(const ivec3& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; v[3 - 1] = a[3 - 1]; }
474  inline vec3(const vec3& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; v[3 - 1] = a[3 - 1]; }
475 
476  inline vec3& operator=(const vec3& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; v[3 - 1] = a[3 - 1]; return *this; }
477  inline Real& operator[] (int i) { return v[i]; }
478  inline const Real& operator[] (int i) const { return v[i]; }
479  inline Real& operator() (int i) { return v[i % 3]; }
480  inline const Real& operator() (int i) const { return v[i % 3]; }
481 
482  inline bool is_zero() const { return (v[0] == 0.0 && v[1] == 0.0 && v[2] == 0.0); }
483  inline bool is_tiny(Real tiny_tol = epsilon) const { return (fabs(v[0]) <= tiny_tol && fabs(v[1]) <= tiny_tol && fabs(v[2]) <= tiny_tol); }
484 
485  bool unit();
486  Real length() const;
487 
488 
489  //
490  // returns 1: this and other vectors are parallel
491  // -1: this and other vectors are anti-parallel
492  // 0: this and other vectors are not parallel
493  // or at least one of the vectors is zero
494  int is_parallel(
495  const vec3&, // other vector
496  Real = angle_tolerance // optional angle tolerance (radians)
497  ) const;
498 
499  //
500  // returns true: this and other vectors are perpendicular
501  // false: this and other vectors are not perpendicular
502  // or at least one of the vectors is zero
503  bool is_perpendicular(
504  const vec3&, // other vector
505  Real = angle_tolerance // optional angle tolerance (radians)
506  ) const;
507 
508  //
509  // set this vector to be perpendicular to another vector
510  bool perpendicular( // Result is not unitized.
511  // returns false if input vector is zero
512  const vec3&
513  );
514 
515  //
516  // set this vector to be perpendicular to a plane defined by 3 points
517  // returns false if points are coincident or colinear
518  bool perpendicular(
519  const vec3&, const vec3&, const vec3&
520  );
521  };
522 
523  //| binary op : componentwise
524 
525 
526  inline vec3 operator + (const vec3& a, const vec3& b)
527  {
528  vec3 c;
529  for (int i = 0; i < 3; ++i) { c[i] = a[i] + b[i]; }
530  return c;
531  }
532 
533  inline vec3 operator + (Real a, const vec3& b)
534  {
535  vec3 c;
536  for (int i = 0; i < 3; ++i) { c[i] = a + b[i]; }
537  return c;
538  }
539 
540  inline vec3 operator + (const vec3& a, Real b)
541  {
542  vec3 c;
543  for (int i = 0; i < 3; ++i) { c[i] = a[i] + b; }
544  return c;
545  }
546 
547 
548 
549  inline vec3 operator - (const vec3& a, const vec3& b)
550  {
551  vec3 c;
552  for (int i = 0; i < 3; ++i) { c[i] = a[i] - b[i]; }
553  return c;
554  }
555 
556  inline vec3 operator - (Real a, const vec3& b)
557  {
558  vec3 c;
559  for (int i = 0; i < 3; ++i) { c[i] = a - b[i]; }
560  return c;
561  }
562 
563  inline vec3 operator - (const vec3& a, Real b)
564  {
565  vec3 c;
566  for (int i = 0; i < 3; ++i) { c[i] = a[i] - b; }
567  return c;
568  }
569 
570 
571 
572  inline vec3 operator * (const vec3& a, const vec3& b)
573  {
574  vec3 c;
575  for (int i = 0; i < 3; ++i) { c[i] = a[i] * b[i]; }
576  return c;
577  }
578 
579  inline vec3 operator * (Real a, const vec3& b)
580  {
581  vec3 c;
582  for (int i = 0; i < 3; ++i) { c[i] = a * b[i]; }
583  return c;
584  }
585 
586  inline vec3 operator * (const vec3& a, Real b)
587  {
588  vec3 c;
589  for (int i = 0; i < 3; ++i) { c[i] = a[i] * b; }
590  return c;
591  }
592 
593 
594 
595  inline vec3 operator / (const vec3& a, const vec3& b)
596  {
597  vec3 c;
598  for (int i = 0; i < 3; ++i) { c[i] = a[i] / b[i]; }
599  return c;
600  }
601 
602  inline vec3 operator / (Real a, const vec3& b)
603  {
604  vec3 c;
605  for (int i = 0; i < 3; ++i) { c[i] = a / b[i]; }
606  return c;
607  }
608 
609  inline vec3 operator / (const vec3& a, Real b)
610  {
611  vec3 c;
612  for (int i = 0; i < 3; ++i) { c[i] = a[i] / b; }
613  return c;
614  }
615 
616  //| cumulative op : componentwise
617 
618 
619  inline vec3 operator += (vec3& a, const vec3& b)
620  {
621  return a = (a + b);
622  }
623 
624  inline vec3 operator += (vec3& a, Real b)
625  {
626  return a = (a + b);
627  }
628 
629 
630 
631  inline vec3 operator -= (vec3& a, const vec3& b)
632  {
633  return a = (a - b);
634  }
635 
636  inline vec3 operator -= (vec3& a, Real b)
637  {
638  return a = (a - b);
639  }
640 
641 
642 
643  inline vec3 operator *= (vec3& a, const vec3& b)
644  {
645  return a = (a * b);
646  }
647 
648  inline vec3 operator *= (vec3& a, Real b)
649  {
650  return a = (a * b);
651  }
652 
653 
654 
655  inline vec3 operator /= (vec3& a, const vec3& b)
656  {
657  return a = (a / b);
658  }
659 
660  inline vec3 operator /= (vec3& a, Real b)
661  {
662  return a = (a / b);
663  }
664 
665 
666 
667  //| logical op : componentwise
668 
669 
670  inline int operator == (const vec3& a, const vec3& b)
671  {
672  int c = 1;
673  for (int i = 0; i < 3; ++i) { c = c && a[i] == b[i]; }
674  return c;
675  }
676 
677  inline int operator == (Real a, const vec3& b)
678  {
679  int c = 1;
680  for (int i = 0; i < 3; ++i) { c = c && a == b[i]; }
681  return c;
682  }
683 
684  inline int operator == (const vec3& a, Real b)
685  {
686  int c = 1;
687  for (int i = 0; i < 3; ++i) { c = c && a[i] == b; }
688  return c;
689  }
690 
691 
692 
693  inline int operator < (const vec3& a, const vec3& b)
694  {
695  int c = 1;
696  for (int i = 0; i < 3; ++i) { c = c && a[i] < b[i]; }
697  return c;
698  }
699 
700  inline int operator < (Real a, const vec3& b)
701  {
702  int c = 1;
703  for (int i = 0; i < 3; ++i) { c = c && a < b[i]; }
704  return c;
705  }
706 
707  inline int operator < (const vec3& a, Real b)
708  {
709  int c = 1;
710  for (int i = 0; i < 3; ++i) { c = c && a[i] < b; }
711  return c;
712  }
713 
714 
715 
716  inline int operator <= (const vec3& a, const vec3& b)
717  {
718  int c = 1;
719  for (int i = 0; i < 3; ++i) { c = c && a[i] <= b[i]; }
720  return c;
721  }
722 
723  inline int operator <= (Real a, const vec3& b)
724  {
725  int c = 1;
726  for (int i = 0; i < 3; ++i) { c = c && a <= b[i]; }
727  return c;
728  }
729 
730  inline int operator <= (const vec3& a, Real b)
731  {
732  int c = 1;
733  for (int i = 0; i < 3; ++i) { c = c && a[i] <= b; }
734  return c;
735  }
736 
737 
738 
739  inline int operator > (const vec3& a, const vec3& b)
740  {
741  int c = 1;
742  for (int i = 0; i < 3; ++i) { c = c && a[i] > b[i]; }
743  return c;
744  }
745 
746  inline int operator > (Real a, const vec3& b)
747  {
748  int c = 1;
749  for (int i = 0; i < 3; ++i) { c = c && a > b[i]; }
750  return c;
751  }
752 
753  inline int operator > (const vec3& a, Real b)
754  {
755  int c = 1;
756  for (int i = 0; i < 3; ++i) { c = c && a[i] > b; }
757  return c;
758  }
759 
760 
761 
762  inline int operator >= (const vec3& a, const vec3& b)
763  {
764  int c = 1;
765  for (int i = 0; i < 3; ++i) { c = c && a[i] >= b[i]; }
766  return c;
767  }
768 
769  inline int operator >= (Real a, const vec3& b)
770  {
771  int c = 1;
772  for (int i = 0; i < 3; ++i) { c = c && a >= b[i]; }
773  return c;
774  }
775 
776  inline int operator >= (const vec3& a, Real b)
777  {
778  int c = 1;
779  for (int i = 0; i < 3; ++i) { c = c && a[i] >= b; }
780  return c;
781  }
782 
783 
784 
785  //| unary op : componentwise
786  inline vec3 operator - (const vec3& a)
787  {
788  vec3 c;
789  for (int i = 0; i < 3; ++i) { c[i] = -a[i]; }
790  return c;
791  }
792 
793  inline vec3 absolute(const vec3& val)
794  {
795  return vec3(fabs(val[0]), fabs(val[1]), fabs(val[2]));
796  }
797 
798 
799 
800  //| R^n -> R
801  inline Real sum(const vec3& a)
802  {
803  Real s = 0;
804 
805  s += a[1 - 1];
806 
807  s += a[2 - 1];
808 
809  s += a[3 - 1];
810 
811  return s;
812  }
813 
814  inline Real inner(const vec3& a, const vec3& b)
815  {
816  vec3 tmp = a * b;
817  return sum(tmp);
818  }
819 
820  inline Real squaredNorm(const vec3& a) {
821  return inner(a, a);
822  }
823 
824  inline Real norm(const vec3& a)
825  {
826  return sqrt(inner(a, a));
827  }
828 
829  inline vec3 unit(const vec3& a)
830  {
831  Real n = norm(a);
832  if (n < zero_epsilon)
833  return 0;
834  else
835  return a / n;
836  }
837 
838  inline Real angle(const vec3& a, const vec3& b)
839  {
840  Real ang = inner(unit(a), unit(b));
841  if (ang > 1 - epsilon)
842  return 0;
843  else if (ang < -1 + epsilon)
844  return M_PI;
845  else
846  return acos(ang);
847  }
848 
849  inline vec3 proj(const vec3& axis, const vec3& a)
850  {
851  vec3 u = unit(axis);
852  return inner(a, u) * u;
853  }
854 
855  void store(FILE* fp, const vec3& v);
856  //int scan(FILE* fp, const vec3& v);
857 
858  int apx_equal(const vec3& a, const vec3& b);
859  int apx_equal(const vec3& a, const vec3& b, Real eps);
860 
864  struct __declspec(dllexport) vec4 {
865  Real v[4];
866  static const int n;
867 
868  inline vec4() { }
869  inline vec4(Real a) { v[1 - 1] = a; v[2 - 1] = a; v[3 - 1] = a; v[4 - 1] = a; }
870  inline vec4(Real v_1, Real v_2, Real v_3, Real v_4) { v[1 - 1] = v_1; v[2 - 1] = v_2; v[3 - 1] = v_3; v[4 - 1] = v_4; }
871  inline vec4(const ivec4& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; v[3 - 1] = a[3 - 1]; v[4 - 1] = a[4 - 1]; }
872  inline vec4(const vec4& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; v[3 - 1] = a[3 - 1]; v[4 - 1] = a[4 - 1]; }
873 
874  inline vec4& operator=(const vec4& a) { v[1 - 1] = a[1 - 1]; v[2 - 1] = a[2 - 1]; v[3 - 1] = a[3 - 1]; v[4 - 1] = a[4 - 1]; return *this; }
875  inline Real& operator[] (int i) { return v[i]; }
876  inline const Real& operator[] (int i) const { return v[i]; }
877  inline Real& operator() (int i) { return v[i % 4]; }
878  inline const Real& operator() (int i) const { return v[i % 4]; }
879 
880  inline bool is_zero() const { return (v[0] == 0.0 && v[1] == 0.0 && v[2] == 0.0 && v[3] == 0.0); }
881  inline bool is_tiny(Real tiny_tol = epsilon) const {
882  return (fabs(v[0]) <= tiny_tol && fabs(v[1]) <= tiny_tol && fabs(v[2]) <= tiny_tol && fabs(v[3]) <= tiny_tol);
883  }
884 
885  bool unit();
886  Real length() const;
887  };
888 
889 
890 
891 
892 
893  //| binary op : componentwise
894 
895 
896  inline vec4 operator + (const vec4& a, const vec4& b)
897  {
898  vec4 c;
899  for (int i = 0; i < 4; ++i) { c[i] = a[i] + b[i]; }
900  return c;
901  }
902 
903  inline vec4 operator + (Real a, const vec4& b)
904  {
905  vec4 c;
906  for (int i = 0; i < 4; ++i) { c[i] = a + b[i]; }
907  return c;
908  }
909 
910  inline vec4 operator + (const vec4& a, Real b)
911  {
912  vec4 c;
913  for (int i = 0; i < 4; ++i) { c[i] = a[i] + b; }
914  return c;
915  }
916 
917 
918 
919  inline vec4 operator - (const vec4& a, const vec4& b)
920  {
921  vec4 c;
922  for (int i = 0; i < 4; ++i) { c[i] = a[i] - b[i]; }
923  return c;
924  }
925 
926  inline vec4 operator - (Real a, const vec4& b)
927  {
928  vec4 c;
929  for (int i = 0; i < 4; ++i) { c[i] = a - b[i]; }
930  return c;
931  }
932 
933  inline vec4 operator - (const vec4& a, Real b)
934  {
935  vec4 c;
936  for (int i = 0; i < 4; ++i) { c[i] = a[i] - b; }
937  return c;
938  }
939 
940 
941 
942  inline vec4 operator * (const vec4& a, const vec4& b)
943  {
944  vec4 c;
945  for (int i = 0; i < 4; ++i) { c[i] = a[i] * b[i]; }
946  return c;
947  }
948 
949  inline vec4 operator * (Real a, const vec4& b)
950  {
951  vec4 c;
952  for (int i = 0; i < 4; ++i) { c[i] = a * b[i]; }
953  return c;
954  }
955 
956  inline vec4 operator * (const vec4& a, Real b)
957  {
958  vec4 c;
959  for (int i = 0; i < 4; ++i) { c[i] = a[i] * b; }
960  return c;
961  }
962 
963 
964 
965  inline vec4 operator / (const vec4& a, const vec4& b)
966  {
967  vec4 c;
968  for (int i = 0; i < 4; ++i) { c[i] = a[i] / b[i]; }
969  return c;
970  }
971 
972  inline vec4 operator / (Real a, const vec4& b)
973  {
974  vec4 c;
975  for (int i = 0; i < 4; ++i) { c[i] = a / b[i]; }
976  return c;
977  }
978 
979  inline vec4 operator / (const vec4& a, Real b)
980  {
981  vec4 c;
982  for (int i = 0; i < 4; ++i) { c[i] = a[i] / b; }
983  return c;
984  }
985 
986 
987 
988  //| cumulative op : componentwise
989 
990 
991  inline vec4 operator += (vec4& a, const vec4& b)
992  {
993  return a = (a + b);
994  }
995 
996  inline vec4 operator += (vec4& a, Real b)
997  {
998  return a = (a + b);
999  }
1000 
1001 
1002 
1003  inline vec4 operator -= (vec4& a, const vec4& b)
1004  {
1005  return a = (a - b);
1006  }
1007 
1008  inline vec4 operator -= (vec4& a, Real b)
1009  {
1010  return a = (a - b);
1011  }
1012 
1013 
1014 
1015  inline vec4 operator *= (vec4& a, const vec4& b)
1016  {
1017  return a = (a * b);
1018  }
1019 
1020  inline vec4 operator *= (vec4& a, Real b)
1021  {
1022  return a = (a * b);
1023  }
1024 
1025 
1026 
1027  inline vec4 operator /= (vec4& a, const vec4& b)
1028  {
1029  return a = (a / b);
1030  }
1031 
1032  inline vec4 operator /= (vec4& a, Real b)
1033  {
1034  return a = (a / b);
1035  }
1036 
1037 
1038 
1039  //| logical op : componentwise
1040 
1041 
1042  inline int operator == (const vec4& a, const vec4& b)
1043  {
1044  int c = 1;
1045  for (int i = 0; i < 4; ++i) { c = c && a[i] == b[i]; }
1046  return c;
1047  }
1048 
1049  inline int operator == (Real a, const vec4& b)
1050  {
1051  int c = 1;
1052  for (int i = 0; i < 4; ++i) { c = c && a == b[i]; }
1053  return c;
1054  }
1055 
1056  inline int operator == (const vec4& a, Real b)
1057  {
1058  int c = 1;
1059  for (int i = 0; i < 4; ++i) { c = c && a[i] == b; }
1060  return c;
1061  }
1062 
1063 
1064 
1065  inline int operator < (const vec4& a, const vec4& b)
1066  {
1067  int c = 1;
1068  for (int i = 0; i < 4; ++i) { c = c && a[i] < b[i]; }
1069  return c;
1070  }
1071 
1072  inline int operator < (Real a, const vec4& b)
1073  {
1074  int c = 1;
1075  for (int i = 0; i < 4; ++i) { c = c && a < b[i]; }
1076  return c;
1077  }
1078 
1079  inline int operator < (const vec4& a, Real b)
1080  {
1081  int c = 1;
1082  for (int i = 0; i < 4; ++i) { c = c && a[i] < b; }
1083  return c;
1084  }
1085 
1086 
1087 
1088  inline int operator <= (const vec4& a, const vec4& b)
1089  {
1090  int c = 1;
1091  for (int i = 0; i < 4; ++i) { c = c && a[i] <= b[i]; }
1092  return c;
1093  }
1094 
1095  inline int operator <= (Real a, const vec4& b)
1096  {
1097  int c = 1;
1098  for (int i = 0; i < 4; ++i) { c = c && a <= b[i]; }
1099  return c;
1100  }
1101 
1102  inline int operator <= (const vec4& a, Real b)
1103  {
1104  int c = 1;
1105  for (int i = 0; i < 4; ++i) { c = c && a[i] <= b; }
1106  return c;
1107  }
1108 
1109 
1110 
1111  inline int operator > (const vec4& a, const vec4& b)
1112  {
1113  int c = 1;
1114  for (int i = 0; i < 4; ++i) { c = c && a[i] > b[i]; }
1115  return c;
1116  }
1117 
1118  inline int operator > (Real a, const vec4& b)
1119  {
1120  int c = 1;
1121  for (int i = 0; i < 4; ++i) { c = c && a > b[i]; }
1122  return c;
1123  }
1124 
1125  inline int operator > (const vec4& a, Real b)
1126  {
1127  int c = 1;
1128  for (int i = 0; i < 4; ++i) { c = c && a[i] > b; }
1129  return c;
1130  }
1131 
1132 
1133 
1134  inline int operator >= (const vec4& a, const vec4& b)
1135  {
1136  int c = 1;
1137  for (int i = 0; i < 4; ++i) { c = c && a[i] >= b[i]; }
1138  return c;
1139  }
1140 
1141  inline int operator >= (Real a, const vec4& b)
1142  {
1143  int c = 1;
1144  for (int i = 0; i < 4; ++i) { c = c && a >= b[i]; }
1145  return c;
1146  }
1147 
1148  inline int operator >= (const vec4& a, Real b)
1149  {
1150  int c = 1;
1151  for (int i = 0; i < 4; ++i) { c = c && a[i] >= b; }
1152  return c;
1153  }
1154 
1155 
1156 
1157  //| unary op : componentwise
1158  inline vec4 operator - (const vec4& a)
1159  {
1160  vec4 c;
1161  for (int i = 0; i < 4; ++i) { c[i] = -a[i]; }
1162  return c;
1163  }
1164 
1165  inline vec4 absolute(const vec4& val)
1166  {
1167  return vec4(fabs(val[0]), fabs(val[1]), fabs(val[2]), fabs(val[3]));
1168  }
1169 
1170 
1171  //| R^n -> R
1172  inline Real sum(const vec4& a)
1173  {
1174  Real s = 0;
1175 
1176  s += a[1 - 1];
1177 
1178  s += a[2 - 1];
1179 
1180  s += a[3 - 1];
1181 
1182  s += a[4 - 1];
1183 
1184  return s;
1185  }
1186 
1187  inline Real inner(const vec4& a, const vec4& b)
1188  {
1189  vec4 tmp = a * b;
1190  return sum(tmp);
1191  }
1192  inline Real squaredNorm(const vec4& a) {
1193  return inner(a, a);
1194  }
1195  inline Real norm(const vec4& a)
1196  {
1197  return sqrt(inner(a, a));
1198  }
1199 
1200  inline vec4 unit(const vec4& a)
1201  {
1202  Real n = norm(a);
1203  if (n < epsilon)
1204  return 0;
1205  else
1206  return a / n;
1207  }
1208 
1209  inline Real angle(const vec4& a, const vec4& b)
1210  {
1211  Real ang = inner(unit(a), unit(b));
1212  if (ang > 1 - epsilon)
1213  return 0;
1214  else if (ang < -1 + epsilon)
1215  return M_PI;
1216  else
1217  return acos(ang);
1218  }
1219 
1220  inline vec4 proj(const vec4& axis, const vec4& a)
1221  {
1222  vec4 u = unit(axis);
1223  return inner(a, u) * u;
1224  }
1225 
1226  void store(FILE* fp, const vec4& v);
1227 
1228  //int scan(FILE* fp, const vec4& v);
1229 
1230  int apx_equal(const vec4& a, const vec4& b);
1231  int apx_equal(const vec4& a, const vec4& b, Real eps);
1232 
1233  vec3 cross(const vec3& a, const vec3& b);
1234 
1235 
1236 }; //namespace oph
1237 
1238 #endif // !__vec_h
Real inner(const vec2 &a, const vec2 &b)
Definition: vec.h:411
ivec2 operator+=(ivec2 &a, const ivec2 &b)
Definition: ivec.h:194
vec2 proj(const vec2 &axis, const vec2 &a)
Definition: vec.h:446
Real norm(const vec2 &a)
Definition: vec.h:417
ivec2 operator-(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:122
vec2 absolute(const vec2 &val)
Definition: vec.h:452
int operator>=(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:338
Real squaredNorm(const vec2 &a)
Definition: vec.h:422
ivec2 operator+(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:99
void store(FILE *fp, const vec2 &v)
Definition: vec.cpp:127
Real epsilon
Definition: epsilon.cpp:53
ivec2 operator*(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:145
int apx_equal(Real x, Real y)
Definition: epsilon.cpp:91
vec2 operator/=(vec2 &a, const vec2 &b)
Definition: vec.h:262
int operator<=(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:291
vec3 cross(const vec3 &a, const vec3 &b)
Definition: vec.cpp:444
vec2 operator/(const vec2 &a, const vec2 &b)
Definition: vec.h:200
ivec2 operator-=(ivec2 &a, const ivec2 &b)
Definition: ivec.h:206
Real zero_epsilon
Definition: epsilon.cpp:59
#define M_PI
Definition: define.h:52
int operator>(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:314
int operator==(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:245
void angle(const std::vector< Complex< T >> &src, std::vector< T > &dst)
Definition: function.h:159
float Real
Definition: typedef.h:55
int operator<(const ivec2 &a, const ivec2 &b)
Definition: ivec.h:268
Real angle_tolerance
Definition: epsilon.cpp:60
Definition: Bitmap.h:49
Real sum(const vec2 &a)
Definition: vec.h:401
ivec2 operator*=(ivec2 &a, const ivec2 &b)
Definition: ivec.h:218
vec2 unit(const vec2 &a)
Definition: vec.h:426