PCL
StarDatabaseFile.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.7.0
6 // ----------------------------------------------------------------------------
7 // pcl/StarDatabaseFile.h - Released 2024-06-18T15:48:54Z
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_StarDatabaseFile_h
53 #define __PCL_StarDatabaseFile_h
54 
56 
57 #include <pcl/Defs.h>
58 
59 #include <pcl/AutoPointer.h>
60 #include <pcl/Compression.h>
61 #include <pcl/ElapsedTime.h>
62 #include <pcl/File.h>
63 #include <pcl/TimePoint.h>
64 #include <pcl/Vector.h>
65 
66 namespace pcl
67 {
68 
69 // ----------------------------------------------------------------------------
70 
75 // ----------------------------------------------------------------------------
76 
77 class PCL_CLASS StarDatabaseFile;
78 
88 class PCL_CLASS XPSD
89 {
90 public:
91 
111  struct Metadata
112  {
124  };
125 
136  struct Statistics
137  {
138  uint64 totalSources = 0;
139  uint32 totalNodes = 0;
140  uint32 totalLeaves = 0;
141  float medianLeafLength = 0;
142  uint32 minimumLeafLength = 0;
143  uint32 maximumLeafLength = 0;
144  };
145 
154  {
155  double centerRA = 0;
156  double centerDec = 0;
157  double radius = 1;
158  float magnitudeLow = -1.5;
160  float magnitudeHigh = 26;
162  uint32 sourceLimit = uint32_max;
164  uint32 requiredFlags = 0u;
166  uint32 inclusionFlags = 0u;
168  uint32 exclusionFlags = 0u;
170  uint32 excessCount = 0u;
172  uint32 rejectCount = 0u;
175  double timeTotal = 0;
176  double timeIO = 0;
177  uint32 countIO = 0u;
178  double timeUncompress = 0;
179  double timeDecode = 0;
180  };
181 
192  template <class StarData>
193  struct SearchData : public SearchDataBase
194  {
196 
201  {
202  stars.Clear();
203  excessCount = rejectCount = 0u;
204  timeTotal = timeIO = 0;
205  countIO = 0u;
206  timeUncompress = timeDecode = 0;
207  }
208  };
209 
210 protected:
211 
212  struct ChildNodeData
213  {
214  // Zero-based quadtree child node positions in an index node array.
215  uint32 nw; // top-left child node
216  uint32 ne; // top-right child node
217  uint32 sw; // bottom-left child node
218  uint32 se; // bottom-right child node
219  };
220 
221 #ifdef _MSC_VER
222  /*
223  * Our favorite brain-damaged thing does not know how to implement bit
224  * fields. Oh well...
225  */
226  struct LeafNodeData
227  {
228  uint64 blockOffsetAndLeafFlag;
229  uint32 blockSize;
230  uint32 compressedBlockSize;
231  };
232 #else
233  struct LeafNodeData
234  {
235  uint64 blockOffset : 63; // position of source data block, byte offset
236  bool leafFlag : 1; // quadtree node type: 0=structural 1=leaf
237  uint32 blockSize; // size of point source data, in bytes
238  uint32 compressedBlockSize; // size of compressed data, in bytes
239  };
240 #endif
241 
242  /*
243  * Quadtree index node (48 bytes).
244  */
245  struct IndexNode
246  {
247  // Projected coordinates of quadtree node rectangle.
248  double x0; // left
249  double y0; // top
250  double x1; // right
251  double y1; // bottom
252 
253  // Quadtree child node indexes or leaf node data.
254  union { ChildNodeData child;
255  LeafNodeData leaf; } index;
256 
257  IndexNode()
258  {
259  static_assert( sizeof( *this ) == 48, "Invalid sizeof( XPSD::IndexNode )" );
260  static_assert( sizeof( ChildNodeData ) == 16, "Invalid sizeof( XPSD::ChildNodeData )" );
261  static_assert( sizeof( LeafNodeData ) == 16, "Invalid sizeof( XPSD::LeafNodeData )" );
262  static_assert( sizeof( index ) == 16, "Invalid sizeof( XPSD::IndexNode::index )" );
263 
264  index.child.nw = index.child.ne = index.child.sw = index.child.se = 0;
265  }
266 
267  bool IsLeaf() const
268  {
269 #ifdef _MSC_VER
270  return (index.leaf.blockOffsetAndLeafFlag & 0x8000000000000000) != 0;
271 #else
272  return index.leaf.leafFlag;
273 #endif
274  }
275 
276  uint64 BlockOffset() const
277  {
278 #ifdef _MSC_VER
279  return index.leaf.blockOffsetAndLeafFlag & 0x7FFFFFFFFFFFFFFF;
280 #else
281  return index.leaf.blockOffset;
282 #endif
283  }
284 
285  uint32 BlockSize() const
286  {
287  return index.leaf.blockSize;
288  }
289 
290  uint32 CompressedBlockSize() const
291  {
292  return index.leaf.compressedBlockSize;
293  }
294  };
295 
296  static double Distance( double lon1, double lat1, double lon2, double lat2 )
297  {
298  return Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) ).Angle3D( Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) ) );
299  }
300 
301  static double CrossTrackDistance( double lon, double lat, double lon1, double lat1, double lon2, double lat2 )
302  {
303  if ( lon == lon1 )
304  if ( lat == lat1 )
305  return 0;
306 
307  Vector p = Vector::FromSpherical( Rad( lon ), Rad( lat ) );
308  Vector c = Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) ).Cross( Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) ) );
309  return c.Angle3D( p ) - Pi()/2;
310  }
311 
312  static bool WithinExtent( double lon, double lat, double lon1, double lat1, double lon2, double lat2 )
313  {
314  if ( lon1 == lon2 )
315  if ( lat1 == lat2 )
316  return lon == lon1 && lat == lat1; // null segment
317 
318  Vector n0 = Vector::FromSpherical( Rad( lon ), Rad( lat ) );
319  Vector n1 = Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) );
320  Vector n2 = Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) );
321 
322  // Get vectors representing p0->p1, p0->p2, p1->p2, p2->p1
323  Vector d10 = n0 - n1, d12 = n2 - n1;
324  Vector d20 = n0 - n2, d21 = n1 - n2;
325 
326  // Dot product d10*d12 tells us if p0 is on p2 side of p1, similarly for d20*d21
327  if ( d10 * d12 >= 0 )
328  if ( d20 * d21 >= 0 )
329  return (n0 * n1) >= 0 && (n0 * n2) >= 0; // same hemisphere
330 
331  return false;
332  }
333 
334  static bool InRegion( double lon, double lat,
335  double lon1, double lat1, double lon2, double lat2,
336  double lon3, double lat3, double lon4, double lat4 )
337  {
338  Vector p = Vector::FromSpherical( Rad( lon ), Rad( lat ) );
339  Vector v1 = p - Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) );
340  Vector v2 = p - Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) );
341  Vector v3 = p - Vector::FromSpherical( Rad( lon3 ), Rad( lat3 ) );
342  Vector v4 = p - Vector::FromSpherical( Rad( lon4 ), Rad( lat4 ) );
343  return Abs( v1.Angle3D( v2, p ) + v2.Angle3D( v3, p ) + v3.Angle3D( v4, p ) + v4.Angle3D( v1, p ) ) > Pi();
344  }
345 
346  enum projection_type { Equirectangular, TransverseEquirectangular, AzimuthalEquidistant };
347 
348  static String ProjectionToAttributeValue( int );
349  static projection_type ProjectionFromAttributeValue( const String& );
350 
351  class IndexTree
352  {
353  public:
354 
355  IndexTree( StarDatabaseFile* parent,
356  projection_type projection, double centerRA, double centerDec,
357  const Array<IndexNode>& nodes )
358  : m_parent( parent )
359  , m_projection( projection )
360  , m_centerRA( centerRA )
361  , m_centerDec( centerDec )
362  , m_nodes( nodes )
363  {
364  }
365 
366  IndexTree() = default;
367  IndexTree( const IndexTree& ) = default;
368 
369  void Project( double& x, double& y, double ra, double dec ) const
370  {
371  switch( m_projection )
372  {
373  default: // ?!
374  case Equirectangular:
375  x = ra - m_centerRA;
376  y = dec;
377  break;
378  case AzimuthalEquidistant:
379  {
380  double sa, ca;
381  SinCos( Rad( ra ), sa, ca );
382  double r = 90 - Abs( dec );
383  x = r*sa;
384  y = r*ca;
385  }
386  break;
387  case TransverseEquirectangular:
388  {
389  double sa, ca;
390  SinCos( Rad( ra ), sa, ca );
391  double sd, cd;
392  SinCos( Rad( Abs( dec ) ), sd, cd );
393  x = Deg( ArcSin( cd*sa ) );
394  y = Deg( ArcTan( sd, cd*ca ) ) - 90;
395  }
396  break;
397  }
398  }
399 
400  void Unproject( double& ra, double& dec, double x, double y ) const
401  {
402  switch( m_projection )
403  {
404  default: // ?!
405  case Equirectangular:
406  ra = x + m_centerRA;
407  dec = y;
408  break;
409  case AzimuthalEquidistant:
410  x = Rad( x );
411  y = Rad( y );
412  ra = Deg( ArcTan( x, y ) );
413  if ( ra < 0 )
414  ra += 360;
415  dec = Deg( ArcSin( Cos( Sqrt( x*x + y*y ) ) ) );
416  if ( m_centerDec < 0 )
417  dec = -dec;
418  break;
419  case TransverseEquirectangular:
420  {
421  double sx, cx;
422  SinCos( Rad( x ), sx, cx );
423  double sy, cy;
424  SinCos( Rad( y + 90 ), sy, cy );
425  ra = Deg( ArcTan( sx, cx*cy ) );
426  if ( ra < 0 )
427  ra += 360;
428  dec = Deg( ArcSin( sy*cx ) );
429  if ( m_centerDec < 0 )
430  dec = -dec;
431  }
432  break;
433  }
434  }
435 
436  void Search( double ra, double dec, double r, void* searchData ) const
437  {
438  SearchRecursive( 0, ra, dec, r, searchData );
439  }
440 
441  private:
442 
443  StarDatabaseFile* m_parent = nullptr;
444  projection_type m_projection = Equirectangular;
445  double m_centerRA = 0;
446  double m_centerDec = 0;
447  Array<IndexNode> m_nodes;
448 
449  void GetNodeBounds( double& ra1, double& dec1, double& ra2, double& dec2,
450  double& ra3, double& dec3, double& ra4, double& dec4, const IndexNode& node ) const
451  {
452  Unproject( ra1, dec1, node.x0, node.y0 );
453  Unproject( ra2, dec2, node.x1, node.y0 );
454  Unproject( ra3, dec3, node.x1, node.y1 );
455  Unproject( ra4, dec4, node.x0, node.y1 );
456  }
457 
458  bool InNodeRegion( double ra, double dec, const IndexNode& node ) const
459  {
460  double ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4;
461  GetNodeBounds( ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4, node );
462  return InRegion( ra, dec, ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4 );
463  }
464 
465  bool IntersectsNodeRegion( double ra, double dec, double r, const IndexNode& node ) const
466  {
467  double ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4;
468  GetNodeBounds( ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4, node );
469  double rr = Rad( r );
470  return InRegion( ra, dec, ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4 )
471  || Distance( ra, dec, ra1, dec1 ) < rr
472  || Distance( ra, dec, ra2, dec2 ) < rr
473  || Distance( ra, dec, ra3, dec3 ) < rr
474  || Distance( ra, dec, ra4, dec4 ) < rr
475  || WithinExtent( ra, dec, ra1, dec1, ra2, dec2 ) && CrossTrackDistance( ra, dec, ra1, dec1, ra2, dec2 ) < rr
476  || WithinExtent( ra, dec, ra2, dec2, ra3, dec3 ) && CrossTrackDistance( ra, dec, ra2, dec2, ra3, dec3 ) < rr
477  || WithinExtent( ra, dec, ra3, dec3, ra4, dec4 ) && CrossTrackDistance( ra, dec, ra3, dec3, ra4, dec4 ) < rr
478  || WithinExtent( ra, dec, ra4, dec4, ra1, dec1 ) && CrossTrackDistance( ra, dec, ra4, dec4, ra1, dec1 ) < rr;
479  }
480 
481  // Defined after StarDatabaseFile declaration.
482  void SearchRecursive( uint32 nodeIndex, double ra, double dec, double r, void* searchData ) const;
483 
484  friend class StarDatabaseFile;
485  };
486 };
487 
488 // ----------------------------------------------------------------------------
489 
506 class PCL_CLASS StarDatabaseFile : public XPSD
507 {
508 public:
509 
516  StarDatabaseFile() = default;
517 
525  StarDatabaseFile( const String& filePath )
526  {
527  Open( filePath );
528  }
529 
534 
538  StarDatabaseFile& operator =( StarDatabaseFile&& ) = default;
539 
544  StarDatabaseFile( const StarDatabaseFile& ) = delete;
545 
550  StarDatabaseFile& operator =( const StarDatabaseFile& ) = delete;
551 
555  virtual ~StarDatabaseFile() noexcept( false )
556  {
557  }
558 
572  void Open( const String& filePath );
573 
583  void Close();
584 
589  bool IsOpen() const
590  {
591  return m_file.IsOpen();
592  }
593 
598  const String& FilePath() const
599  {
600  return m_file.FilePath();
601  }
602 
608  float MagnitudeLow() const
609  {
610  return m_magnitudeLow;
611  }
612 
618  float MagnitudeHigh() const
619  {
620  return m_magnitudeHigh;
621  }
622 
627  const XPSD::Metadata& Metadata() const
628  {
629  return m_metadata;
630  }
631 
638  {
639  return m_statistics;
640  }
641 
691  static void Serialize( const String& filePath,
692  const XPSD::Metadata& metadata,
693  const XPSD::Statistics& statistics,
694  float magnitudeLow, float magnitudeHigh,
695  const Array<XPSD::IndexTree>& index,
696  const ByteArray& data,
697  const Compression* compression = nullptr,
698  const String& parameters = String() );
699 
700 protected:
701 
702  mutable File m_file;
703  XPSD::Metadata m_metadata;
704  XPSD::Statistics m_statistics;
705  float m_magnitudeLow = 0;
706  float m_magnitudeHigh = 0;
707  Array<XPSD::IndexTree> m_index;
708  uint64 m_dataPosition = 0;
709  AutoPointer<Compression> m_compression;
710  String m_parameters;
711 
712  virtual void LoadData( void* block, uint64 offset, uint32 size, void* searchData ) const
713  {
714  ElapsedTime T;
715  m_file.SetPosition( m_dataPosition + offset );
716  m_file.Read( block, size );
717  reinterpret_cast<SearchDataBase*>( searchData )->timeIO += T();
718  ++reinterpret_cast<SearchDataBase*>( searchData )->countIO;
719  }
720 
721  virtual void Uncompress( ByteArray& block, uint32 uncompressedSize, void* searchData ) const
722  {
723  if ( m_compression )
724  {
725  ElapsedTime T;
726  block = m_compression->Uncompress( block, uncompressedSize );
727  reinterpret_cast<SearchDataBase*>( searchData )->timeUncompress += T();
728  }
729  }
730 
731  virtual void GetEncodedData( const ByteArray&, const XPSD::IndexTree&, const XPSD::IndexNode&, void* ) const = 0;
732 
733  friend class XPSD::IndexTree;
734 };
735 
736 // ----------------------------------------------------------------------------
737 
738 inline void
739 XPSD::IndexTree::SearchRecursive( uint32 nodeIndex, double ra, double dec, double r, void* searchData ) const
740 {
741  const IndexNode& node = m_nodes[nodeIndex];
742  if ( IntersectsNodeRegion( ra, dec, r, node ) )
743  {
744  if ( node.IsLeaf() )
745  {
746  ByteArray block( size_type( node.CompressedBlockSize() ) );
747  m_parent->LoadData( block.Begin(), node.BlockOffset(), node.CompressedBlockSize(), searchData );
748  if ( node.CompressedBlockSize() != node.BlockSize() )
749  m_parent->Uncompress( block, node.BlockSize(), searchData );
750  m_parent->GetEncodedData( block, *this, node, searchData );
751  }
752  else
753  {
754  if ( node.index.child.nw != 0 )
755  SearchRecursive( node.index.child.nw, ra, dec, r, searchData );
756  if ( node.index.child.ne != 0 )
757  SearchRecursive( node.index.child.ne, ra, dec, r, searchData );
758  if ( node.index.child.sw != 0 )
759  SearchRecursive( node.index.child.sw, ra, dec, r, searchData );
760  if ( node.index.child.se != 0 )
761  SearchRecursive( node.index.child.se, ra, dec, r, searchData );
762  }
763  }
764 }
765 
766 // ----------------------------------------------------------------------------
767 
768 } // pcl
769 
770 #endif // __PCL_StarDatabaseFile_h
771 
772 // ----------------------------------------------------------------------------
773 // EOF pcl/StarDatabaseFile.h - Released 2024-06-18T15:48:54Z
iterator Begin()
Definition: Array.h:426
void Clear()
Definition: Array.h:1173
A smart pointer with exclusive object ownership and optional automatic object destruction.
Definition: AutoPointer.h:241
Dynamic array of 8-bit unsigned integers.
Abstract base class of data compression algorithm implementations.
Definition: Compression.h:84
High-resolution time stamp.
Definition: ElapsedTime.h:131
A platform-independent interface to the local file system.
Definition: File.h:344
virtual void Read(void *buffer, fsize_type len)
virtual void SetPosition(fpos_type pos)
static GenericVector FromSpherical(double slon, double clon, double slat, double clat)
Definition: Vector.h:2143
double Angle3D(const GenericVector &v) const noexcept
Definition: Vector.h:2221
GenericVector Cross(const GenericVector &v2) const
Definition: Vector.h:665
Point source and star catalog database files (XPSD format).
virtual ~StarDatabaseFile() noexcept(false)
const XPSD::Metadata & Metadata() const
StarDatabaseFile()=default
static void Serialize(const String &filePath, const XPSD::Metadata &metadata, const XPSD::Statistics &statistics, float magnitudeLow, float magnitudeHigh, const Array< XPSD::IndexTree > &index, const ByteArray &data, const Compression *compression=nullptr, const String &parameters=String())
StarDatabaseFile(StarDatabaseFile &&)=default
const XPSD::Statistics & Statistics() const
StarDatabaseFile(const StarDatabaseFile &)=delete
void Open(const String &filePath)
float MagnitudeLow() const
const String & FilePath() const
float MagnitudeHigh() const
StarDatabaseFile(const String &filePath)
Unicode (UTF-16) string.
Definition: String.h:8113
An instant in any timescale.
Definition: TimePoint.h:103
64-bit floating point real vector.
Base class of point source database implementations.
Complex< T > Sqrt(const Complex< T > &c) noexcept
Definition: Complex.h:674
T Abs(const Complex< T > &c) noexcept
Definition: Complex.h:429
Complex< T > Cos(const Complex< T > &c) noexcept
Definition: Complex.h:806
#define uint32_max
Definition: Defs.h:875
void SinCos(T x, T &sx, T &cx) noexcept
Definition: Math.h:1030
constexpr T ArcTan(T x) noexcept
Definition: Math.h:526
constexpr T ArcSin(T x) noexcept
Definition: Math.h:514
constexpr T Rad(T x) noexcept
Definition: Math.h:1894
constexpr T Deg(T x) noexcept
Definition: Math.h:606
constexpr long double Pi() noexcept
Definition: Math.h:461
unsigned long long uint64
Definition: Defs.h:682
unsigned int uint32
Definition: Defs.h:666
FI1 Search(FI1 i1, FI1 j1, FI2 i2, FI2 j2) noexcept
Definition: Search.h:397
size_t size_type
Definition: Defs.h:609
PCL root namespace.
Definition: AbstractImage.h:77
distance_type Distance(FI i, FI j)
Definition: Iterator.h:161
Metadata items available in point source database files.
String briefDescription
A brief (single-line) description of this XPSD file.
String title
A title that represents or identifies this XPSD file.
String copyright
Copyright information applicable to the data stored in this XPSD file.
String organizationName
The name of the organization responsible for this file.
String authors
The names of one or more persons or groups that have created the data in this file.
String databaseIdentifier
The unique identifier of the database this file belongs to.
String creatorApplication
The software application or program that created this file.
String description
A full description of the data stored in this XPSD file.
String creatorOS
The operating system on which this file was created.
String databaseVersion
The version of the database this file belongs to.
TimePoint creationTime
The date this file was created.
Parameters and output instrumentation data for catalog search operations.
Data items and parameters for catalog search operations.
Array< StarData > stars
The list of stars found by the search operation (output data).
Structural and statistical data about an XPSD database file.