File indexing completed on 2024-04-14 14:11:45

 /*
     SPDX-FileCopyrightText: 2010 Akarsh Simha <akarshsimha@gmail.com>
 
     SPDX-License-Identifier: GPL-2.0-or-later
 */
 
 #include "starhopper.h"
 
 #include "kstarsdata.h"
 #include "ksutils.h"
 #include "starcomponent.h"
 #include "skyobjects/starobject.h"
 
 #include <kstars_debug.h>
 
 QList<StarObject *> *StarHopper::computePath(const SkyPoint &src, const SkyPoint &dest, float fov__, float maglim__,
                                              QStringList *metadata_)
 {
     QList<const StarObject *> starHopList_const = computePath_const(src, dest, fov__, maglim__, metadata_);
     QList<StarObject *> *starHopList_unconst    = new QList<StarObject *>();
     foreach (const StarObject *so, starHopList_const)
     {
         starHopList_unconst->append(const_cast<StarObject *>(so));
     }
     return starHopList_unconst;
 }
 
 QList<const StarObject *> StarHopper::computePath_const(const SkyPoint &src, const SkyPoint &dest, float fov_,
                                                         float maglim_, QStringList *metadata)
 {
     fov    = fov_;
     maglim = maglim_;
     start  = &src;
     end    = &dest;
 
     came_from.clear();
     result_path.clear();
 
     // Implements the A* search algorithm
 
     QList<SkyPoint const *> cSet;
     QList<SkyPoint const *> oSet;
     QHash<SkyPoint const *, double> g_score;
     QHash<SkyPoint const *, double> f_score;
     QHash<SkyPoint const *, double> h_score;
 
     qCDebug(KSTARS) << "StarHopper is trying to compute a path from source: " << src.ra().toHMSString()
              << src.dec().toDMSString() << " to destination: " << dest.ra().toHMSString() << dest.dec().toDMSString()
              << "; a starhop of " << src.angularDistanceTo(&dest).Degrees() << " degrees!";
 
     oSet.append(&src);
     g_score[&src] = 0;
     h_score[&src] = src.angularDistanceTo(&dest).Degrees() / fov;
     f_score[&src] = h_score[&src];
 
     while (!oSet.isEmpty())
     {
         qCDebug(KSTARS) << "Next step";
         // Find the node with the lowest f_score value
         SkyPoint const *curr_node = nullptr;
         double lowfscore          = 1.0e8;
 
         foreach (const SkyPoint *sp, oSet)
         {
             if (f_score[sp] < lowfscore)
             {
                 lowfscore = f_score[sp];
                 curr_node = sp;
             }
         }
         if (curr_node == nullptr)
             continue;
 
         qCDebug(KSTARS) << "Lowest fscore (vertex distance-plus-cost score) is " << lowfscore
                  << " with coords: " << curr_node->ra().toHMSString() << curr_node->dec().toDMSString()
                  << ". Considering this node now.";
         if (curr_node == &dest || (curr_node != &src && h_score[curr_node] < 0.5))
         {
             // We are at destination
             reconstructPath(came_from[curr_node]);
             qCDebug(KSTARS) << "We've arrived at the destination! Yay! Result path count: " << result_path.count();
 
             // Just a test -- try to print out useful instructions to the debug console. Once we make star hopper unexperimental, we should move this to some sort of a display
             qCDebug(KSTARS) << "Star Hopping Directions: ";
             const SkyPoint *prevHop = start;
 
             for (auto &hopStar : result_path)
             {
                 QString direction;
                 QString spectralChar = "";
                 double pa; // should be 0 to 2pi
 
                 dms angDist = prevHop->angularDistanceTo(hopStar, &pa);
 
                 dms dmsPA;
                 dmsPA.setRadians(pa);
                 direction = KSUtils::toDirectionString(dmsPA);
 
                 if (metadata)
                 {
                     if (!patternNames.contains(hopStar))
                     {
                         spectralChar += hopStar->spchar();
                         starHopDirections = i18n(" Slew %1 degrees %2 to find an %3 star of mag %4 ", QString::number(angDist.Degrees(), 'f', 2),
                                                  direction, spectralChar, QString::number(hopStar->mag()));
                         qCDebug(KSTARS) << starHopDirections;
                     }
                     else
                     {
                         starHopDirections = i18n(" Slew %1 degrees %2 to find a(n) %3", QString::number(angDist.Degrees(), 'f', 2), direction,
                                                      patternNames.value(hopStar));
                         qCDebug(KSTARS) << starHopDirections;
                     }
                     metadata->append(starHopDirections);
                     qCDebug(KSTARS) << "Appended starHopDirections to metadata";
                 }
                 prevHop = hopStar;
             }
             qCDebug(KSTARS) << "  The destination is within a field-of-view";
 
             return result_path;
         }
 
         oSet.removeOne(curr_node);
         cSet.append(curr_node);
 
         // FIXME: Make sense. If current node ---> dest distance is
         // larger than src --> dest distance by more than 20%, don't
         // even bother considering it.
 
         if (h_score[curr_node] > h_score[&src] * 1.2)
         {
             qCDebug(KSTARS) << "Node under consideration has larger distance to destination (h-score) than start node! "
                         "Ignoring it.";
             continue;
         }
 
         // Get the list of stars that are neighbours of this node
         QList<StarObject *> neighbors;
 
         // FIXME: Actually, this should be done in
         // HorizontalToEquatorial, but we do it here because SkyPoint
         // needs a lot of fixing to handle unprecessed and precessed,
         // equatorial and horizontal coordinates nicely
         SkyPoint *CurrentNode = const_cast<SkyPoint *>(curr_node);
         //CurrentNode->deprecess(KStarsData::Instance()->updateNum());
         CurrentNode->catalogueCoord(KStarsData::Instance()->updateNum()->julianDay());
         //qCDebug(KSTARS) << "Calling starsInAperture";
         StarComponent::Instance()->starsInAperture(neighbors, *curr_node, fov, maglim);
         qCDebug(KSTARS) << "Choosing next node from a set of " << neighbors.count();
         // Look for the potential next node
         double curr_g_score = g_score[curr_node];
 
         for (auto &nhd_node : neighbors)
         {
             if (cSet.contains(nhd_node))
                 continue;
 
             // Compute the tentative g_score
             double tentative_g_score = curr_g_score + cost(curr_node, nhd_node);
             bool tentative_better;
             if (!oSet.contains(nhd_node))
             {
                 oSet.append(nhd_node);
                 tentative_better = true;
             }
             else if (tentative_g_score < g_score[nhd_node])
                 tentative_better = true;
             else
                 tentative_better = false;
 
             if (tentative_better)
             {
                 came_from[nhd_node] = curr_node;
                 g_score[nhd_node]   = tentative_g_score;
                 h_score[nhd_node]   = nhd_node->angularDistanceTo(&dest).Degrees() / fov;
                 f_score[nhd_node]   = g_score[nhd_node] + h_score[nhd_node];
             }
         }
     }
     qCDebug(KSTARS) << "REGRET! Returning empty list!";
     return QList<StarObject const *>(); // Return an empty QList
 }
 
 void StarHopper::reconstructPath(SkyPoint const *curr_node)
 {
     if (curr_node != start)
     {
         StarObject const *s = dynamic_cast<StarObject const *>(curr_node);
         Q_ASSERT(s);
         result_path.prepend(s);
         reconstructPath(came_from[s]);
     }
 }
 
 float StarHopper::cost(const SkyPoint *curr, const SkyPoint *next)
 {
     // This is a very heuristic method, that tries to produce a cost
     // for each hop.
 
     float netcost = 0;
 
     // Convert 'next' into a StarObject
     if (next == start)
     {
         // If the next hop is back to square one, junk it
         return 1e8;
     }
     bool isThisTheEnd = (next == end);
 
     float magcost, speccost;
     magcost = speccost = 0;
 
     if (!isThisTheEnd)
     {
         // We ought to be dealing with a star
         StarObject const *nextstar = dynamic_cast<StarObject const *>(next);
         Q_ASSERT(nextstar);
 
         // Test 1: How bright is the star?
         magcost =
             nextstar->mag() - 7.0 +
             5 * log10(
                     fov); // The brighter, the better. FIXME: 8.0 is now an arbitrary reference to the average faint star. Should actually depend on FOV, something like log( FOV ).
 
         // Test 2: Is the star strikingly red / yellow coloured?
         QString SpType = nextstar->sptype();
         char spclass   = SpType.at(0).toLatin1();
         speccost       = (spclass == 'G' || spclass == 'K' || spclass == 'M') ? -0.3 : 0;
         /*
         // Test 3: Is the star in the general direction of the object?
         // We use the cosine rule to find the angle between the hop direction, and the direction to destination
         // a = side joining curr to end
         // b = side joining curr to next
         // c = side joining next to end
         // C = angle between curr-next and curr-end
         double sina, sinb, cosa, cosb;
         curr->angularDistanceTo(dest).SinCos( &sina, &cosa );
         curr->angularDistanceTo(next).SinCos( &sinb, &cosb );
         double cosc = cos(next->angularDistanceTo(end).radians());
         double cosC = ( cosc - cosa * cosb ) / (sina * sinb);
         float dircost;
         if( cosC < 0 ) // Wrong direction!
         dircost = 1e8; // Some humongous number;
         else
         dircost = sqrt( 1 - cosC * cosC ) / cosC; // tan( C )
         */
     }
 
     // Test 4: How far is the hop?
     double distcost =
         (curr->angularDistanceTo(next).Degrees() /
          fov); // 1 "magnitude" incremental cost for 1 FOV. Is this even required, or is it just equivalent to halving our distance unit? I think it is required since the hop is not necessarily in the direction of the object -- asimha
 
     // Test 5: How effective is the hop? [Might not be required with A*]
     //    double distredcost = -((src->angularDistanceTo( dest ).Degrees() - next->angularDistanceTo( dest ).Degrees()) * 60 / fov)*3; // 3 "magnitudes" for 1 FOV closer
 
     // Test 6: Is the destination an asterism? Are there bright stars clustered nearby?
     QList<StarObject *> localNeighbors;
     StarComponent::Instance()->starsInAperture(localNeighbors, *next, fov / 10, maglim + 1.0);
     double stardensitycost = 1 - localNeighbors.count(); // -1 "magnitude" for every neighbouring star
 
 // Test 7: Identify star patterns
 
 #define RIGHT_ANGLE_THRESHOLD 0.05
 #define EQUAL_EDGE_THRESHOLD  0.025
 
     double patterncost = 0;
     QString patternName;
     if (!isThisTheEnd)
     {
         // We ought to be dealing with a star
         StarObject const *nextstar = dynamic_cast<StarObject const *>(next);
         Q_ASSERT(nextstar);
 
         float factor = 1.0;
         while (factor <= 10.0)
         {
             localNeighbors.clear();
             StarComponent::Instance()->starsInAperture(
                 localNeighbors, *next, fov / factor,
                 nextstar->mag() + 1.0); // Use a larger aperture for pattern identification; max 1.0 mag difference
             foreach (StarObject *star, localNeighbors)
             {
                 if (star == nextstar)
                     localNeighbors.removeOne(star);
                 else if (fabs(star->mag() - nextstar->mag()) > 1.0)
                     localNeighbors.removeOne(star);
             } // Now, we should have a pruned list
             factor += 1.0;
             if (localNeighbors.size() == 2)
                 break;
         }
         factor -= 1.0;
         if (localNeighbors.size() == 2)
         {
             patternName = i18n("triangle (of similar magnitudes)"); // any three stars form a triangle!
             // Try to find triangles. Note that we assume that the standard Euclidian metric works on a sphere for small angles, i.e. the celestial sphere is nearly flat over our FOV.
             StarObject *star1 = localNeighbors[0];
             double dRA1       = nextstar->ra().radians() - star1->ra().radians();
             double dDec1      = nextstar->dec().radians() - star1->dec().radians();
             double dist1sqr   = dRA1 * dRA1 + dDec1 * dDec1;
 
             StarObject *star2 = localNeighbors[1];
             double dRA2       = nextstar->ra().radians() - star2->ra().radians();
             double dDec2      = nextstar->dec().radians() - star2->dec().radians();
             double dist2sqr   = dRA2 * dRA2 + dDec2 * dDec2;
 
             // Check for right-angled triangles (without loss of generality, right angle is at this vertex)
             if (fabs((dRA1 * dRA2 - dDec1 * dDec2) / sqrt(dist1sqr * dist2sqr)) < RIGHT_ANGLE_THRESHOLD)
             {
                 // We have a right angled triangle! Give -3 magnitudes!
                 patterncost += -3;
                 patternName = i18n("right-angled triangle");
             }
 
             // Check for isosceles triangles (without loss of generality, this is the vertex)
             if (fabs((dist1sqr - dist2sqr) / (dist1sqr)) < EQUAL_EDGE_THRESHOLD)
             {
                 patterncost += -1;
                 patternName = i18n("isosceles triangle");
                 if (fabs((dRA2 * dDec1 - dRA1 * dDec2) / sqrt(dist1sqr * dist2sqr)) < RIGHT_ANGLE_THRESHOLD)
                 {
                     patterncost += -1;
                     patternName = i18n("straight line of 3 stars");
                 }
                 // Check for equilateral triangles
                 double dist3    = star1->angularDistanceTo(star2).radians();
                 double dist3sqr = dist3 * dist3;
                 if (fabs((dist3sqr - dist1sqr) / dist1sqr) < EQUAL_EDGE_THRESHOLD)
                 {
                     patterncost += -1;
                     patternName = i18n("equilateral triangle");
                 }
             }
         }
         // TODO: Identify squares.
         if (!patternName.isEmpty())
         {
             patternName += i18n(" within %1% of FOV of the marked star", (int)(100.0 / factor));
             patternNames.insert(nextstar, patternName);
         }
     }
 
     netcost = magcost + speccost + distcost + stardensitycost + patterncost;
     if (netcost < 0)
         netcost = 0.1; // FIXME: Heuristics aren't supposed to be entirely random. This one is.
     qCDebug(KSTARS) << "Mag cost: " << magcost << "; Spec Cost: " << speccost << "; Dist Cost: " << distcost
              << "; Density cost: " << stardensitycost << "; Pattern cost: " << patterncost << "; Net cost: " << netcost
              << "; Pattern: " << patternName;
     return netcost;
 }