GaiaDatabaseFile.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/GaiaDatabaseFile.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// 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_GaiaDatabaseFile_h
#define __PCL_GaiaDatabaseFile_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/Flags.h>
#include <pcl/MetaParameter.h> // pcl_bool
#include <pcl/StarDatabaseFile.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
namespace GaiaStarFlag
{
   enum mask_type
   {
      NoPM             = 0x00000001,
      NoGBPMag         = 0x00000002,
      NoGRPMag         = 0x00000004,
      LackingData      = 0x00000007,
 
      GoldRA           = 0x00000010,
      GoldDec          = 0x00000020,
      GoldPMRA         = 0x00000040,
      GoldPMDec        = 0x00000080,
 
      SilverRA         = 0x00000100,
      SilverDec        = 0x00000200,
      SilverPMRA       = 0x00000400,
      SilverPMDec      = 0x00000800,
 
      BronzeRA         = 0x00001000,
      BronzeDec        = 0x00002000,
      BronzePMRA       = 0x00004000,
      BronzePMDec      = 0x00008000,
 
      GoldGMag         = 0x00010000,
      GoldGBPMag       = 0x00020000,
      GoldGRPMag       = 0x00040000,
      GoldParx         = 0x00080000,
 
      SilverGMag       = 0x00100000,
      SilverGBPMag     = 0x00200000,
      SilverGRPMag     = 0x00400000,
      SilverParx       = 0x00800000,
 
      BronzeGMag       = 0x01000000,
      BronzeGBPMag     = 0x02000000,
      BronzeGRPMag     = 0x04000000,
      BronzeParx       = 0x08000000,
 
      BPRPExcess       = 0x10000008,
      BPRPExcessHigh   = 0x20000000,
 
      GoldAstrometry   = 0x000800F0,
      SilverAstrometry = 0x00800F00,
      BronzeAstrometry = 0x0800F000,
 
      GoldPhotometry   = 0x00070000,
      SilverPhotometry = 0x00700000,
      BronzePhotometry = 0x07000000
   };
}
 
// ----------------------------------------------------------------------------
 
struct PCL_CLASS GaiaStarDataBase
{
   double  ra = 0;     
   double  dec = 0;    
   float   parx = 0;   
   float   pmra = 0;   
   float   pmdec = 0;  
   float   magG = 0;   
   float   magBP = 0;  
   float   magRP = 0;  
   uint32  flags = 0u; 
};
 
struct PCL_CLASS GaiaStarData : public GaiaStarDataBase
{
   FVector flux; 
};
 
// ----------------------------------------------------------------------------
 
struct GaiaSearchData : public XPSD::SearchData<GaiaStarData>
{
   pcl_bool normalizeSpectrum = false; 
   pcl_bool photonFluxUnits = false; 
};
 
// ----------------------------------------------------------------------------
 
class GaiaDatabaseFile : public StarDatabaseFile
{
public:
 
   GaiaDatabaseFile() = default;
 
   GaiaDatabaseFile( const String& filePath )
      : StarDatabaseFile( filePath )
   {
      static_assert( sizeof( EncodedStarData ) == 32, "Invalid sizeof( GaiaDatabaseFile::EncodedStarData )" );
      static_assert( sizeof( EncodedStarSPData ) == 41, "Invalid sizeof( GaiaDatabaseFile::EncodedStarSPData )" );
      if ( Metadata().databaseIdentifier == "GaiaDR3" )
         m_dr = "DR3";
      else if ( Metadata().databaseIdentifier == "GaiaDR3SP" )
      {
         m_dr = "DR3";
         m_hasSpectrumData = true;
         // spectrumStart=<f>,spectrumStep=<f>,spectrumCount=<n>,spectrumBits=<n>
         if ( m_parameters.IsEmpty() )
            throw Error( "Missing parameters Data element attribute: " + filePath );
         StringList tokens;
         m_parameters.Break( tokens, ',' );
         for ( const String& token : tokens )
         {
            StringList items;
            token.Break( items, '=' );
            if ( items.Length() != 2 )
               throw Error( "Invalid parameter specification '" + token + "': " + filePath );
            if ( items[0] == "spectrumStart" )
               m_spectrumStart = items[1].ToFloat();
            else if ( items[0] == "spectrumStep" )
               m_spectrumStep = items[1].ToFloat();
            else if ( items[0] == "spectrumCount" )
               m_spectrumCount = items[1].ToInt();
            else if ( items[0] == "spectrumBits" )
               m_spectrumBits = items[1].ToInt();
         }
         if ( m_spectrumStart <= 0 || m_spectrumStep <= 0 || m_spectrumCount <= 0 ||
              m_spectrumBits != 8 && m_spectrumBits != 16 )
            throw Error( "Invalid parameters attribute value '" + m_parameters + "': " + filePath );
         m_spectrumRange = (1 << m_spectrumBits) - 1;
      }
      else if ( Metadata().databaseIdentifier == "GaiaEDR3" )
         m_dr = "EDR3";
      else if ( Metadata().databaseIdentifier == "GaiaDR2" )
      {
         m_dr = "DR2";
         StringList tokens;
         Metadata().databaseVersion.Break( tokens, '.' );
         // Make sure we reject an unsupported DR2 version older than 1.0.2
         if ( tokens.Length() < 3 || tokens[0].ToInt() < 1 || tokens[1].ToInt() < 0 || tokens[2].ToInt() < 2 )
            throw Error( "Unsupported Gaia DR2 database version '"
                        + Metadata().databaseVersion  + "': " + filePath );
      }
      else
         throw Error( "Invalid or unsupported Gaia database file with unknown identifier '"
                     + Metadata().databaseIdentifier  + "': " + filePath );
   }
 
   GaiaDatabaseFile& operator =( GaiaDatabaseFile&& ) = default;
 
   GaiaDatabaseFile( const GaiaDatabaseFile& ) = delete;
 
   GaiaDatabaseFile& operator =( const GaiaDatabaseFile& ) = delete;
 
   const IsoString& DataRelease() const
   {
      return m_dr;
   }
 
   bool HasMeanSpectrumData() const
   {
      return m_hasSpectrumData;
   }
 
   float SpectrumStart() const
   {
      return m_spectrumStart;
   }
 
   float SpectrumStep() const
   {
      return m_spectrumStep;
   }
 
   int SpectrumCount() const
   {
      return m_spectrumCount;
   }
 
   int SpectrumBits() const
   {
      return m_spectrumBits;
   }
 
   void Search( GaiaSearchData& data ) const
   {
      ElapsedTime T;
      for ( const XPSD::IndexTree& tree : m_index )
         tree.Search( data.centerRA, data.centerDec, data.radius, &data );
      data.timeTotal += T();
   }
 
private:
 
   IsoString m_dr;            // data release, one of "DR2", "EDR3", "DR3"
   bool      m_hasSpectrumData = false;   // with BP/RP mean spectrum data?
   float     m_spectrumStart = 0;         // nm
   float     m_spectrumStep = 0;          // nm
   int       m_spectrumCount = 0;
   int       m_spectrumBits = 0;          // 8 | 16
   float     m_spectrumRange = 0;         // (1 << m_spectrumBits) - 1
 
#pragma pack(push, 1)
 
   /*
    * Encoded star record (32 bytes uncompressed).
    */
   struct EncodedStarData
   {
      // Projected coordinates relative to the origin of the parent quadtree
      // node, in 0.002 mas units.
      uint32 dx;
      uint32 dy;
      // Parallax in mas units.
      float  parx;
      // Proper motions, mas/yr.
      float  pmra;
      float  pmdec;
      // Mean magnitudes in 0.001 mag units, encoded as (mag + 1.5)*1000.
      uint16 magG;
      uint16 magBP;
      uint16 magRP;
      // Right ascension correction for high declinations, in 0.01 mas units.
      int16  dra;
      // Data availability and quality flags.
      uint32 flags;
   };
 
   /*
    * Encoded star record with mean spectrum data.
    * Uncompressed size in bytes: 40 + (m_spectrumBits/8)*(m_spectrumCount + (m_spectrumCount & 1)).
    */
   struct EncodedStarSPData : public EncodedStarData
   {
      float fluxMin;
      float fluxMul; // (fluxMax - fluxMin)/(2**m_spectrumBits - 1)
      uint8 flux[ 1 ];
   };
 
   constexpr size_type EncodedStarSPDataSize() const
   {
      return sizeof( EncodedStarSPData ) - 1 + (m_spectrumBits >> 3)*(m_spectrumCount + (m_spectrumCount & 1));
   }
 
#pragma pack(pop)
 
   void GetEncodedData( const ByteArray& data, const XPSD::IndexTree& tree, const XPSD::IndexNode& node, void* searchData ) const override
   {
      ElapsedTime T;
      GaiaSearchData* search = reinterpret_cast<GaiaSearchData*>( searchData );
      double searchRadius = Rad( search->radius );
      int itemSize = m_hasSpectrumData ? EncodedStarSPDataSize() : sizeof( EncodedStarData );
      int count = int( data.Size() / itemSize );
      int matchCount = 0;
      for ( int i = 0; i < count; ++i )
      {
         const EncodedStarData* S = reinterpret_cast<const EncodedStarData*>( data.Begin() + i*itemSize );
         if ( search->requiredFlags == 0 || (S->flags & search->requiredFlags) == search->requiredFlags )
            if ( search->inclusionFlags == 0 || (S->flags & search->inclusionFlags) != 0 )
               if ( search->exclusionFlags == 0 || (S->flags & search->exclusionFlags) == 0 )
               {
                  float magG = 0.001*S->magG - 1.5;
                  if ( magG >= search->magnitudeLow )
                     if ( magG <= search->magnitudeHigh )
                     {
                        double x = node.x0 + double( S->dx )/3600/1000/500;
                        double y = node.y0 + double( S->dy )/3600/1000/500;
                        GaiaStarData star;
                        tree.Unproject( star.ra, star.dec, x, y );
                        if ( unlikely( S->dra != 0 ) )
                        {
                           star.ra += double( S->dra )/3600/1000/100;
                           if ( star.ra < 0 )
                              star.ra += 360;
                           else if ( star.ra >= 360 )
                              star.ra -= 360;
                        }
                        if ( Distance( search->centerRA, search->centerDec, star.ra, star.dec ) < searchRadius )
                        {
                           ++matchCount;
 
                           if ( search->stars.Length() < size_type( search->sourceLimit ) )
                           {
                              star.parx = S->parx;
                              star.pmra = S->pmra;
                              star.pmdec = S->pmdec;
                              star.magG = magG;
                              star.magBP = 0.001*S->magBP - 1.5;
                              star.magRP = 0.001*S->magRP - 1.5;
                              star.flags = S->flags;
 
                              if ( m_hasSpectrumData )
                              {
                                 const EncodedStarSPData* SS = static_cast<const EncodedStarSPData*>( S );
                                 star.flux = FVector( m_spectrumCount );
                                 if ( search->photonFluxUnits )
                                 {
                                    for ( int j = 0; j < m_spectrumCount; ++j )
                                       star.flux[j] = (SS->flux[j]*SS->fluxMul + SS->fluxMin) * (m_spectrumStart + j*m_spectrumStep)/1.602e-19/1.239979e-3;
                                    if ( search->normalizeSpectrum )
                                       star.flux /= star.flux.MaxComponent();
                                 }
                                 else
                                 {
                                    if ( search->normalizeSpectrum )
                                       for ( int j = 0; j < m_spectrumCount; ++j )
                                          star.flux[j] = SS->flux[j]/m_spectrumRange;
                                    else
                                       for ( int j = 0; j < m_spectrumCount; ++j )
                                          star.flux[j] = SS->flux[j]*SS->fluxMul + SS->fluxMin;
                                 }
                              }
 
                              search->stars << star;
                           }
                           else
                              ++search->excessCount;
                        }
                     }
               }
      }
 
      search->rejectCount += count - matchCount;
      search->timeDecode += T();
   }
 
   friend class GaiaDR2DatabaseFileGenerator;
   friend class GaiaEDR3DatabaseFileGenerator;
   friend class GaiaDR3DatabaseFileGenerator;
   friend class GaiaDR3SPDatabaseFileGenerator;
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif  // __PCL_GaiaDatabaseFile_h
 
// ----------------------------------------------------------------------------
// EOF pcl/GaiaDatabaseFile.h - Released 2024-06-18T15:48:54Z