Openholo  v1.0
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
float Real
Definition: typedef.h:55
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
int i
Definition: ophTriMesh.cpp:413
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
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
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
#define M_PI
Definition: define.h:52
vec2 unit(const vec2 &a)
Definition: vec.h:426