PCL
GridInterpolation.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.8.5
6 // ----------------------------------------------------------------------------
7 // pcl/GridInterpolation.h - Released 2024-12-28T16:53:48Z
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-2024 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 (https://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_GridInterpolation_h
53 #define __PCL_GridInterpolation_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
60 #include <pcl/AbstractImage.h>
61 #include <pcl/BicubicInterpolation.h>
62 #include <pcl/Matrix.h>
63 #include <pcl/ParallelProcess.h>
64 #include <pcl/Rectangle.h>
65 #include <pcl/ReferenceArray.h>
66 #include <pcl/Thread.h>
67 
68 
69 #include <pcl/File.h>
70 #include <pcl/TimePoint.h>
71 
72 
73 
74 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
75 # include <pcl/Console.h>
76 # include <pcl/StandardStatus.h>
77 #endif
78 
79 namespace pcl
80 {
81 
82 // ----------------------------------------------------------------------------
83 
97 class PCL_CLASS GridInterpolation : public ParallelProcess
98 {
99 public:
100 
105  GridInterpolation() = default;
106 
110  GridInterpolation( const GridInterpolation& ) = default;
111 
115  GridInterpolation( GridInterpolation&& ) = default;
116 
120  GridInterpolation& operator =( const GridInterpolation& ) = default;
121 
125  GridInterpolation& operator =( GridInterpolation&& ) = default;
126 
188  template <class R, class SI>
189  void Initialize( const R& rect, int delta, const SI& S, bool verbose = true )
190  {
191  PCL_PRECONDITION( rect.IsRect() )
192  PCL_PRECONDITION( delta > 0 )
193 
194  m_rect = rect.Ordered();
195  double width = m_rect.Width();
196  double height = m_rect.Height();
197  if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
198  throw Error( "GridInterpolation::Initialize(): Empty interpolation space." );
199 
200  m_delta = Abs( delta );
201  if ( 1 + m_delta == 1 )
202  throw Error( "GridInterpolation::Initialize(): Zero or insignificant grid distance." );
203 
204  int rows = 1 + int( Ceil( height/m_delta ) );
205  int cols = 1 + int( Ceil( width/m_delta ) );
206  m_G = DMatrix( rows, cols );
207 
208  StatusMonitor monitor;
209 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
210  StandardStatus status;
211  if ( verbose )
212  {
213  monitor.SetCallback( &status );
214  monitor.Initialize( "Building surface interpolation grid", rows );
215  }
216 #endif
217 
218  Array<size_type> L = Thread::OptimalThreadLoads( rows,
219  1/*overheadLimit*/,
220  m_parallel ? m_maxProcessors : 1 );
221  AbstractImage::ThreadData data( monitor, rows );
222  ReferenceArray<GridInitializationThread<SI> > threads;
223  for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
224  threads.Add( new GridInitializationThread<SI>( data, *this, S, n, n + int( L[i] ) ) );
225  AbstractImage::RunThreads( threads, data );
226  threads.Destroy();
227 
228  m_I.Initialize( m_G.Begin(), cols, rows );
229  }
230 
261  template <class R>
262  void Initialize( const R& rect, double delta, const DMatrix& G )
263  {
264  PCL_PRECONDITION( rect.IsRect() )
265  PCL_PRECONDITION( delta > 0 )
266 
267  m_rect = DRect( rect.Ordered() );
268  double width = m_rect.Width();
269  double height = m_rect.Height();
270  if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
271  throw Error( "GridInterpolation::Initialize(): Empty interpolation space." );
272 
273  m_delta = Abs( delta );
274  if ( 1 + m_delta == 1 )
275  throw Error( "GridInterpolation::Initialize(): Zero or insignificant grid distance." );
276 
277  int rows = 1 + int( Ceil( height/m_delta ) );
278  int cols = 1 + int( Ceil( width/m_delta ) );
279  if ( G.Rows() != rows || G.Cols() != cols )
280  throw Error( "GridInterpolation::Initialize(): Invalid matrix dimensions." );
281 
282  m_G = G;
283  m_I.Initialize( m_G.Begin(), cols, rows );
284  }
285 
290  bool IsValid() const
291  {
292  return !m_G.IsEmpty();
293  }
294 
299  void Clear()
300  {
301  m_I.Clear();
302  m_G.Clear();
303  }
304 
309  const DRect& ReferenceRect() const
310  {
311  return m_rect;
312  }
313 
318  double Delta() const
319  {
320  return m_delta;
321  }
322 
330  const DMatrix& InterpolationMatrix() const
331  {
332  return m_G;
333  }
334 
338  template <typename T>
339  double operator ()( T x, T y ) const
340  {
341  double fx = (double( x ) - m_rect.x0)/m_delta;
342  double fy = (double( y ) - m_rect.y0)/m_delta;
343  return m_I( fx, fy );
344  }
345 
349  template <typename T>
350  double operator ()( const GenericPoint<T>& p ) const
351  {
352  return operator ()( p.x, p.y );
353  }
354 
355 private:
356 
362  using grid_interpolation = BicubicBSplineInterpolation<double>;
363 
364  DRect m_rect;
365  double m_delta;
366  DMatrix m_G;
367  grid_interpolation m_I;
368 
369  template <class SI>
370  class GridInitializationThread : public Thread
371  {
372  public:
373 
374  GridInitializationThread( const AbstractImage::ThreadData& data,
375  GridInterpolation& grid, const SI& surface, int startRow, int endRow )
376  : m_data( data )
377  , m_grid( grid )
378  , m_surface( surface )
379  , m_startRow( startRow )
380  , m_endRow( endRow )
381  {
382  if ( m_surface.HasFastVectorEvaluation() )
383  {
384  m_N = int( Min( size_type( 4096 ), size_type( m_endRow - m_startRow )*m_grid.m_G.Cols() ) );
385  m_X = DVector( m_N );
386  m_Y = DVector( m_N );
387  m_Z = DVector( m_N );
388  }
389  }
390 
391  PCL_HOT_FUNCTION void Run() override
392  {
393  INIT_THREAD_MONITOR()
394 
395  if ( m_surface.HasFastVectorEvaluation() )
396  {
397  int i0 = m_startRow, j0 = 0, n = 0;
398  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
399  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
400  {
401  double x = m_grid.m_rect.x0;
402  for ( int j = 0; j < m_grid.m_G.Cols(); ++j, ++n, x += m_grid.m_delta )
403  {
404  if ( n == m_N )
405  {
406  m_surface.Evaluate( m_Z.Begin(), m_X.Begin(), m_Y.Begin(), n );
407  for ( int i1 = i0, k1 = 0; i1 <= i && k1 < n; ++i1, j0 = 0 )
408  for ( int j1 = j0; j1 < m_grid.m_G.Cols() && k1 < n; ++j1, ++k1 )
409  m_grid.m_G[i1][j1] = m_Z[k1];
410  i0 = i;
411  j0 = j;
412  n = 0;
413  }
414 
415  m_X[n] = x;
416  m_Y[n] = y;
417  }
418 
419  UPDATE_THREAD_MONITOR( 1 )
420  }
421 
422  if ( n > 0 )
423  {
424  m_surface.Evaluate( m_Z.Begin(), m_X.Begin(), m_Y.Begin(), n );
425  for ( int i1 = i0, k1 = 0; i1 < m_endRow; ++i1, j0 = 0 )
426  for ( int j1 = j0; j1 < m_grid.m_G.Cols(); ++j1, ++k1 )
427  m_grid.m_G[i1][j1] = m_Z[k1];
428 
429  UPDATE_THREAD_MONITOR( 1 )
430  }
431  }
432  else
433  {
434  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
435  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
436  {
437  double x = m_grid.m_rect.x0;
438  for ( int j = 0; j < m_grid.m_G.Cols(); ++j, x += m_grid.m_delta )
439  m_grid.m_G[i][j] = m_surface( x, y );
440 
441  UPDATE_THREAD_MONITOR( 32 )
442  }
443  }
444  }
445 
446  private:
447 
448  const AbstractImage::ThreadData& m_data;
449  GridInterpolation& m_grid;
450  const SI& m_surface;
451  int m_startRow, m_endRow, m_N;
452  DVector m_X, m_Y, m_Z;
453  };
454 };
455 
456 // ----------------------------------------------------------------------------
457 
471 class PCL_CLASS PointGridInterpolation : public ParallelProcess
472 {
473 public:
474 
479  PointGridInterpolation() = default;
480 
484  PointGridInterpolation( const PointGridInterpolation& ) = default;
485 
489  PointGridInterpolation( PointGridInterpolation&& ) = default;
490 
494  PointGridInterpolation& operator =( const PointGridInterpolation& ) = default;
495 
499  PointGridInterpolation& operator =( PointGridInterpolation&& ) = default;
500 
563  template <class R, class PSI>
564  void Initialize( const R& rect, double delta, const PSI& PS, bool verbose = true )
565  {
566  PCL_PRECONDITION( rect.IsRect() )
567  PCL_PRECONDITION( delta > 0 )
568 
569  m_rect = DRect( rect.Ordered() );
570  double width = m_rect.Width();
571  double height = m_rect.Height();
572  if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
573  throw Error( "PointGridInterpolation::Initialize(): Empty interpolation space." );
574 
575  m_delta = Abs( delta );
576  if ( 1 + m_delta == 1 )
577  throw Error( "PointGridInterpolation::Initialize(): Zero or insignificant grid distance." );
578 
579  int rows = 1 + int( Ceil( height/m_delta ) );
580  int cols = 1 + int( Ceil( width/m_delta ) );
581  m_Gx = DMatrix( rows, cols );
582  m_Gy = DMatrix( rows, cols );
583 
584  StatusMonitor monitor;
585 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
586  StandardStatus status;
587  if ( verbose )
588  {
589  monitor.SetCallback( &status );
590  monitor.Initialize( "Building surface interpolation grid", rows );
591  }
592 #endif
593 
594  Array<size_type> L = Thread::OptimalThreadLoads( rows,
595  1/*overheadLimit*/,
596  m_parallel ? m_maxProcessors : 1 );
597  AbstractImage::ThreadData data( monitor, rows );
598  ReferenceArray<GridInitializationThread<PSI> > threads;
599  for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
600  threads.Add( new GridInitializationThread<PSI>( data, *this, PS, n, n + int( L[i] ) ) );
601  AbstractImage::RunThreads( threads, data );
602  threads.Destroy();
603 
604  m_Ix.Initialize( m_Gx.Begin(), cols, rows );
605  m_Iy.Initialize( m_Gy.Begin(), cols, rows );
606  }
607 
676  template <class R, class SI>
677  void Initialize( const R& rect, double delta, const SI& Sx, const SI& Sy, bool verbose = true )
678  {
679  PCL_PRECONDITION( rect.IsRect() )
680  PCL_PRECONDITION( delta > 0 )
681 
682  m_rect = DRect( rect.Ordered() );
683  double width = m_rect.Width();
684  double height = m_rect.Height();
685  if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
686  throw Error( "PointGridInterpolation::Initialize(): Empty interpolation space." );
687 
688  m_delta = Abs( delta );
689  if ( 1 + m_delta == 1 )
690  throw Error( "PointGridInterpolation::Initialize(): Zero or insignificant grid distance." );
691 
692  int rows = 1 + int( Ceil( height/m_delta ) );
693  int cols = 1 + int( Ceil( width/m_delta ) );
694  m_Gx = DMatrix( rows, cols );
695  m_Gy = DMatrix( rows, cols );
696 
697  StatusMonitor monitor;
698 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
699  StandardStatus status;
700  if ( verbose )
701  {
702  monitor.SetCallback( &status );
703  monitor.Initialize( "Building surface interpolation grid", rows );
704  }
705 #endif
706 
707  Array<size_type> L = Thread::OptimalThreadLoads( rows,
708  1/*overheadLimit*/,
709  m_parallel ? m_maxProcessors : 1 );
710  AbstractImage::ThreadData data( monitor, rows );
711  ReferenceArray<GridInitializationXYThread<SI> > threads;
712  for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
713  threads.Add( new GridInitializationXYThread<SI>( data, *this, Sx, Sy, n, n + int( L[i] ) ) );
714  AbstractImage::RunThreads( threads, data );
715  threads.Destroy();
716 
717  m_Ix.Initialize( m_Gx.Begin(), cols, rows );
718  m_Iy.Initialize( m_Gy.Begin(), cols, rows );
719  }
720 
753  template <class R>
754  void Initialize( const R& rect, int delta, const DMatrix& Gx, const DMatrix& Gy )
755  {
756  PCL_PRECONDITION( rect.IsRect() )
757  PCL_PRECONDITION( delta > 0 )
758 
759  m_rect = DRect( rect.Ordered() );
760  double width = m_rect.Width();
761  double height = m_rect.Height();
762  if ( !m_rect.IsRect() || 1 + width == 1 || 1 + height == 1 )
763  throw Error( "PointGridInterpolation::Initialize(): Empty interpolation region." );
764 
765  m_delta = Abs( delta );
766  if ( 1 + m_delta == 1 )
767  throw Error( "PointGridInterpolation::Initialize(): Zero or insignificant grid distance." );
768 
769  int rows = 1 + int( Ceil( height/m_delta ) );
770  int cols = 1 + int( Ceil( width/m_delta ) );
771  if ( Gx.Rows() != rows || Gx.Cols() != cols || Gy.Rows() != rows || Gy.Cols() != cols )
772  throw Error( "PointGridInterpolation::Initialize(): Invalid matrix dimensions." );
773 
774  m_Gx = Gx;
775  m_Gy = Gy;
776  m_Ix.Initialize( m_Gx.Begin(), cols, rows );
777  m_Iy.Initialize( m_Gy.Begin(), cols, rows );
778  }
779 
783  template <class PSI>
784  void ApplyLocalModel( const PSI& PS,
785  const String& message = "Applying local interpolation model",
786  bool verbose = true )
787  {
788  PCL_PRECONDITION( IsValid() )
789  if ( !IsValid() )
790  throw Error( "PointGridInterpolation::ApplyLocalModel(): Uninitialized interpolation." );
791 
792  StatusMonitor monitor;
793 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
794  StandardStatus status;
795  if ( verbose )
796  {
797  monitor.SetCallback( &status );
798  monitor.Initialize( message, m_Gx.Rows() );
799  }
800 #endif
801 
802  Array<size_type> L = Thread::OptimalThreadLoads( m_Gx.Rows(),
803  1/*overheadLimit*/,
804  m_parallel ? m_maxProcessors : 1 );
805  AbstractImage::ThreadData data( monitor, m_Gx.Rows() );
806  ReferenceArray<LocalModelThread<PSI> > threads;
807  for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
808  threads.Add( new LocalModelThread<PSI>( data, *this, PS, n, n + int( L[i] ) ) );
809  AbstractImage::RunThreads( threads, data );
810  threads.Destroy();
811  }
812 
816  template <class SI>
817  void ApplyLocalModel( const SI& Sx, const SI& Sy,
818  const String& message = "Applying local interpolation model",
819  bool verbose = true )
820  {
821  PCL_PRECONDITION( IsValid() )
822  if ( !IsValid() )
823  throw Error( "PointGridInterpolation::ApplyLocalModel(): Uninitialized interpolation." );
824 
825  StatusMonitor monitor;
826 #ifndef __PCL_BUILDING_PIXINSIGHT_APPLICATION
827  StandardStatus status;
828  if ( verbose )
829  {
830  monitor.SetCallback( &status );
831  monitor.Initialize( message, m_Gx.Rows() );
832  }
833 #endif
834 
835  Array<size_type> L = Thread::OptimalThreadLoads( m_Gx.Rows(),
836  1/*overheadLimit*/,
837  m_parallel ? m_maxProcessors : 1 );
838  AbstractImage::ThreadData data( monitor, m_Gx.Rows() );
839  ReferenceArray<LocalModelThread2<SI> > threads;
840  for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
841  threads.Add( new LocalModelThread2<SI>( data, *this, Sx, Sy, n, n + int( L[i] ) ) );
842  AbstractImage::RunThreads( threads, data );
843  threads.Destroy();
844  }
845 
850  bool IsValid() const
851  {
852  return !m_Gx.IsEmpty() && !m_Gy.IsEmpty();
853  }
854 
859  void Clear()
860  {
861  m_Ix.Clear();
862  m_Iy.Clear();
863  m_Gx.Clear();
864  m_Gy.Clear();
865  }
866 
873  const DRect& ReferenceRect() const
874  {
875  return m_rect;
876  }
877 
882  double Delta() const
883  {
884  return m_delta;
885  }
886 
894  const DMatrix& XInterpolationMatrix() const
895  {
896  return m_Gx;
897  }
898 
906  const DMatrix& YInterpolationMatrix() const
907  {
908  return m_Gy;
909  }
910 
914  template <typename T>
915  DPoint operator ()( T x, T y ) const
916  {
917  PCL_PRECONDITION( IsValid() )
918  double fx = (double( x ) - m_rect.x0)/m_delta;
919  double fy = (double( y ) - m_rect.y0)/m_delta;
920  return DPoint( m_Ix( fx, fy ), m_Iy( fx, fy ) );
921  }
922 
926  template <typename T>
927  DPoint operator ()( const GenericPoint<T>& p ) const
928  {
929  return operator ()( p.x, p.y );
930  }
931 
932 private:
933 
939  using grid_interpolation = BicubicBSplineInterpolation<double>;
940 
941  DRect m_rect;
942  double m_delta;
943  DMatrix m_Gx, m_Gy;
944  grid_interpolation m_Ix, m_Iy;
945 
946  template <class PSI>
947  class GridInitializationThread : public Thread
948  {
949  public:
950 
951  GridInitializationThread( const AbstractImage::ThreadData& data,
952  PointGridInterpolation& grid, const PSI& surface, int startRow, int endRow )
953  : m_data( data )
954  , m_grid( grid )
955  , m_surface( surface )
956  , m_startRow( startRow )
957  , m_endRow( endRow )
958  {
959  if ( m_surface.HasFastVectorEvaluation() )
960  {
961  m_N = int( Min( size_type( 4096 ), size_type( m_endRow - m_startRow )*m_grid.m_Gx.Cols() ) );
962  m_X = DVector( m_N );
963  m_Y = DVector( m_N );
964  m_ZX = DVector( m_N );
965  m_ZY = DVector( m_N );
966  }
967  }
968 
969  PCL_HOT_FUNCTION void Run() override
970  {
971  INIT_THREAD_MONITOR()
972 
973  if ( m_surface.HasFastVectorEvaluation() )
974  {
975  int i0 = m_startRow, j0 = 0, n = 0;
976  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
977  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
978  {
979  double x = m_grid.m_rect.x0;
980  for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, ++n, x += m_grid.m_delta )
981  {
982  if ( n == m_N )
983  {
984  m_surface.Evaluate( m_ZX.Begin(), m_ZY.Begin(), m_X.Begin(), m_Y.Begin(), n );
985  for ( int i1 = i0, k1 = 0; i1 <= i && k1 < n; ++i1, j0 = 0 )
986  for ( int j1 = j0; j1 < m_grid.m_Gx.Cols() && k1 < n; ++j1, ++k1 )
987  {
988  m_grid.m_Gx[i1][j1] = m_ZX[k1];
989  m_grid.m_Gy[i1][j1] = m_ZY[k1];
990  }
991  i0 = i;
992  j0 = j;
993  n = 0;
994  }
995 
996  m_X[n] = x;
997  m_Y[n] = y;
998  }
999 
1000  UPDATE_THREAD_MONITOR( 1 )
1001  }
1002 
1003  if ( n > 0 )
1004  {
1005  m_surface.Evaluate( m_ZX.Begin(), m_ZY.Begin(), m_X.Begin(), m_Y.Begin(), n );
1006  for ( int i1 = i0, k1 = 0; i1 < m_endRow; ++i1, j0 = 0 )
1007  for ( int j1 = j0; j1 < m_grid.m_Gx.Cols(); ++j1, ++k1 )
1008  {
1009  m_grid.m_Gx[i1][j1] = m_ZX[k1];
1010  m_grid.m_Gy[i1][j1] = m_ZY[k1];
1011  }
1012 
1013  UPDATE_THREAD_MONITOR( 1 )
1014  }
1015  }
1016  else
1017  {
1018  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
1019  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
1020  {
1021  double x = m_grid.m_rect.x0;
1022  for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
1023  {
1024  DPoint p = m_surface( x, y );
1025  m_grid.m_Gx[i][j] = p.x;
1026  m_grid.m_Gy[i][j] = p.y;
1027  }
1028 
1029  UPDATE_THREAD_MONITOR( 32 )
1030  }
1031  }
1032  }
1033 
1034  private:
1035 
1036  const AbstractImage::ThreadData& m_data;
1037  PointGridInterpolation& m_grid;
1038  const PSI& m_surface;
1039  int m_startRow, m_endRow, m_N;
1040  DVector m_X, m_Y, m_ZX, m_ZY;
1041  };
1042 
1043  template <class SI>
1044  class GridInitializationXYThread : public Thread
1045  {
1046  public:
1047 
1048  GridInitializationXYThread( const AbstractImage::ThreadData& data,
1049  PointGridInterpolation& grid,
1050  const SI& surfaceX, const SI& surfaceY, int startRow, int endRow )
1051  : m_data( data )
1052  , m_grid( grid )
1053  , m_surfaceX( surfaceX )
1054  , m_surfaceY( surfaceY )
1055  , m_startRow( startRow )
1056  , m_endRow( endRow )
1057  {
1058  if ( m_surfaceX.HasFastVectorEvaluation() && m_surfaceY.HasFastVectorEvaluation() )
1059  {
1060  m_N = int( Min( size_type( 4096 ), size_type( m_endRow - m_startRow )*m_grid.m_Gx.Cols() ) );
1061  m_X = DVector( m_N );
1062  m_Y = DVector( m_N );
1063  m_Z = DVector( m_N );
1064  }
1065  }
1066 
1067  PCL_HOT_FUNCTION void Run() override
1068  {
1069  INIT_THREAD_MONITOR()
1070 
1071  if ( m_surfaceX.HasFastVectorEvaluation() && m_surfaceY.HasFastVectorEvaluation() )
1072  {
1073  int i0 = m_startRow, j0 = 0, n = 0;
1074  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
1075  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
1076  {
1077  double x = m_grid.m_rect.x0;
1078  for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, ++n, x += m_grid.m_delta )
1079  {
1080  if ( n == m_N )
1081  {
1082  m_surfaceX.Evaluate( m_Z.Begin(), m_X.Begin(), m_Y.Begin(), n );
1083  for ( int i1 = i0, jj0 = j0, k1 = 0; i1 <= i && k1 < n; ++i1, jj0 = 0 )
1084  for ( int j1 = jj0; j1 < m_grid.m_Gx.Cols() && k1 < n; ++j1, ++k1 )
1085  m_grid.m_Gx[i1][j1] = m_Z[k1];
1086  m_surfaceY.Evaluate( m_Z.Begin(), m_X.Begin(), m_Y.Begin(), n );
1087  for ( int i1 = i0, jj0 = j0, k1 = 0; i1 <= i && k1 < n; ++i1, jj0 = 0 )
1088  for ( int j1 = jj0; j1 < m_grid.m_Gx.Cols() && k1 < n; ++j1, ++k1 )
1089  m_grid.m_Gy[i1][j1] = m_Z[k1];
1090  i0 = i;
1091  j0 = j;
1092  n = 0;
1093  }
1094 
1095  m_X[n] = x;
1096  m_Y[n] = y;
1097  }
1098 
1099  UPDATE_THREAD_MONITOR( 1 )
1100  }
1101 
1102  if ( n > 0 )
1103  {
1104  m_surfaceX.Evaluate( m_Z.Begin(), m_X.Begin(), m_Y.Begin(), n );
1105  for ( int i1 = i0, jj0 = j0, k1 = 0; i1 < m_endRow; ++i1, jj0 = 0 )
1106  for ( int j1 = jj0; j1 < m_grid.m_Gx.Cols(); ++j1, ++k1 )
1107  m_grid.m_Gx[i1][j1] = m_Z[k1];
1108  m_surfaceY.Evaluate( m_Z.Begin(), m_X.Begin(), m_Y.Begin(), n );
1109  for ( int i1 = i0, jj0 = j0, k1 = 0; i1 < m_endRow; ++i1, jj0 = 0 )
1110  for ( int j1 = jj0; j1 < m_grid.m_Gx.Cols(); ++j1, ++k1 )
1111  m_grid.m_Gy[i1][j1] = m_Z[k1];
1112 
1113  UPDATE_THREAD_MONITOR( 1 )
1114  }
1115  }
1116  else
1117  {
1118  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
1119  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
1120  {
1121  double x = m_grid.m_rect.x0;
1122  for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
1123  {
1124  m_grid.m_Gx[i][j] = m_surfaceX( x, y );
1125  m_grid.m_Gy[i][j] = m_surfaceY( x, y );
1126  }
1127 
1128  UPDATE_THREAD_MONITOR( 32 )
1129  }
1130  }
1131  }
1132 
1133  private:
1134 
1135  const AbstractImage::ThreadData& m_data;
1136  PointGridInterpolation& m_grid;
1137  const SI& m_surfaceX;
1138  const SI& m_surfaceY;
1139  int m_startRow, m_endRow, m_N;
1140  DVector m_X, m_Y, m_Z;
1141  };
1142 
1143  template <class PSI>
1144  class LocalModelThread : public Thread
1145  {
1146  public:
1147 
1148  LocalModelThread( const AbstractImage::ThreadData& data,
1149  PointGridInterpolation& grid, const PSI& model, int startRow, int endRow )
1150  : m_data( data )
1151  , m_grid( grid )
1152  , m_model( model )
1153  , m_startRow( startRow )
1154  , m_endRow( endRow )
1155  {
1156  }
1157 
1158  PCL_HOT_FUNCTION void Run() override
1159  {
1160  INIT_THREAD_MONITOR()
1161 
1162  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
1163  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
1164  {
1165  double x = m_grid.m_rect.x0;
1166  for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
1167  {
1168  DPoint d = m_model( x, y );
1169  m_grid.m_Gx[i][j] += d.x;
1170  m_grid.m_Gy[i][j] += d.y;
1171  }
1172 
1173  UPDATE_THREAD_MONITOR( 32 )
1174  }
1175  }
1176 
1177  private:
1178 
1179  const AbstractImage::ThreadData& m_data;
1180  PointGridInterpolation& m_grid;
1181  const PSI& m_model;
1182  int m_startRow, m_endRow;
1183  };
1184 
1185  template <class SI>
1186  class LocalModelThread2 : public Thread
1187  {
1188  public:
1189 
1190  LocalModelThread2( const AbstractImage::ThreadData& data,
1191  PointGridInterpolation& grid, const SI& modelX, const SI& modelY, int startRow, int endRow )
1192  : m_data( data )
1193  , m_grid( grid )
1194  , m_modelX( modelX )
1195  , m_modelY( modelY )
1196  , m_startRow( startRow )
1197  , m_endRow( endRow )
1198  {
1199  }
1200 
1201  PCL_HOT_FUNCTION void Run() override
1202  {
1203  INIT_THREAD_MONITOR()
1204 
1205  double y = m_grid.m_rect.y0 + m_startRow*m_grid.m_delta;
1206  for ( int i = m_startRow; i < m_endRow; ++i, y += m_grid.m_delta )
1207  {
1208  double x = m_grid.m_rect.x0;
1209  for ( int j = 0; j < m_grid.m_Gx.Cols(); ++j, x += m_grid.m_delta )
1210  {
1211  m_grid.m_Gx[i][j] += m_modelX( x, y );
1212  m_grid.m_Gy[i][j] += m_modelY( x, y );
1213  }
1214 
1215  UPDATE_THREAD_MONITOR( 32 )
1216  }
1217  }
1218 
1219  private:
1220 
1221  const AbstractImage::ThreadData& m_data;
1222  PointGridInterpolation& m_grid;
1223  const SI& m_modelX, m_modelY;
1224  int m_startRow, m_endRow;
1225  };
1226 
1227  friend class SplineWorldTransformation;
1228 };
1229 
1230 // ----------------------------------------------------------------------------
1231 
1232 } // pcl
1233 
1234 #endif // __PCL_GridInterpolation_h
1235 
1236 // ----------------------------------------------------------------------------
1237 // EOF pcl/GridInterpolation.h - Released 2024-12-28T16:53:48Z
pcl::AbstractImage::RunThreads
static void RunThreads(ReferenceArray< thread > &threads, ThreadData &data, bool useAffinity=true)
Definition: AbstractImage.h:1072
pcl::Array
Generic dynamic array.
Definition: Array.h:100
DPoint
64-bit floating-point point in the R^2 space.
DVector
64-bit floating point real vector.
pcl::Error
A simple exception with an associated error message.
Definition: Exception.h:239
pcl::GenericMatrix
Generic dynamic matrix of arbitrary dimensions.
Definition: Matrix.h:123
pcl::GenericMatrix::Cols
int Cols() const noexcept
Definition: Matrix.h:987
pcl::GenericMatrix::Rows
int Rows() const noexcept
Definition: Matrix.h:978
pcl::GenericPoint
A generic point in the two-dimensional space.
Definition: Point.h:100
pcl::GenericPoint::x
component x
Abscissa (horizontal, or X-axis coordinate).
Definition: Point.h:111
pcl::GenericPoint::y
component y
Ordinate (vertical, or Y-axis coordinate).
Definition: Point.h:112
pcl::GenericRectangle
A generic rectangle in the two-dimensional space.
Definition: Rectangle.h:316
pcl::GenericRectangle::IsRect
bool IsRect() const noexcept
Definition: Rectangle.h:764
pcl::GenericRectangle::y0
component y0
Vertical coordinate of the upper left corner.
Definition: Rectangle.h:335
pcl::GenericRectangle::x0
component x0
Horizontal coordinate of the upper left corner.
Definition: Rectangle.h:334
pcl::GenericRectangle::Width
component Width() const noexcept
Definition: Rectangle.h:637
pcl::GenericRectangle::Height
component Height() const noexcept
Definition: Rectangle.h:646
pcl::GenericVector
Generic vector of arbitrary length.
Definition: Vector.h:107
pcl::GridInterpolation
Discretized scalar surface interpolation/approximation in two dimensions.
Definition: GridInterpolation.h:98
pcl::GridInterpolation::GridInterpolation
GridInterpolation(GridInterpolation &&)=default
pcl::GridInterpolation::Initialize
void Initialize(const R &rect, double delta, const DMatrix &G)
Definition: GridInterpolation.h:262
pcl::GridInterpolation::GridInterpolation
GridInterpolation(const GridInterpolation &)=default
pcl::GridInterpolation::InterpolationMatrix
const DMatrix & InterpolationMatrix() const
Definition: GridInterpolation.h:330
pcl::GridInterpolation::ReferenceRect
const DRect & ReferenceRect() const
Definition: GridInterpolation.h:309
pcl::GridInterpolation::Initialize
void Initialize(const R &rect, int delta, const SI &S, bool verbose=true)
Definition: GridInterpolation.h:189
pcl::ParallelProcess
A process using multiple concurrent execution threads.
Definition: ParallelProcess.h:73
pcl::PointGridInterpolation
Discretized vector surface interpolation/approximation in two dimensions.
Definition: GridInterpolation.h:472
pcl::PointGridInterpolation::Initialize
void Initialize(const R &rect, double delta, const SI &Sx, const SI &Sy, bool verbose=true)
Definition: GridInterpolation.h:677
pcl::PointGridInterpolation::XInterpolationMatrix
const DMatrix & XInterpolationMatrix() const
Definition: GridInterpolation.h:894
pcl::PointGridInterpolation::Initialize
void Initialize(const R &rect, double delta, const PSI &PS, bool verbose=true)
Definition: GridInterpolation.h:564
pcl::PointGridInterpolation::PointGridInterpolation
PointGridInterpolation(const PointGridInterpolation &)=default
pcl::PointGridInterpolation::YInterpolationMatrix
const DMatrix & YInterpolationMatrix() const
Definition: GridInterpolation.h:906
pcl::PointGridInterpolation::Initialize
void Initialize(const R &rect, int delta, const DMatrix &Gx, const DMatrix &Gy)
Definition: GridInterpolation.h:754
pcl::PointGridInterpolation::PointGridInterpolation
PointGridInterpolation(PointGridInterpolation &&)=default
pcl::PointGridInterpolation::ReferenceRect
const DRect & ReferenceRect() const
Definition: GridInterpolation.h:873
pcl::ReferenceArray
Dynamic array of pointers to objects providing direct iteration and element access by reference.
Definition: ReferenceArray.h:96
pcl::ReferenceArray::Destroy
void Destroy(iterator i, size_type n=1)
Definition: ReferenceArray.h:1518
pcl::ReferenceArray::Add
void Add(const ReferenceArray &x)
Definition: ReferenceArray.h:1345
pcl::StandardStatus
A status monitoring callback that sends progress information to the process console.
Definition: StandardStatus.h:100
pcl::StatusMonitor
An asynchronous status monitoring system.
Definition: StatusMonitor.h:223
pcl::StatusMonitor::SetCallback
void SetCallback(StatusCallback *callback)
Definition: StatusMonitor.h:499
pcl::StatusMonitor::Initialize
void Initialize(const String &info, size_type count=0)
pcl::String
Unicode (UTF-16) string.
Definition: String.h:8146
pcl::Thread
Client-side interface to a PixInsight thread.
Definition: Thread.h:126
pcl::Thread::OptimalThreadLoads
static Array< size_type > OptimalThreadLoads(size_type count, size_type overheadLimit=1u, int maxThreads=PCL_MAX_PROCESSORS)
pcl::Abs
T Abs(const Complex< T > &c) noexcept
Definition: Complex.h:443
pcl::Ceil
constexpr T Ceil(T x) noexcept
Definition: Math.h:633
pcl::size_type
size_t size_type
Definition: Defs.h:609
INIT_THREAD_MONITOR
#define INIT_THREAD_MONITOR()
Declares and initializes local variables used for synchronization of thread status monitoring.
Definition: AbstractImage.h:1226
UPDATE_THREAD_MONITOR
#define UPDATE_THREAD_MONITOR(N)
Synchronized status monitoring of a set of image processing threads.
Definition: AbstractImage.h:1314
pcl::Min
constexpr const T & Min(const T &a, const T &b) noexcept
Definition: Utility.h:90
pcl
PCL root namespace.
Definition: AbstractImage.h:77
pcl::AbstractImage::ThreadData
Thread synchronization data for status monitoring of parallel image processing tasks.
Definition: AbstractImage.h:965