Kstars

starhopper.cpp
1 /*
2  SPDX-FileCopyrightText: 2010 Akarsh Simha <[email protected]>
3 
4  SPDX-License-Identifier: GPL-2.0-or-later
5 */
6 
7 #include "starhopper.h"
8 
9 #include "kstarsdata.h"
10 #include "ksutils.h"
11 #include "starcomponent.h"
12 #include "skyobjects/starobject.h"
13 
14 #include <kstars_debug.h>
15 
16 QList<StarObject *> *StarHopper::computePath(const SkyPoint &src, const SkyPoint &dest, float fov__, float maglim__,
17  QStringList *metadata_)
18 {
19  QList<const StarObject *> starHopList_const = computePath_const(src, dest, fov__, maglim__, metadata_);
20  QList<StarObject *> *starHopList_unconst = new QList<StarObject *>();
21  foreach (const StarObject *so, starHopList_const)
22  {
23  starHopList_unconst->append(const_cast<StarObject *>(so));
24  }
25  return starHopList_unconst;
26 }
27 
28 QList<const StarObject *> StarHopper::computePath_const(const SkyPoint &src, const SkyPoint &dest, float fov_,
29  float maglim_, QStringList *metadata)
30 {
31  fov = fov_;
32  maglim = maglim_;
33  start = &src;
34  end = &dest;
35 
36  came_from.clear();
37  result_path.clear();
38 
39  // Implements the A* search algorithm
40 
46 
47  qCDebug(KSTARS) << "StarHopper is trying to compute a path from source: " << src.ra().toHMSString()
48  << src.dec().toDMSString() << " to destination: " << dest.ra().toHMSString() << dest.dec().toDMSString()
49  << "; a starhop of " << src.angularDistanceTo(&dest).Degrees() << " degrees!";
50 
51  oSet.append(&src);
52  g_score[&src] = 0;
53  h_score[&src] = src.angularDistanceTo(&dest).Degrees() / fov;
54  f_score[&src] = h_score[&src];
55 
56  while (!oSet.isEmpty())
57  {
58  qCDebug(KSTARS) << "Next step";
59  // Find the node with the lowest f_score value
60  SkyPoint const *curr_node = nullptr;
61  double lowfscore = 1.0e8;
62 
63  foreach (const SkyPoint *sp, oSet)
64  {
65  if (f_score[sp] < lowfscore)
66  {
67  lowfscore = f_score[sp];
68  curr_node = sp;
69  }
70  }
71  if (curr_node == nullptr)
72  continue;
73 
74  qCDebug(KSTARS) << "Lowest fscore (vertex distance-plus-cost score) is " << lowfscore
75  << " with coords: " << curr_node->ra().toHMSString() << curr_node->dec().toDMSString()
76  << ". Considering this node now.";
77  if (curr_node == &dest || (curr_node != &src && h_score[curr_node] < 0.5))
78  {
79  // We are at destination
80  reconstructPath(came_from[curr_node]);
81  qCDebug(KSTARS) << "We've arrived at the destination! Yay! Result path count: " << result_path.count();
82 
83  // 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
84  qCDebug(KSTARS) << "Star Hopping Directions: ";
85  const SkyPoint *prevHop = start;
86 
87  for (auto &hopStar : result_path)
88  {
89  QString direction;
90  QString spectralChar = "";
91  double pa; // should be 0 to 2pi
92 
93  dms angDist = prevHop->angularDistanceTo(hopStar, &pa);
94 
95  dms dmsPA;
96  dmsPA.setRadians(pa);
97  direction = KSUtils::toDirectionString(dmsPA);
98 
99  if (metadata)
100  {
101  if (!patternNames.contains(hopStar))
102  {
103  spectralChar += hopStar->spchar();
104  starHopDirections = i18n(" Slew %1 degrees %2 to find an %3 star of mag %4 ", QString::number(angDist.Degrees(), 'f', 2),
105  direction, spectralChar, QString::number(hopStar->mag()));
106  qCDebug(KSTARS) << starHopDirections;
107  }
108  else
109  {
110  starHopDirections = i18n(" Slew %1 degrees %2 to find a(n) %3", QString::number(angDist.Degrees(), 'f', 2), direction,
111  patternNames.value(hopStar));
112  qCDebug(KSTARS) << starHopDirections;
113  }
114  metadata->append(starHopDirections);
115  qCDebug(KSTARS) << "Appended starHopDirections to metadata";
116  }
117  prevHop = hopStar;
118  }
119  qCDebug(KSTARS) << " The destination is within a field-of-view";
120 
121  return result_path;
122  }
123 
124  oSet.removeOne(curr_node);
125  cSet.append(curr_node);
126 
127  // FIXME: Make sense. If current node ---> dest distance is
128  // larger than src --> dest distance by more than 20%, don't
129  // even bother considering it.
130 
131  if (h_score[curr_node] > h_score[&src] * 1.2)
132  {
133  qCDebug(KSTARS) << "Node under consideration has larger distance to destination (h-score) than start node! "
134  "Ignoring it.";
135  continue;
136  }
137 
138  // Get the list of stars that are neighbours of this node
139  QList<StarObject *> neighbors;
140 
141  // FIXME: Actually, this should be done in
142  // HorizontalToEquatorial, but we do it here because SkyPoint
143  // needs a lot of fixing to handle unprecessed and precessed,
144  // equatorial and horizontal coordinates nicely
145  SkyPoint *CurrentNode = const_cast<SkyPoint *>(curr_node);
146  //CurrentNode->deprecess(KStarsData::Instance()->updateNum());
147  CurrentNode->catalogueCoord(KStarsData::Instance()->updateNum()->julianDay());
148  //qCDebug(KSTARS) << "Calling starsInAperture";
149  StarComponent::Instance()->starsInAperture(neighbors, *curr_node, fov, maglim);
150  qCDebug(KSTARS) << "Choosing next node from a set of " << neighbors.count();
151  // Look for the potential next node
152  double curr_g_score = g_score[curr_node];
153 
154  for (auto &nhd_node : neighbors)
155  {
156  if (cSet.contains(nhd_node))
157  continue;
158 
159  // Compute the tentative g_score
160  double tentative_g_score = curr_g_score + cost(curr_node, nhd_node);
161  bool tentative_better;
162  if (!oSet.contains(nhd_node))
163  {
164  oSet.append(nhd_node);
165  tentative_better = true;
166  }
167  else if (tentative_g_score < g_score[nhd_node])
168  tentative_better = true;
169  else
170  tentative_better = false;
171 
172  if (tentative_better)
173  {
174  came_from[nhd_node] = curr_node;
175  g_score[nhd_node] = tentative_g_score;
176  h_score[nhd_node] = nhd_node->angularDistanceTo(&dest).Degrees() / fov;
177  f_score[nhd_node] = g_score[nhd_node] + h_score[nhd_node];
178  }
179  }
180  }
181  qCDebug(KSTARS) << "REGRET! Returning empty list!";
182  return QList<StarObject const *>(); // Return an empty QList
183 }
184 
185 void StarHopper::reconstructPath(SkyPoint const *curr_node)
186 {
187  if (curr_node != start)
188  {
189  StarObject const *s = dynamic_cast<StarObject const *>(curr_node);
190  Q_ASSERT(s);
191  result_path.prepend(s);
192  reconstructPath(came_from[s]);
193  }
194 }
195 
196 float StarHopper::cost(const SkyPoint *curr, const SkyPoint *next)
197 {
198  // This is a very heuristic method, that tries to produce a cost
199  // for each hop.
200 
201  float netcost = 0;
202 
203  // Convert 'next' into a StarObject
204  if (next == start)
205  {
206  // If the next hop is back to square one, junk it
207  return 1e8;
208  }
209  bool isThisTheEnd = (next == end);
210 
211  float magcost, speccost;
212  magcost = speccost = 0;
213 
214  if (!isThisTheEnd)
215  {
216  // We ought to be dealing with a star
217  StarObject const *nextstar = dynamic_cast<StarObject const *>(next);
218  Q_ASSERT(nextstar);
219 
220  // Test 1: How bright is the star?
221  magcost =
222  nextstar->mag() - 7.0 +
223  5 * log10(
224  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 ).
225 
226  // Test 2: Is the star strikingly red / yellow coloured?
227  QString SpType = nextstar->sptype();
228  char spclass = SpType.at(0).toLatin1();
229  speccost = (spclass == 'G' || spclass == 'K' || spclass == 'M') ? -0.3 : 0;
230  /*
231  // Test 3: Is the star in the general direction of the object?
232  // We use the cosine rule to find the angle between the hop direction, and the direction to destination
233  // a = side joining curr to end
234  // b = side joining curr to next
235  // c = side joining next to end
236  // C = angle between curr-next and curr-end
237  double sina, sinb, cosa, cosb;
238  curr->angularDistanceTo(dest).SinCos( &sina, &cosa );
239  curr->angularDistanceTo(next).SinCos( &sinb, &cosb );
240  double cosc = cos(next->angularDistanceTo(end).radians());
241  double cosC = ( cosc - cosa * cosb ) / (sina * sinb);
242  float dircost;
243  if( cosC < 0 ) // Wrong direction!
244  dircost = 1e8; // Some humongous number;
245  else
246  dircost = sqrt( 1 - cosC * cosC ) / cosC; // tan( C )
247  */
248  }
249 
250  // Test 4: How far is the hop?
251  double distcost =
252  (curr->angularDistanceTo(next).Degrees() /
253  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
254 
255  // Test 5: How effective is the hop? [Might not be required with A*]
256  // double distredcost = -((src->angularDistanceTo( dest ).Degrees() - next->angularDistanceTo( dest ).Degrees()) * 60 / fov)*3; // 3 "magnitudes" for 1 FOV closer
257 
258  // Test 6: Is the destination an asterism? Are there bright stars clustered nearby?
259  QList<StarObject *> localNeighbors;
260  StarComponent::Instance()->starsInAperture(localNeighbors, *next, fov / 10, maglim + 1.0);
261  double stardensitycost = 1 - localNeighbors.count(); // -1 "magnitude" for every neighbouring star
262 
263 // Test 7: Identify star patterns
264 
265 #define RIGHT_ANGLE_THRESHOLD 0.05
266 #define EQUAL_EDGE_THRESHOLD 0.025
267 
268  double patterncost = 0;
269  QString patternName;
270  if (!isThisTheEnd)
271  {
272  // We ought to be dealing with a star
273  StarObject const *nextstar = dynamic_cast<StarObject const *>(next);
274  Q_ASSERT(nextstar);
275 
276  float factor = 1.0;
277  while (factor <= 10.0)
278  {
279  localNeighbors.clear();
281  localNeighbors, *next, fov / factor,
282  nextstar->mag() + 1.0); // Use a larger aperture for pattern identification; max 1.0 mag difference
283  foreach (StarObject *star, localNeighbors)
284  {
285  if (star == nextstar)
286  localNeighbors.removeOne(star);
287  else if (fabs(star->mag() - nextstar->mag()) > 1.0)
288  localNeighbors.removeOne(star);
289  } // Now, we should have a pruned list
290  factor += 1.0;
291  if (localNeighbors.size() == 2)
292  break;
293  }
294  factor -= 1.0;
295  if (localNeighbors.size() == 2)
296  {
297  patternName = i18n("triangle (of similar magnitudes)"); // any three stars form a triangle!
298  // 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.
299  StarObject *star1 = localNeighbors[0];
300  double dRA1 = nextstar->ra().radians() - star1->ra().radians();
301  double dDec1 = nextstar->dec().radians() - star1->dec().radians();
302  double dist1sqr = dRA1 * dRA1 + dDec1 * dDec1;
303 
304  StarObject *star2 = localNeighbors[1];
305  double dRA2 = nextstar->ra().radians() - star2->ra().radians();
306  double dDec2 = nextstar->dec().radians() - star2->dec().radians();
307  double dist2sqr = dRA2 * dRA2 + dDec2 * dDec2;
308 
309  // Check for right-angled triangles (without loss of generality, right angle is at this vertex)
310  if (fabs((dRA1 * dRA2 - dDec1 * dDec2) / sqrt(dist1sqr * dist2sqr)) < RIGHT_ANGLE_THRESHOLD)
311  {
312  // We have a right angled triangle! Give -3 magnitudes!
313  patterncost += -3;
314  patternName = i18n("right-angled triangle");
315  }
316 
317  // Check for isosceles triangles (without loss of generality, this is the vertex)
318  if (fabs((dist1sqr - dist2sqr) / (dist1sqr)) < EQUAL_EDGE_THRESHOLD)
319  {
320  patterncost += -1;
321  patternName = i18n("isosceles triangle");
322  if (fabs((dRA2 * dDec1 - dRA1 * dDec2) / sqrt(dist1sqr * dist2sqr)) < RIGHT_ANGLE_THRESHOLD)
323  {
324  patterncost += -1;
325  patternName = i18n("straight line of 3 stars");
326  }
327  // Check for equilateral triangles
328  double dist3 = star1->angularDistanceTo(star2).radians();
329  double dist3sqr = dist3 * dist3;
330  if (fabs((dist3sqr - dist1sqr) / dist1sqr) < EQUAL_EDGE_THRESHOLD)
331  {
332  patterncost += -1;
333  patternName = i18n("equilateral triangle");
334  }
335  }
336  }
337  // TODO: Identify squares.
338  if (!patternName.isEmpty())
339  {
340  patternName += i18n(" within %1% of FOV of the marked star", (int)(100.0 / factor));
341  patternNames.insert(nextstar, patternName);
342  }
343  }
344 
345  netcost = magcost + speccost + distcost + stardensitycost + patterncost;
346  if (netcost < 0)
347  netcost = 0.1; // FIXME: Heuristics aren't supposed to be entirely random. This one is.
348  qCDebug(KSTARS) << "Mag cost: " << magcost << "; Spec Cost: " << speccost << "; Dist Cost: " << distcost
349  << "; Density cost: " << stardensitycost << "; Pattern cost: " << patterncost << "; Net cost: " << netcost
350  << "; Pattern: " << patternName;
351  return netcost;
352 }
void append(const T &value)
const T value(const Key &key) const const
QString number(int n, int base)
Stores dms coordinates for a point in the sky. for converting between coordinate systems.
Definition: skypoint.h:44
int count(const T &value) const const
void clear()
static StarComponent * Instance()
Definition: starcomponent.h:64
float mag() const
Definition: skyobject.h:206
bool contains(const T &value) const const
const QString toHMSString(const bool machineReadable=false, const bool highPrecision=false) const
Definition: dms.cpp:370
QHash::iterator insert(const Key &key, const T &value)
QList< StarObject * > * computePath(const SkyPoint &src, const SkyPoint &dest, float fov__, float maglim__, QStringList *metadata_=nullptr)
Computes path for Star Hop.
Definition: starhopper.cpp:16
int size() const const
QString i18n(const char *text, const TYPE &arg...)
const QString toDMSString(const bool forceSign=false, const bool machineReadable=false, const bool highPrecision=false) const
Definition: dms.cpp:279
const CachingDms & dec() const
Definition: skypoint.h:269
subclass of SkyObject specialized for stars.
Definition: starobject.h:32
bool removeOne(const T &value)
bool isEmpty() const const
QString sptype(void) const
Returns entire spectral type string.
Definition: starobject.cpp:549
bool isEmpty() const const
dms angularDistanceTo(const SkyPoint *sp, double *const positionAngle=nullptr) const
Computes the angular distance between two SkyObjects.
Definition: skypoint.cpp:899
An angle, stored as degrees, but expressible in many ways.
Definition: dms.h:37
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition: dms.h:333
const CachingDms & ra() const
Definition: skypoint.h:263
double radians() const
Express the angle in radians.
Definition: dms.h:325
const double & Degrees() const
Definition: dms.h:141
const QChar at(int position) const const
void clear()
bool contains(const Key &key) const const
void starsInAperture(QList< StarObject * > &list, const SkyPoint &center, float radius, float maglim=-29)
Add to the given list, the stars from this component, that lie within the specified circular aperture...
char toLatin1() const const
QAction * next(const QObject *recvr, const char *slot, QObject *parent)
SkyPoint catalogueCoord(long double jdf)
Computes the J2000.0 catalogue coordinates for this SkyPoint using the epoch removing aberration,...
Definition: skypoint.cpp:710
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Fri Jun 9 2023 04:02:26 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.