1 /* -----------------------------------------------------------------------------
2
3 SPDX-FileCopyrightText: 1997-2016 Martin Reinecke,
4 SPDX-FileCopyrightText: 1997-2016 Krzysztof M. Gorski
5 SPDX-FileCopyrightText: 1997-2016 Eric Hivon
6 SPDX-FileCopyrightText: Benjamin D. Wandelt
7 SPDX-FileCopyrightText: Anthony J. Banday,
8 SPDX-FileCopyrightText: Matthias Bartelmann,
9 SPDX-FileCopyrightText: Reza Ansari
10 SPDX-FileCopyrightText: Kenneth M. Ganga
11
12 This file is part of HEALPix.
13
14 Based on work by Pavel Mraz from SkyTechX.
15
16 Modified by Jasem Mutlaq for KStars:
17 SPDX-FileCopyrightText: Jasem Mutlaq <mutlaqja@ikarustech.com>
18
19 SPDX-License-Identifier: GPL-2.0-or-later
20
21 For more information about HEALPix see https://healpix.sourceforge.io/
22
23 ---------------------------------------------------------------------------*/
24
25 #include "healpix.h"
26
27 #include <QMatrix4x4>
28
29 #include "hipsmanager.h"
30 #include "skypoint.h"
31 #include "kstarsdata.h"
32 #include "geolocation.h"
33
34 static const double twothird = 2.0 / 3.0;
35 static const double twopi = 6.283185307179586476925286766559005768394;
36 static const double inv_halfpi = 0.6366197723675813430755350534900574;
37
38 static const int jrll[] = { 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 };
39 static const int jpll[] = { 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 };
40
41 static const short ctab[] =
42 {
43 #define Z(a) a,a+1,a+256,a+257
44 #define Y(a) Z(a),Z(a+2),Z(a+512),Z(a+514)
45 #define X(a) Y(a),Y(a+4),Y(a+1024),Y(a+1028)
46 X(0), X(8), X(2048), X(2056)
47 #undef X
48 #undef Y
49 #undef Z
50 };
51
52 static const short utab[] =
53 {
54 #define Z(a) 0x##a##0, 0x##a##1, 0x##a##4, 0x##a##5
55 #define Y(a) Z(a##0), Z(a##1), Z(a##4), Z(a##5)
56 #define X(a) Y(a##0), Y(a##1), Y(a##4), Y(a##5)
57 X(0), X(1), X(4), X(5)
58 #undef X
59 #undef Y
60 #undef Z
61 };
62
63 static short xoffset[] = { -1, -1, 0, 1, 1, 1, 0, -1 };
64 static short yoffset[] = { 0, 1, 1, 1, 0, -1, -1, -1 };
65
66 static short facearray[9][12] =
67 {
68 { 8, 9, 10, 11, -1, -1, -1, -1, 10, 11, 8, 9 }, // S
69 { 5, 6, 7, 4, 8, 9, 10, 11, 9, 10, 11, 8 }, // SE
70 { -1, -1, -1, -1, 5, 6, 7, 4, -1, -1, -1, -1 }, // E
71 { 4, 5, 6, 7, 11, 8, 9, 10, 11, 8, 9, 10 }, // SW
72 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 }, // center
73 { 1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 }, // NE
74 { -1, -1, -1, -1, 7, 4, 5, 6, -1, -1, -1, -1 }, // W
75 { 3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 }, // NW
76 { 2, 3, 0, 1, -1, -1, -1, -1, 0, 1, 2, 3 }
77 }; // N
78
79 static uchar swaparray[9][12] =
80 {
81 { 0, 0, 3 }, // S
82 { 0, 0, 6 }, // SE
83 { 0, 0, 0 }, // E
84 { 0, 0, 5 }, // SW
85 { 0, 0, 0 }, // center
86 { 5, 0, 0 }, // NE
87 { 0, 0, 0 }, // W
88 { 6, 0, 0 }, // NW
89 { 3, 0, 0 }
90 }; // N
91
fmodulo(double v1,double v2)92 static double fmodulo(double v1, double v2)
93 {
94 if (v1 >= 0)
95 {
96 return (v1 < v2) ? v1 : fmod(v1, v2);
97 }
98
99 double tmp = fmod(v1, v2) + v2;
100 return (tmp == v2) ? 0. : tmp;
101 }
102
103
getCornerPoints(int level,int pix,SkyPoint * skyCoords)104 void HEALPix::getCornerPoints(int level, int pix, SkyPoint *skyCoords)
105 {
106 QVector3D v[4];
107 // Transform from HealPIX convention to KStars
108 QVector3D transformed[4];
109
110 int nside = 1 << level;
111 boundaries(nside, pix, 1, v);
112
113 // From rectangular coordinates to Sky coordinates
114 for (int i = 0; i < 4; i++)
115 {
116 transformed[i].setX(v[i].x());
117 transformed[i].setY(v[i].y());
118 transformed[i].setZ(v[i].z());
119
120 double ra = 0, de = 0;
121 xyz2sph(transformed[i], ra, de);
122 de /= dms::DegToRad;
123 ra /= dms::DegToRad;
124
125 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME)
126 {
127 dms galacticLong(ra);
128 dms galacticLat(de);
129 skyCoords[i].GalacticToEquatorial1950(&galacticLong, &galacticLat);
130 skyCoords[i].B1950ToJ2000();
131 skyCoords[i].setRA0(skyCoords[i].ra());
132 skyCoords[i].setDec0(skyCoords[i].dec());
133 }
134 else
135 {
136 skyCoords[i].setRA0(ra / 15.0);
137 skyCoords[i].setDec0(de);
138 }
139
140 skyCoords[i].updateCoords(KStarsData::Instance()->updateNum(), false);
141 skyCoords[i].EquatorialToHorizontal(KStarsData::Instance()->lst(), KStarsData::Instance()->geo()->lat());
142 }
143
144 }
145
boundaries(qint32 nside,qint32 pix,int step,QVector3D * out)146 void HEALPix::boundaries(qint32 nside, qint32 pix, int step, QVector3D *out)
147 {
148 int ix, iy, fn;
149
150 nest2xyf(nside, pix, &ix, &iy, &fn);
151
152 double dc = 0.5 / nside;
153 double xc = (ix + 0.5) / nside;
154 double yc = (iy + 0.5) / nside;
155 double d = 1. / (step * nside);
156
157 for (int i = 0; i < step; ++i)
158 {
159 out[0] = toVec3(xc + dc - i * d, yc + dc, fn);
160 out[1] = toVec3(xc - dc, yc + dc - i * d, fn);
161 out[2] = toVec3(xc - dc + i * d, yc - dc, fn);
162 out[3] = toVec3(xc + dc, yc - dc + i * d, fn);
163 }
164 }
165
toVec3(double fx,double fy,int face)166 QVector3D HEALPix::toVec3(double fx, double fy, int face)
167 {
168 double jr = jrll[face] - fx - fy;
169
170 double locz;
171 double locsth;
172 double locphi;
173 bool lochave_sth = false;
174
175 double nr;
176 if (jr < 1)
177 {
178 nr = jr;
179 double tmp = nr * nr / 3.;
180 locz = 1 - tmp;
181 if (locz > 0.99)
182 {
183 locsth = sqrt(tmp * (2. - tmp));
184 lochave_sth = true;
185 }
186 }
187 else if (jr > 3)
188 {
189 nr = 4 - jr;
190 double tmp = nr * nr / 3.;
191 locz = tmp - 1;
192 if (locz < -0.99)
193 {
194 locsth = sqrt(tmp * (2. - tmp));
195 lochave_sth = true;
196 }
197 }
198 else
199 {
200 nr = 1;
201 locz = (2 - jr) * 2. / 3.;
202 }
203
204 double tmp = jpll[face] * nr + fx - fy;
205 if (tmp < 0) tmp += 8;
206 if (tmp >= 8) tmp -= 8;
207 locphi = (nr < 1e-15) ? 0 : (0.5 * dms::PI / 2.0 * tmp) / nr;
208
209 double st = lochave_sth ? locsth : sqrt((1.0 - locz) * (1.0 + locz));
210
211 QVector3D out;
212
213 out.setX(st * cos(locphi));
214 out.setY(st * sin(locphi));
215 out.setZ(locz);
216
217 return out;
218 }
219
nest2xyf(int nside,int pix,int * ix,int * iy,int * face_num)220 void HEALPix::nest2xyf(int nside, int pix, int *ix, int *iy, int *face_num)
221 {
222 int npface_ = nside * nside, raw;
223 *face_num = pix / npface_;
224 pix &= (npface_ - 1);
225 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
226 *ix = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
227 pix >>= 1;
228 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
229 *iy = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
230 }
231
spread_bits64(int v)232 static int64_t spread_bits64 (int v)
233 {
234 return (int64_t)(utab[ v & 0xff])
235 | ((int64_t)(utab[(v >> 8) & 0xff]) << 16)
236 | ((int64_t)(utab[(v >> 16) & 0xff]) << 32)
237 | ((int64_t)(utab[(v >> 24) & 0xff]) << 48);
238 }
239
240
xyf2nest(int nside,int ix,int iy,int face_num)241 static int xyf2nest (int nside, int ix, int iy, int face_num)
242 {
243 return (face_num * nside * nside) +
244 (utab[ix & 0xff] | (utab[ix >> 8] << 16)
245 | (utab[iy & 0xff] << 1) | (utab[iy >> 8] << 17));
246 }
247
248
249 /*
250 int static leadingZeros(qint32 value)
251 {
252 int leadingZeros = 0;
253 while(value != 0)
254 {
255 value = value >> 1;
256 leadingZeros++;
257 }
258
259 return (32 - leadingZeros);
260 }
261 */
262
263 /*
264 static int ilog2(qint32 arg)
265 {
266 return 32 - leadingZeros(qMax(arg, 1));
267 }
268 */
269
nside2order(qint32 nside)270 static int nside2order(qint32 nside)
271 {
272 {
273 int i = 0;
274 while((nside >> (++i)) > 0);
275 return --i;
276 }
277 }
278
279 /** Returns the neighboring pixels of ipix.
280 This method works in both RING and NEST schemes, but is
281 considerably faster in the NEST scheme.
282 @param nside defines the size of the neighboring area
283 @param ipix the requested pixel number.
284 @param result the result array
285 @return array with indices of the neighboring pixels.
286 The returned array contains (in this order)
287 the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor
288 of ipix. If a neighbor does not exist (this can only happen
289 for the W, N, E and S neighbors), its entry is set to -1. */
290
neighbours(int nside,qint32 ipix,int * result)291 void HEALPix::neighbours(int nside, qint32 ipix, int *result)
292 {
293 int ix, iy, face_num;
294
295 int order = nside2order(nside);
296
297 nest2xyf(nside, ipix, &ix, &iy, &face_num);
298
299 qint32 nsm1 = nside - 1;
300 if ((ix > 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1))
301 {
302 qint32 fpix = (qint32)(face_num) << (2 * order),
303 px0 = spread_bits64(ix ), py0 = spread_bits64(iy ) << 1,
304 pxp = spread_bits64(ix + 1), pyp = spread_bits64(iy + 1) << 1,
305 pxm = spread_bits64(ix - 1), pym = spread_bits64(iy - 1) << 1;
306
307 result[0] = fpix + pxm + py0;
308 result[1] = fpix + pxm + pyp;
309 result[2] = fpix + px0 + pyp;
310 result[3] = fpix + pxp + pyp;
311 result[4] = fpix + pxp + py0;
312 result[5] = fpix + pxp + pym;
313 result[6] = fpix + px0 + pym;
314 result[7] = fpix + pxm + pym;
315 }
316 else
317 {
318 for (int i = 0; i < 8; ++i)
319 {
320 int x = ix + xoffset[i];
321 int y = iy + yoffset[i];
322 int nbnum = 4;
323 if (x < 0)
324 {
325 x += nside;
326 nbnum -= 1;
327 }
328 else if (x >= nside)
329 {
330 x -= nside;
331 nbnum += 1;
332 }
333 if (y < 0)
334 {
335 y += nside;
336 nbnum -= 3;
337 }
338 else if (y >= nside)
339 {
340 y -= nside;
341 nbnum += 3;
342 }
343
344 int f = facearray[nbnum][face_num];
345
346 if (f >= 0)
347 {
348 int bits = swaparray[nbnum][face_num >> 2];
349 if ((bits & 1) > 0) x = (int)(nside - x - 1);
350 if ((bits & 2) > 0) y = (int)(nside - y - 1);
351 if ((bits & 4) > 0)
352 {
353 int tint = x;
354 x = y;
355 y = tint;
356 }
357 result[i] = xyf2nest(nside, x, y, f);
358 }
359 else
360 result[i] = -1;
361 }
362 }
363 }
364
ang2pix_nest_z_phi(qint32 nside_,double z,double phi)365 int HEALPix::ang2pix_nest_z_phi (qint32 nside_, double z, double phi)
366 {
367 double za = fabs(z);
368 double tt = fmodulo(phi, twopi) * inv_halfpi; /* in [0,4) */
369 int face_num, ix, iy;
370
371 if (za <= twothird) /* Equatorial region */
372 {
373 double temp1 = nside_ * (0.5 + tt);
374 double temp2 = nside_ * (z * 0.75);
375 int jp = (int)(temp1 - temp2); /* index of ascending edge line */
376 int jm = (int)(temp1 + temp2); /* index of descending edge line */
377 int ifp = jp / nside_; /* in {0,4} */
378 int ifm = jm / nside_;
379 face_num = (ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8));
380
381 ix = jm & (nside_ - 1);
382 iy = nside_ - (jp & (nside_ - 1)) - 1;
383 }
384 else /* polar region, za > 2/3 */
385 {
386 int ntt = (int)tt, jp, jm;
387 double tp, tmp;
388 if (ntt >= 4) ntt = 3;
389 tp = tt - ntt;
390 tmp = nside_ * sqrt(3 * (1 - za));
391
392 jp = (int)(tp * tmp); /* increasing edge line index */
393 jm = (int)((1.0 - tp) * tmp); /* decreasing edge line index */
394 if (jp >= nside_) jp = nside_ - 1; /* for points too close to the boundary */
395 if (jm >= nside_) jm = nside_ - 1;
396 if (z >= 0)
397 {
398 face_num = ntt; /* in {0,3} */
399 ix = nside_ - jm - 1;
400 iy = nside_ - jp - 1;
401 }
402 else
403 {
404 face_num = ntt + 8; /* in {8,11} */
405 ix = jp;
406 iy = jm;
407 }
408 }
409
410 return xyf2nest(nside_, ix, iy, face_num);
411 }
412
getPix(int level,double ra,double dec)413 int HEALPix::getPix(int level, double ra, double dec)
414 {
415 int nside = 1 << level;
416 double polar[2] = { 0, 0 };
417
418 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME)
419 {
420 static QMatrix4x4 gl(-0.0548762f, -0.873437f, -0.483835f, 0,
421 0.4941100f, -0.444830f, 0.746982f, 0,
422 -0.8676660f, -0.198076f, 0.455984f, 0,
423 0, 0, 0, 1);
424
425 double rcb = cos(dec);
426 QVector3D xyz = QVector3D(rcb * cos(ra), rcb * sin(ra), sin(dec));
427
428 xyz = gl.mapVector(xyz);
429 xyz2sph(xyz, polar[1], polar[0]);
430 }
431 else
432 {
433 polar[0] = dec;
434 polar[1] = ra;
435 }
436
437 return ang2pix_nest_z_phi(nside, sin(polar[0]), polar[1]);
438 }
439
getPixChilds(int pix,int * childs)440 void HEALPix::getPixChilds(int pix, int *childs)
441 {
442 childs[0] = pix * 4 + 0;
443 childs[1] = pix * 4 + 1;
444 childs[2] = pix * 4 + 2;
445 childs[3] = pix * 4 + 3;
446 }
447
xyz2sph(const QVector3D & vec,double & l,double & b)448 void HEALPix::xyz2sph(const QVector3D &vec, double &l, double &b)
449 {
450 double rho = vec.x() * vec.x() + vec.y() * vec.y();
451
452 if (rho > 0)
453 {
454 l = atan2(vec.y(), vec.x());
455 l -= floor(l / (M_PI * 2)) * (M_PI * 2);
456 b = atan2(vec.z(), sqrt(rho));
457 }
458 else
459 {
460 l = 0.0;
461 if (vec.z() == 0.0)
462 {
463 b = 0.0;
464 }
465 else
466 {
467 b = (vec.z() > 0.0) ? M_PI / 2. : -dms::PI / 2.;
468 }
469 }
470 }
471
472
473
474