StarDatabaseFile.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.6.5
// ----------------------------------------------------------------------------
// pcl/StarDatabaseFile.h - Released 2024-01-13T15:47:58Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_StarDatabaseFile_h
#define __PCL_StarDatabaseFile_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/AutoPointer.h>
#include <pcl/Compression.h>
#include <pcl/ElapsedTime.h>
#include <pcl/File.h>
#include <pcl/TimePoint.h>
#include <pcl/Vector.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS StarDatabaseFile;
 
class PCL_CLASS XPSD
{
public:
 
   struct Metadata
   {
      TimePoint creationTime;       
      String    creatorOS;          
      String    creatorApplication; 
      String    databaseIdentifier; 
      String    databaseVersion;    
      String    title;              
      String    briefDescription;   
      String    description;        
      String    organizationName;   
      String    authors;            
      String    copyright;          
   };
 
   struct Statistics
   {
      uint64   totalSources = 0;      
      uint32   totalNodes = 0;        
      uint32   totalLeaves = 0;       
      float    medianLeafLength = 0;  
      uint32   minimumLeafLength = 0; 
      uint32   maximumLeafLength = 0; 
   };
 
   struct SearchDataBase
   {
      double   centerRA = 0;             
      double   centerDec = 0;            
      double   radius = 1;               
      float    magnitudeLow = -1.5;      
      float    magnitudeHigh = 26;       
      uint32   sourceLimit = uint32_max; 
      uint32   requiredFlags = 0u;       
      uint32   inclusionFlags = 0u;      
      uint32   exclusionFlags = 0u;      
      uint32   excessCount = 0u;         
      uint32   rejectCount = 0u;         
      double   timeTotal = 0;            
      double   timeIO = 0;               
      uint32   countIO = 0u;             
      double   timeUncompress = 0;       
      double   timeDecode = 0;           
   };
 
   template <class StarData>
   struct SearchData : public SearchDataBase
   {
      Array<StarData> stars; 
 
      void ResetSearchResults()
      {
         stars.Clear();
         excessCount = rejectCount = 0u;
         timeTotal = timeIO = 0;
         countIO = 0u;
         timeUncompress = timeDecode = 0;
      }
   };
 
protected:
 
   struct ChildNodeData
   {
      // Zero-based quadtree child node positions in an index node array.
      uint32 nw; // top-left child node
      uint32 ne; // top-right child node
      uint32 sw; // bottom-left child node
      uint32 se; // bottom-right child node
   };
 
#ifdef _MSC_VER
   /*
    * Our favorite brain-damaged thing does not know how to implement bit
    * fields. Oh well...
    */
   struct LeafNodeData
   {
      uint64 blockOffsetAndLeafFlag;
      uint32 blockSize;
      uint32 compressedBlockSize;
   };
#else
   struct LeafNodeData
   {
      uint64 blockOffset : 63;    // position of source data block, byte offset
      bool   leafFlag    :  1;    // quadtree node type: 0=structural 1=leaf
      uint32 blockSize;           // size of point source data, in bytes
      uint32 compressedBlockSize; // size of compressed data, in bytes
   };
#endif
 
   /*
    * Quadtree index node (48 bytes).
    */
   struct IndexNode
   {
      // Projected coordinates of quadtree node rectangle.
      double x0; // left
      double y0; // top
      double x1; // right
      double y1; // bottom
 
      // Quadtree child node indexes or leaf node data.
      union { ChildNodeData child;
              LeafNodeData  leaf; } index;
 
      IndexNode()
      {
         static_assert( sizeof( *this ) == 48, "Invalid sizeof( XPSD::IndexNode )" );
         static_assert( sizeof( ChildNodeData ) == 16, "Invalid sizeof( XPSD::ChildNodeData )" );
         static_assert( sizeof( LeafNodeData ) == 16, "Invalid sizeof( XPSD::LeafNodeData )" );
         static_assert( sizeof( index ) == 16, "Invalid sizeof( XPSD::IndexNode::index )" );
 
         index.child.nw = index.child.ne = index.child.sw = index.child.se = 0;
      }
 
      bool IsLeaf() const
      {
#ifdef _MSC_VER
         return (index.leaf.blockOffsetAndLeafFlag & 0x8000000000000000) != 0;
#else
         return index.leaf.leafFlag;
#endif
      }
 
      uint64 BlockOffset() const
      {
#ifdef _MSC_VER
         return index.leaf.blockOffsetAndLeafFlag & 0x7FFFFFFFFFFFFFFF;
#else
         return index.leaf.blockOffset;
#endif
      }
 
      uint32 BlockSize() const
      {
         return index.leaf.blockSize;
      }
 
      uint32 CompressedBlockSize() const
      {
         return index.leaf.compressedBlockSize;
      }
   };
 
   static double Distance( double lon1, double lat1, double lon2, double lat2 )
   {
      return Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) ).Angle3D( Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) ) );
   }
 
   static double CrossTrackDistance( double lon, double lat, double lon1, double lat1, double lon2, double lat2 )
   {
      if ( lon == lon1 )
         if ( lat == lat1 )
            return 0;
 
      Vector p = Vector::FromSpherical( Rad( lon ), Rad( lat ) );
      Vector c = Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) ).Cross( Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) ) );
      return c.Angle3D( p ) - Pi()/2;
   }
 
   static bool WithinExtent( double lon, double lat, double lon1, double lat1, double lon2, double lat2 )
   {
      if ( lon1 == lon2 )
         if ( lat1 == lat2 )
            return lon == lon1 && lat == lat1; // null segment
 
      Vector n0 = Vector::FromSpherical( Rad( lon ), Rad( lat ) );
      Vector n1 = Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) );
      Vector n2 = Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) );
 
      // Get vectors representing p0->p1, p0->p2, p1->p2, p2->p1
      Vector d10 = n0 - n1, d12 = n2 - n1;
      Vector d20 = n0 - n2, d21 = n1 - n2;
 
      // Dot product d10*d12 tells us if p0 is on p2 side of p1, similarly for d20*d21
      if ( d10 * d12 >= 0 )
         if ( d20 * d21 >= 0 )
            return (n0 * n1) >= 0 && (n0 * n2) >= 0; // same hemisphere
 
      return false;
   }
 
   static bool InRegion( double lon, double lat,
                         double lon1, double lat1, double lon2, double lat2,
                         double lon3, double lat3, double lon4, double lat4 )
   {
      Vector p = Vector::FromSpherical( Rad( lon ), Rad( lat ) );
      Vector v1 = p - Vector::FromSpherical( Rad( lon1 ), Rad( lat1 ) );
      Vector v2 = p - Vector::FromSpherical( Rad( lon2 ), Rad( lat2 ) );
      Vector v3 = p - Vector::FromSpherical( Rad( lon3 ), Rad( lat3 ) );
      Vector v4 = p - Vector::FromSpherical( Rad( lon4 ), Rad( lat4 ) );
      return Abs( v1.Angle3D( v2, p ) + v2.Angle3D( v3, p ) + v3.Angle3D( v4, p ) + v4.Angle3D( v1, p ) ) > Pi();
   }
 
   enum projection_type { Equirectangular, TransverseEquirectangular, AzimuthalEquidistant };
 
   static String ProjectionToAttributeValue( int );
   static projection_type ProjectionFromAttributeValue( const String& );
 
   class IndexTree
   {
   public:
 
      IndexTree( StarDatabaseFile* parent,
                 projection_type projection, double centerRA, double centerDec,
                 const Array<IndexNode>& nodes )
         : m_parent( parent )
         , m_projection( projection )
         , m_centerRA( centerRA )
         , m_centerDec( centerDec )
         , m_nodes( nodes )
      {
      }
 
      IndexTree() = default;
      IndexTree( const IndexTree& ) = default;
 
      void Project( double& x, double& y, double ra, double dec ) const
      {
         switch( m_projection )
         {
         default: // ?!
         case Equirectangular:
            x = ra - m_centerRA;
            y = dec;
            break;
         case AzimuthalEquidistant:
            {
               double sa, ca;
               SinCos( Rad( ra ), sa, ca );
               double r = 90 - Abs( dec );
               x = r*sa;
               y = r*ca;
            }
            break;
         case TransverseEquirectangular:
            {
               double sa, ca;
               SinCos( Rad( ra ), sa, ca );
               double sd, cd;
               SinCos( Rad( Abs( dec ) ), sd, cd );
               x = Deg( ArcSin( cd*sa ) );
               y = Deg( ArcTan( sd, cd*ca ) ) - 90;
            }
            break;
         }
      }
 
      void Unproject( double& ra, double& dec, double x, double y ) const
      {
         switch( m_projection )
         {
         default: // ?!
         case Equirectangular:
            ra = x + m_centerRA;
            dec = y;
            break;
         case AzimuthalEquidistant:
            x = Rad( x );
            y = Rad( y );
            ra = Deg( ArcTan( x, y ) );
            if ( ra < 0 )
               ra += 360;
            dec = Deg( ArcSin( Cos( Sqrt( x*x + y*y ) ) ) );
            if ( m_centerDec < 0 )
               dec = -dec;
            break;
         case TransverseEquirectangular:
            {
               double sx, cx;
               SinCos( Rad( x ), sx, cx );
               double sy, cy;
               SinCos( Rad( y + 90 ), sy, cy );
               ra = Deg( ArcTan( sx, cx*cy ) );
               if ( ra < 0 )
                  ra += 360;
               dec = Deg( ArcSin( sy*cx ) );
               if ( m_centerDec < 0 )
                  dec = -dec;
            }
            break;
         }
      }
 
      void Search( double ra, double dec, double r, void* searchData ) const
      {
         SearchRecursive( 0, ra, dec, r, searchData );
      }
 
   private:
 
      StarDatabaseFile* m_parent = nullptr;
      projection_type   m_projection = Equirectangular;
      double            m_centerRA = 0;
      double            m_centerDec = 0;
      Array<IndexNode>  m_nodes;
 
      void GetNodeBounds( double& ra1, double& dec1, double& ra2, double& dec2,
                          double& ra3, double& dec3, double& ra4, double& dec4, const IndexNode& node ) const
      {
         Unproject( ra1, dec1, node.x0, node.y0 );
         Unproject( ra2, dec2, node.x1, node.y0 );
         Unproject( ra3, dec3, node.x1, node.y1 );
         Unproject( ra4, dec4, node.x0, node.y1 );
      }
 
      bool InNodeRegion( double ra, double dec, const IndexNode& node ) const
      {
         double ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4;
         GetNodeBounds( ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4, node );
         return InRegion( ra, dec, ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4 );
      }
 
      bool IntersectsNodeRegion( double ra, double dec, double r, const IndexNode& node ) const
      {
         double ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4;
         GetNodeBounds( ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4, node );
         double rr = Rad( r );
         return InRegion( ra, dec, ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4 )
             || Distance( ra, dec, ra1, dec1 ) < rr
             || Distance( ra, dec, ra2, dec2 ) < rr
             || Distance( ra, dec, ra3, dec3 ) < rr
             || Distance( ra, dec, ra4, dec4 ) < rr
             || WithinExtent( ra, dec, ra1, dec1, ra2, dec2 ) && CrossTrackDistance( ra, dec, ra1, dec1, ra2, dec2 ) < rr
             || WithinExtent( ra, dec, ra2, dec2, ra3, dec3 ) && CrossTrackDistance( ra, dec, ra2, dec2, ra3, dec3 ) < rr
             || WithinExtent( ra, dec, ra3, dec3, ra4, dec4 ) && CrossTrackDistance( ra, dec, ra3, dec3, ra4, dec4 ) < rr
             || WithinExtent( ra, dec, ra4, dec4, ra1, dec1 ) && CrossTrackDistance( ra, dec, ra4, dec4, ra1, dec1 ) < rr;
      }
 
      // Defined after StarDatabaseFile declaration.
      void SearchRecursive( uint32 nodeIndex, double ra, double dec, double r, void* searchData ) const;
 
      friend class StarDatabaseFile;
   };
};
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS StarDatabaseFile : public XPSD
{
public:
 
   StarDatabaseFile() = default;
 
   StarDatabaseFile( const String& filePath )
   {
      Open( filePath );
   }
 
   StarDatabaseFile( StarDatabaseFile&& ) = default;
 
   StarDatabaseFile& operator =( StarDatabaseFile&& ) = default;
 
   StarDatabaseFile( const StarDatabaseFile& ) = delete;
 
   StarDatabaseFile& operator =( const StarDatabaseFile& ) = delete;
 
   virtual ~StarDatabaseFile() noexcept( false )
   {
   }
 
   void Open( const String& filePath );
 
   void Close();
 
   bool IsOpen() const
   {
      return m_file.IsOpen();
   }
 
   const String& FilePath() const
   {
      return m_file.FilePath();
   }
 
   float MagnitudeLow() const
   {
      return m_magnitudeLow;
   }
 
   float MagnitudeHigh() const
   {
      return m_magnitudeHigh;
   }
 
   const XPSD::Metadata& Metadata() const
   {
      return m_metadata;
   }
 
   const XPSD::Statistics& Statistics() const
   {
      return m_statistics;
   }
 
   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() );
 
protected:
 
   mutable File                     m_file;
           XPSD::Metadata           m_metadata;
           XPSD::Statistics         m_statistics;
           float                    m_magnitudeLow = 0;
           float                    m_magnitudeHigh = 0;
           Array<XPSD::IndexTree>   m_index;
           uint64                   m_dataPosition = 0;
           AutoPointer<Compression> m_compression;
           String                   m_parameters;
 
   virtual void LoadData( void* block, uint64 offset, uint32 size, void* searchData ) const
   {
      ElapsedTime T;
      m_file.SetPosition( m_dataPosition + offset );
      m_file.Read( block, size );
      reinterpret_cast<SearchDataBase*>( searchData )->timeIO += T();
      ++reinterpret_cast<SearchDataBase*>( searchData )->countIO;
   }
 
   virtual void Uncompress( ByteArray& block, uint32 uncompressedSize, void* searchData ) const
   {
      if ( m_compression )
      {
         ElapsedTime T;
         block = m_compression->Uncompress( block, uncompressedSize );
         reinterpret_cast<SearchDataBase*>( searchData )->timeUncompress += T();
      }
   }
 
   virtual void GetEncodedData( const ByteArray&, const XPSD::IndexTree&, const XPSD::IndexNode&, void* ) const = 0;
 
   friend class XPSD::IndexTree;
};
 
// ----------------------------------------------------------------------------
 
inline void
XPSD::IndexTree::SearchRecursive( uint32 nodeIndex, double ra, double dec, double r, void* searchData ) const
{
   const IndexNode& node = m_nodes[nodeIndex];
   if ( IntersectsNodeRegion( ra, dec, r, node ) )
   {
      if ( node.IsLeaf() )
      {
         ByteArray block( size_type( node.CompressedBlockSize() ) );
         m_parent->LoadData( block.Begin(), node.BlockOffset(), node.CompressedBlockSize(), searchData );
         if ( node.CompressedBlockSize() != node.BlockSize() )
            m_parent->Uncompress( block, node.BlockSize(), searchData );
         m_parent->GetEncodedData( block, *this, node, searchData );
      }
      else
      {
         if ( node.index.child.nw != 0 )
            SearchRecursive( node.index.child.nw, ra, dec, r, searchData );
         if ( node.index.child.ne != 0 )
            SearchRecursive( node.index.child.ne, ra, dec, r, searchData );
         if ( node.index.child.sw != 0 )
            SearchRecursive( node.index.child.sw, ra, dec, r, searchData );
         if ( node.index.child.se != 0 )
            SearchRecursive( node.index.child.se, ra, dec, r, searchData );
      }
   }
}
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif  // __PCL_StarDatabaseFile_h
 
// ----------------------------------------------------------------------------
// EOF pcl/StarDatabaseFile.h - Released 2024-01-13T15:47:58Z