Kstars

starhopper.cpp
1/*
2 SPDX-FileCopyrightText: 2010 Akarsh Simha <akarshsimha@gmail.com>
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
16QList<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
28QList<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
185void 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
196float 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}
float mag() const
Definition skyobject.h:207
The sky coordinates of a point in the sky.
Definition skypoint.h:45
const CachingDms & dec() const
Definition skypoint.h:269
const CachingDms & ra() const
Definition skypoint.h:263
dms angularDistanceTo(const SkyPoint *sp, double *const positionAngle=nullptr) const
Computes the angular distance between two SkyObjects.
Definition skypoint.cpp:899
SkyPoint catalogueCoord(long double jdf)
Computes the J2000.0 catalogue coordinates for this SkyPoint using the epoch removing aberration,...
Definition skypoint.cpp:710
static StarComponent * Instance()
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...
QList< StarObject * > * computePath(const SkyPoint &src, const SkyPoint &dest, float fov__, float maglim__, QStringList *metadata_=nullptr)
Computes path for Star Hop.
This is a subclass of SkyObject.
Definition starobject.h:33
QString sptype(void) const
Returns entire spectral type string.
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
const QString toDMSString(const bool forceSign=false, const bool machineReadable=false, const bool highPrecision=false) const
Definition dms.cpp:287
const QString toHMSString(const bool machineReadable=false, const bool highPrecision=false) const
Definition dms.cpp:378
double radians() const
Express the angle in radians.
Definition dms.h:325
virtual void setRadians(const double &Rad)
Set angle according to the argument, in radians.
Definition dms.h:333
const double & Degrees() const
Definition dms.h:141
QString i18n(const char *text, const TYPE &arg...)
const QList< QKeySequence > & next()
char toLatin1() const const
void clear()
bool contains(const Key &key) const const
iterator insert(const Key &key, const T &value)
T value(const Key &key) const const
void append(QList< T > &&value)
void clear()
bool contains(const AT &value) const const
qsizetype count() const const
bool isEmpty() const const
bool removeOne(const AT &t)
qsizetype size() const const
const QChar at(qsizetype position) const const
bool isEmpty() const const
QString number(double n, char format, int precision)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Sat Dec 21 2024 17:04:47 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.