PCL
Complex.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 02.01.12.0947
6 // ----------------------------------------------------------------------------
7 // pcl/Complex.h - Released 2019-04-30T16:30:41Z
8 // ----------------------------------------------------------------------------
9 // This file is part of the PixInsight Class Library (PCL).
10 // PCL is a multiplatform C++ framework for development of PixInsight modules.
11 //
12 // Copyright (c) 2003-2019 Pleiades Astrophoto S.L. All Rights Reserved.
13 //
14 // Redistribution and use in both source and binary forms, with or without
15 // modification, is permitted provided that the following conditions are met:
16 //
17 // 1. All redistributions of source code must retain the above copyright
18 // notice, this list of conditions and the following disclaimer.
19 //
20 // 2. All redistributions in binary form must reproduce the above copyright
21 // notice, this list of conditions and the following disclaimer in the
22 // documentation and/or other materials provided with the distribution.
23 //
24 // 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
25 // of their contributors, may be used to endorse or promote products derived
26 // from this software without specific prior written permission. For written
27 // permission, please contact info@pixinsight.com.
28 //
29 // 4. All products derived from this software, in any form whatsoever, must
30 // reproduce the following acknowledgment in the end-user documentation
31 // and/or other materials provided with the product:
32 //
33 // "This product is based on software from the PixInsight project, developed
34 // by Pleiades Astrophoto and its contributors (http://pixinsight.com/)."
35 //
36 // Alternatively, if that is where third-party acknowledgments normally
37 // appear, this acknowledgment must be reproduced in the product itself.
38 //
39 // THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
40 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
41 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
42 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
43 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44 // EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
45 // INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
46 // DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
47 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 // POSSIBILITY OF SUCH DAMAGE.
50 // ----------------------------------------------------------------------------
51 
52 #ifndef __PCL_Complex_h
53 #define __PCL_Complex_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
60 #include <pcl/Constants.h>
61 #include <pcl/Math.h>
62 #include <pcl/Relational.h>
63 
64 namespace pcl
65 {
66 
67 // ----------------------------------------------------------------------------
68 
69 #define real C[0]
70 #define imag C[1]
71 
82 template <typename T>
83 class PCL_CLASS Complex
84 {
85 public:
86 
90  typedef T component;
91 
97  {
98  }
99 
103  Complex( T r, T i = 0 )
104  {
105  real = r;
106  imag = i;
107  }
108 
112  template <typename T1>
113  Complex( const Complex<T1>& c )
114  {
115  real = T( c.Real() );
116  imag = T( c.Imag() );
117  }
118 
122  constexpr T Real() const
123  {
124  return real;
125  }
126 
130  T& Real()
131  {
132  return real;
133  }
134 
138  constexpr T Imag() const
139  {
140  return imag;
141  }
142 
146  T& Imag()
147  {
148  return imag;
149  }
150 
154  constexpr bool IsReal() const
155  {
156  return imag == 0;
157  }
158 
162  template <typename T1>
163  Complex<T>& operator =( const Complex<T1>& c )
164  {
165  real = T( c.Real() ), imag = T( c.Imag() );
166  return *this;
167  }
168 
173  template <typename T1>
174  Complex<T>& operator +=( const Complex<T1>& c )
175  {
176  real += c.Real(), imag += c.Imag();
177  return *this;
178  }
179 
184  template <typename T1>
185  Complex<T>& operator -=( const Complex<T1>& c )
186  {
187  real -= c.Real(), imag -= c.Imag();
188  return *this;
189  }
190 
195  template <typename T1>
196  Complex<T>& operator *=( const Complex<T1>& c )
197  {
198  T t = T( real*c.Real() - imag*c.Imag() );
199  imag = T( imag*c.Real() + real*c.Imag() );
200  real = t;
201  return *this;
202  }
203 
208  template <typename T1>
209  Complex<T>& operator /=( const Complex<T1>& c )
210  {
211  T r, d, t;
212  if ( pcl::Abs( c.Real() ) >= pcl::Abs( c.Imag() ) )
213  {
214  PCL_PRECONDITION( c.Real() != 0 )
215  r = T( c.Imag()/c.Real() );
216  d = T( c.Real() + r*c.Imag() );
217  t = T( (real + r*imag)/d );
218  imag = (imag - r*real)/d;
219  }
220  else
221  {
222  PCL_PRECONDITION( c.Imag() != 0 )
223  r = T( c.Real()/c.Imag() );
224  d = T( c.Imag() + r*c.Real() );
225  t = T( (real*r + imag)/d );
226  imag = (imag*r - real)/d;
227  }
228  real = t;
229  return *this;
230  }
231 
237  template <typename T1>
238  Complex<T>& operator =( T1 x )
239  {
240  real = x, imag = 0;
241  return *this;
242  }
243 
248  template <typename T1>
249  Complex<T>& operator +=( T1 x )
250  {
251  real += x;
252  return *this;
253  }
254 
260  template <typename T1>
261  Complex<T>& operator -=( T1 x )
262  {
263  real -= x;
264  return *this;
265  }
266 
272  template <typename T1>
273  Complex<T>& operator *=( T1 x )
274  {
275  real *= x, imag *= x;
276  return *this;
277  }
278 
284  template <typename T1>
285  Complex<T>& operator /=( T1 x )
286  {
287  PCL_PRECONDITION( x != 0 )
288  real /= x, imag /= x;
289  return *this;
290  }
291 
296  {
297  return *this;
298  }
299 
307  {
308  return Complex<T>( -real, -imag );
309  }
310 
315  Complex<T> Conj() const
316  {
317  return Complex<T>( real, -imag );
318  }
319 
323  Complex<T> operator ~() const
324  {
325  return Conj();
326  }
327 
332  void SetConj()
333  {
334  imag = -imag;
335  }
336 
343  T Mag() const
344  {
345  T r = pcl::Abs( real );
346  T i = pcl::Abs( imag );
347  T m;
348  if ( r == 0 )
349  m = i;
350  else if ( i == 0 )
351  m = r;
352  else
353  {
354  bool q = r < i;
355  m = q ? r/i : i/r;
356  m = (q ? i : r) * pcl::Sqrt( 1 + m*m );
357  }
358  return m;
359  }
360 
365  constexpr T Norm() const
366  {
367  return real*real + imag*imag;
368  }
369 
375  constexpr T Arg() const
376  {
377  // Degenerate cases (real=0) are correctly handled by real ArcTan(). For
378  // the undefined case real=imag=0, we silently return zero. Should we
379  // throw an exception instead?
380  return (real != 0 || imag != 0) ? pcl::ArcTan( imag, real ) : 0;
381  }
382 
383 private:
384 
385  /*
386  * C99 binary-compatible complex components.
387  */
388  T C[ 2 ]; // C[0] = real, C[1] = imaginary
389 };
390 
391 #undef real
392 #undef imag
393 
394 /*
395  * ### N.B.: Template class Complex<T> cannot have virtual member functions.
396  * This is because sizeof( Complex<T> ) _must_ be equal to 2*sizeof( T ).
397  */
398 template <typename T>
399 struct PCL_AssertComplexSize
400 {
401  static_assert( sizeof( Complex<T> ) == 2*sizeof( T ), "Invalid sizeof( Complex<> )" );
402 };
403 
404 // ----------------------------------------------------------------------------
405 
419 template <typename T> inline
420 T Abs( const Complex<T>& c )
421 {
422  return c.Mag();
423 }
424 
431 template <typename T> inline
432 Complex<T> Polar( T r, T stheta, T ctheta )
433 {
434  return Complex<T>( r*ctheta, r*stheta );
435 }
436 
442 template <typename T> inline
443 Complex<T> Polar( T r, T theta )
444 {
445  return Polar( r, pcl::Sin( theta ), pcl::Cos( theta ) );
446 }
447 
448 // ----------------------------------------------------------------------------
449 
454 template <typename T1, class T2> inline
456 {
457  return Complex<T1>( T1( c1.Real() + c2.Real() ),
458  T1( c1.Imag() + c2.Imag() ) );
459 }
460 
469 template <typename T1, class T2> inline
471 {
472  return Complex<T1>( T1( c.Real()+x ), c.Imag() );
473 }
474 
483 template <typename T1, class T2> inline
485 {
486  return c + x;
487 }
488 
494 template <typename T1, class T2> inline
496 {
497  return Complex<T1>( T1( c1.Real() - c2.Real() ),
498  T1( c1.Imag() - c2.Imag() ) );
499 }
500 
511 template <typename T1, class T2> inline
513 {
514  return Complex<T1>( T1( c.Real()-x ), c.Imag() );
515 }
516 
527 template <typename T1, class T2> inline
529 {
530  return Complex<T2>( T2( x-c.Real() ), -c.Imag() );
531 }
532 
538 template <typename T1, class T2> inline
540 {
541  return Complex<T1>( T1( c1.Real()*c2.Real() - c1.Imag()*c2.Imag() ),
542  T1( c1.Imag()*c2.Real() + c1.Real()*c2.Imag() ) );
543 }
544 
555 template <typename T1, class T2> inline
557 {
558  return Complex<T1>( T1( c.Real()*x ), c.Imag()*x );
559 }
560 
571 template <typename T1, class T2> inline
573 {
574  return c * x;
575 }
576 
582 template <typename T1, class T2> inline
584 {
585  Complex<T1> c;
586  T2 r, d;
587  if ( pcl::Abs( c2.Real() ) >= pcl::Abs( c2.Imag() ) )
588  {
589  PCL_PRECONDITION( c2.Real() != 0 )
590  r = c2.Imag() / c2.Real();
591  d = c2.Real() + r*c2.Imag();
592  c.Real() = T1( (c1.Real() + r*c1.Imag())/d );
593  c.Imag() = T1( (c1.Imag() - r*c1.Real())/d );
594  }
595  else
596  {
597  PCL_PRECONDITION( c2.Imag() != 0 )
598  r = c2.Real() / c2.Imag();
599  d = c2.Imag() + r*c2.Real();
600  c.Real() = T1( (c1.Real()*r + c1.Imag())/d );
601  c.Imag() = T1( (c1.Imag()*r - c1.Real())/d );
602  }
603  return c;
604 }
605 
616 template <typename T1, class T2> inline
618 {
619  PCL_PRECONDITION( x != 0 )
620  return Complex<T1>( T1( c.Real()/x ), T1( c.Imag()/x ) );
621 }
622 
633 template <typename T1, class T2> inline
635 {
636  // return Complex( x, 0 )/c;
637  Complex<T2> c3;
638  T2 r, d;
639  if ( pcl::Abs( c.Real() ) >= pcl::Abs( c.Imag() ) )
640  {
641  PCL_PRECONDITION( c.Real() != 0 )
642  r = c.Imag()/c.Real();
643  d = c.Real() + r*c.Imag();
644  c3.Real() = T2( x/d );
645  c3.Imag() = T2( -((r*x)/d) );
646  }
647  else
648  {
649  PCL_PRECONDITION( c.Imag() != 0 )
650  r = c.Real()/c.Imag();
651  d = c.Imag() + r*c.Real();
652  c3.Real() = T2( (r*x)/d );
653  c3.Imag() = T2( -(x/d) );
654  }
655  return c3;
656 }
657 
658 // ----------------------------------------------------------------------------
659 
664 template <typename T> inline
666 {
667  if ( c.Real() == 0 && c.Imag() == 0 )
668  return Complex<T>( 0 );
669 
670  Complex<T> c1;
671  T m, r = pcl::Abs( c.Real() ), i = pcl::Abs( c.Imag() );
672 
673  if ( r >= i )
674  {
675  PCL_PRECONDITION( r != 0 )
676  T t = i/r;
677  m = pcl::Sqrt( r ) * pcl::Sqrt( (1 + pcl::Sqrt( 1 + t*t ))/2 );
678  }
679  else
680  {
681  PCL_PRECONDITION( i != 0 )
682  T t = r/i;
683  m = pcl::Sqrt( i ) * pcl::Sqrt( (t + pcl::Sqrt( 1 + t*t ))/2 );
684  }
685 
686  if ( c.Real() >= 0 )
687  {
688  c1.Real() = m;
689  c1.Imag() = c.Imag()/(m+m);
690  }
691  else
692  {
693  c1.Imag() = (c.Imag() >= 0) ? m : -m;
694  c1.Real() = c.Imag()/(c1.Imag()+c1.Imag());
695  }
696 
697  return c1;
698 }
699 
704 template <typename T> inline
706 {
707  T x = pcl::Exp( c.Real() );
708  return Complex<T>( x*pcl::Cos( c.Imag() ), x*pcl::Sin( c.Imag() ) );
709 }
710 
715 template <typename T> inline
717 {
718  return Complex<T>( pcl::Ln( c.Mag() ), c.Arg() );
719 }
720 
725 template <typename T> inline
727 {
728  return pcl::Const<T>::log10e() * pcl::Ln( c );
729 }
730 
731 // ----------------------------------------------------------------------------
732 
737 template <typename T1, class T2> inline
738 Complex<T1> Pow( const Complex<T1>& c, T2 x )
739 {
740  if ( c.Imag() == 0 )
741  return Complex<T1>( pcl::Pow( c.Real(), T1( x ) ) );
742  else
743  return pcl::Exp( T1( x )*pcl::Ln( c ) );
744 }
745 
750 template <typename T1, class T2> inline
751 Complex<T2> Pow( T1 x, const Complex<T2>& c )
752 {
753  if ( c.Imag() == 0 )
754  return Complex<T2>( pcl::Pow( T2( x ), c.Real() ) );
755  else
756  return pcl::Exp( c*pcl::Ln( T2( x ) ) );
757 }
758 
764 template <typename T> inline
765 Complex<T> Pow( const Complex<T>& c1, const Complex<T>& c2 )
766 {
767  if ( c2.Imag() == 0 )
768  return pcl::Pow( c1, c2.Real() );
769  else if ( c1.Imag() == 0 )
770  return Complex<T>( pcl::Pow( c1.Real(), c2 ) );
771  else
772  return pcl::Exp( c2*pcl::Ln( c1 ) );
773 }
774 
775 // ----------------------------------------------------------------------------
776 
785 template <typename T> inline
787 {
788  return Complex<T>( pcl::Sin( c.Real() )*pcl::Cosh( c.Imag() ),
789  pcl::Cos( c.Real() )*pcl::Sinh( c.Imag() ) );
790 }
791 
796 template <typename T> inline
798 {
799  return Complex<T>( pcl::Cos( c.Real() )*pcl::Cosh( c.Imag() ),
800  -pcl::Sin( c.Real() )*pcl::Sinh( c.Imag() ) );
801 }
802 
807 template <typename T> inline
809 {
810  return pcl::Sin( c )/pcl::Cos( c );
811 }
812 
817 template <typename T> inline
819 {
820  return Complex<T>( pcl::Sinh( c.Real() )*pcl::Cos( c.Imag() ),
821  pcl::Cosh( c.Real() )*pcl::Sin( c.Imag() ) );
822 }
823 
828 template <typename T> inline
830 {
831  return Complex<T>( pcl::Cosh( c.Real() )*pcl::Cos( c.Imag() ),
832  pcl::Sinh( c.Real() )*pcl::Sin( c.Imag() ) );
833 }
834 
839 template <typename T> inline
841 {
842  return pcl::Sinh( c )/pcl::Cosh( c );
843 }
844 
845 // ----------------------------------------------------------------------------
846 
861 template <typename T1, class T2> inline
862 bool operator ==( const Complex<T1>& c1, const Complex<T2>& c2 )
863 {
864  return c1.Real() == c2.Real() && c1.Imag() == c2.Imag();
865 }
866 
871 template <typename T1, class T2> inline
872 bool operator ==( const Complex<T1>& c, T2 x )
873 {
874  return c.Real() == x && c.Imag() == T1( 0 );
875 }
876 
881 template <typename T1, class T2> inline
882 bool operator ==( T1 x, const Complex<T2>& c )
883 {
884  return c == x;
885 }
886 
891 template <typename T1, class T2> inline
892 bool operator <( const Complex<T1>& c1, const Complex<T2>& c2 )
893 {
894  return c1.Mag() < c2.Mag();
895 }
896 
901 template <typename T1, class T2> inline
902 bool operator <( const Complex<T1>& c, T2 x )
903 {
904  return c.Mag() < pcl::Abs( x );
905 }
906 
911 template <typename T1, class T2> inline
912 bool operator <( T1 x, const Complex<T2>& c )
913 {
914  return pcl::Abs( x ) < c.Mag();
915 }
916 
917 // ----------------------------------------------------------------------------
918 
928 template <typename T> inline
930 {
931  return Complex<T>( pcl::Round( c.Real() ), pcl::Round( c.Imag() ) );
932 }
933 
940 template <typename T> inline
941 Complex<T> Round( const Complex<T>& c, int n )
942 {
943  PCL_PRECONDITION( n >= 0 )
944  return Complex<T>( pcl::Round( c.Real(), n ), pcl::Round( c.Imag(), n ) );
945 }
946 
947 // ----------------------------------------------------------------------------
948 
973 template <typename T> inline
974 void PhaseCorrelationMatrix( Complex<T>* i, const Complex<T>* j, const Complex<T>* a, const Complex<T>* b )
975 {
976  const T tiny = T( 1.0e-20 );
977  for ( ; i < j; ++i, ++a, ++b )
978  {
979  Complex<T> n = *a * ~*b;
980  *i = n/Max( tiny, Abs( n ) );
981  }
982 }
983 
1004 template <typename T> inline
1005 void CrossPowerSpectrumMatrix( Complex<T>* i, const Complex<T>* j, const Complex<T>* a, const Complex<T>* b )
1006 {
1007  const T tiny = T( 1.0e-20 );
1008  for ( ; i < j; ++i, ++a, ++b )
1009  *i = (*b * ~*a)/Max( tiny, Abs( *a ) * Abs( *b ) );
1010 }
1011 
1012 //
1013 // Stream extraction/insertion.
1014 //
1015 /* ### PCL 2.x: Add these along with pcl::DataStream, etc...
1016 template <class S, typename T> inline
1017 S& operator >>( S& s, Complex<T>& c )
1018 {
1019  return s >> c.Real() >> c.Imag();
1020 }
1021 
1022 template <class S, typename T> inline
1023 S& operator <<( S& s, const Complex<T>& c )
1024 {
1025  return s << c.Real() << c.Imag();
1026 }
1027  */
1028 
1029 // ----------------------------------------------------------------------------
1030 
1031 #ifndef __PCL_NO_COMPLEX_INSTANTIATE
1032 
1045 typedef Complex<float> Complex32;
1046 
1056 typedef Complex32 fcomplex;
1057 
1066 typedef Complex<double> Complex64;
1067 
1077 typedef Complex64 dcomplex;
1078 
1088 typedef dcomplex complex;
1089 
1090 #endif // __PCL_NO_COMPLEX_INSTANTIATE
1091 
1092 // ----------------------------------------------------------------------------
1093 
1094 } // pcl
1095 
1096 #endif // __PCL_Complex_h
1097 
1098 // ----------------------------------------------------------------------------
1099 // EOF pcl/Complex.h - Released 2019-04-30T16:30:41Z
Complex< T > Log(const Complex< T > &c)
Definition: Complex.h:726
static constexpr T log10e()
Definition: Constants.h:159
void CrossPowerSpectrumMatrix(Complex< T > *i, const Complex< T > *j, const Complex< T > *a, const Complex< T > *b)
Definition: Complex.h:1005
Complex< T1 > operator*(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:539
constexpr T ArcTan(T x)
Definition: Math.h:481
Complex< T1 > operator+(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:455
A complex number whose components are 32-bit floating point real numbers.
PCL root namespace.
Definition: AbstractImage.h:76
double Norm(const T *i, const T *j, double p)
Definition: Math.h:1922
Complex< T > Polar(T r, T stheta, T ctheta)
Definition: Complex.h:432
T Mag() const
Definition: Complex.h:343
Complex< T > Sqrt(const Complex< T > &c)
Definition: Complex.h:665
Complex< T > Round(const Complex< T > &c)
Definition: Complex.h:929
Complex(const Complex< T1 > &c)
Definition: Complex.h:113
A complex number whose components are 64-bit floating point real numbers.
T Abs(const Complex< T > &c)
Definition: Complex.h:420
Complex< T > Conj() const
Definition: Complex.h:315
constexpr T Arg() const
Definition: Complex.h:375
constexpr const T & Max(const T &a, const T &b)
Definition: Utility.h:119
Complex(T r, T i=0)
Definition: Complex.h:103
Generic complex number.
Definition: Complex.h:83
void SetConj()
Definition: Complex.h:332
Complex< T > Sinh(const Complex< T > &c)
Definition: Complex.h:818
Complex< T1 > operator-(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:495
constexpr T Imag() const
Definition: Complex.h:138
A complex number whose components are 64-bit floating point real numbers.
Complex< T > Cosh(const Complex< T > &c)
Definition: Complex.h:829
constexpr bool IsReal() const
Definition: Complex.h:154
Complex< T > Ln(const Complex< T > &c)
Definition: Complex.h:716
void PhaseCorrelationMatrix(Complex< T > *i, const Complex< T > *j, const Complex< T > *a, const Complex< T > *b)
Definition: Complex.h:974
A complex number whose components are 64-bit floating point real numbers.
bool operator==(const Array< T, A > &x1, const Array< T, A > &x2)
Definition: Array.h:2075
Complex< T1 > operator/(const Complex< T1 > &c1, const Complex< T2 > &c2)
Definition: Complex.h:583
constexpr T Real() const
Definition: Complex.h:122
Complex< T > Sin(const Complex< T > &c)
Definition: Complex.h:786
Complex< T > Cos(const Complex< T > &c)
Definition: Complex.h:797
Complex< T > Exp(const Complex< T > &c)
Definition: Complex.h:705
A complex number whose components are 32-bit floating point real numbers.
Complex< T > Tan(const Complex< T > &c)
Definition: Complex.h:808
Complex< T1 > Pow(const Complex< T1 > &c, T2 x)
Definition: Complex.h:738
Complex< T > Tanh(const Complex< T > &c)
Definition: Complex.h:840