diff --git a/.gitignore b/.gitignore index 0d545c9..8e04d76 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ node_modules/ .npmrc *.log .DS_Store +dist \ No newline at end of file diff --git a/package-lock.json b/package-lock.json index 0825074..12e8023 100644 --- a/package-lock.json +++ b/package-lock.json @@ -1,12 +1,12 @@ { "name": "@gkucmierz/healpixjs-bigint", - "version": "1.0.0", + "version": "1.0.2", "lockfileVersion": 3, "requires": true, "packages": { "": { "name": "@gkucmierz/healpixjs-bigint", - "version": "1.0.0", + "version": "1.0.2", "license": "MIT" } } diff --git a/package.json b/package.json index 64fe9eb..7dfe745 100644 --- a/package.json +++ b/package.json @@ -1,6 +1,6 @@ { "name": "@gkucmierz/healpixjs-bigint", - "version": "1.0.1", + "version": "1.0.2", "description": "True BigInt port of healpixjs supporting orders beyond the 32-bit (<=14) limit", "keywords": [ "healpix", @@ -31,4 +31,4 @@ }, "author": "Grzegorz Kućmierz", "license": "MIT" -} \ No newline at end of file +} diff --git a/src/Healpix.js b/src/Healpix.js index 157b592..7b72d15 100644 --- a/src/Healpix.js +++ b/src/Healpix.js @@ -23,136 +23,136 @@ import { Zphi } from "./Zphi.js"; // import Constants from './Constants.js'; // import RangeSet from './RangeSet.js'; export class Healpix { - constructor(nside_in) { - this.order_max = 29; - this.inv_halfpi = 2.0 / Math.PI; - this.twothird = 2.0 / 3.; - // console.log("twothird "+this.twothird); - // this.ns_max=1L< 0) { - this.nside = nside_in; - this.npface = BigInt(this.nside) * BigInt(this.nside); - this.npix = 12n * this.npface; - this.order = this.nside2order(this.nside); - this.nl2 = 2 * this.nside; - this.nl3 = 3 * this.nside; - this.nl4 = 4 * this.nside; - this.fact2 = 4.0 / Number(this.npix); - this.fact1 = (this.nside << 1) * this.fact2; - this.ncap = 2 * this.nside * (this.nside - 1); // pixels in each polar cap - // console.log("order: "+this.order); - // console.log("nside: "+this.nside); - } - this.bn = []; - this.mpr = []; - this.cmpr = []; - this.smpr = []; - // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT - // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT - // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT - // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT - // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT - // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT - // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT - // Uncaught RangeError: Maximum call stack size exceeded - // MOVED TO computeBn() - // for (let i=0; i <= this.order_max; ++i) { - // this.bn[i]=new Healpix(1< 0) { + this.nside = nside_in; + this.npface = BigInt(this.nside) * BigInt(this.nside); + this.npix = 12n * this.npface; + this.order = this.nside2order(this.nside); + this.nl2 = 2 * this.nside; + this.nl3 = 3 * this.nside; + this.nl4 = 4 * this.nside; + this.fact2 = 4.0 / Number(this.npix); + this.fact1 = (this.nside << 1) * this.fact2; + this.ncap = 2 * this.nside * (this.nside - 1); // pixels in each polar cap + // console.log("order: "+this.order); + // console.log("nside: "+this.nside); } - computeBn() { - for (let i = 0; i <= this.order_max; ++i) { - this.bn[i] = new Healpix(1 << i); - this.mpr[i] = this.bn[i].maxPixrad(); - this.cmpr[i] = Hploc.cos(this.mpr[i]); - this.smpr[i] = Hploc.sin(this.mpr[i]); - } + this.bn = []; + this.mpr = []; + this.cmpr = []; + this.smpr = []; + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // Uncaught RangeError: Maximum call stack size exceeded + // MOVED TO computeBn() + // for (let i=0; i <= this.order_max; ++i) { + // this.bn[i]=new Healpix(1< 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1)) { - let fpix = Math.floor(face_num << (2 * this.order)); - let px0 = this.spread_bits(ix); - let py0 = this.spread_bits(iy) << 1; - let pxp = this.spread_bits(ix + 1); - let pyp = this.spread_bits(iy + 1) << 1; - let pxm = this.spread_bits(ix - 1); - let pym = this.spread_bits(iy - 1) << 1; - result[0] = fpix + pxm + py0; - result[1] = fpix + pxm + pyp; - result[2] = fpix + px0 + pyp; - result[3] = fpix + pxp + pyp; - result[4] = fpix + pxp + py0; - result[5] = fpix + pxp + pym; - result[6] = fpix + px0 + pym; - result[7] = fpix + pxm + pym; + neighbours(ipix) { + let result = new BigInt64Array(8); + let xyf = this.nest2xyf(ipix); + let ix = xyf.ix; + let iy = xyf.iy; + let face_num = xyf.face; + var nsm1 = this.nside - 1; + if ((ix > 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1)) { + let fpix = BigInt(face_num) << (2n * BigInt(this.order)); + let px0 = this.spread_bits(ix); + let py0 = this.spread_bits(iy) << 1n; + let pxp = this.spread_bits(ix + 1); + let pyp = this.spread_bits(iy + 1) << 1n; + let pxm = this.spread_bits(ix - 1); + let pym = this.spread_bits(iy - 1) << 1n; + result[0] = fpix + pxm + py0; + result[1] = fpix + pxm + pyp; + result[2] = fpix + px0 + pyp; + result[3] = fpix + pxp + pyp; + result[4] = fpix + pxp + py0; + result[5] = fpix + pxp + pym; + result[6] = fpix + px0 + pym; + result[7] = fpix + pxm + pym; + } + else { + for (let i = 0; i < 8; ++i) { + let x = ix + this.xoffset[i]; + let y = iy + this.yoffset[i]; + let nbnum = 4; + if (x < 0) { + x += this.nside; + nbnum -= 1; + } + else if (x >= this.nside) { + x -= this.nside; + nbnum += 1; + } + if (y < 0) { + y += this.nside; + nbnum -= 3; + } + else if (y >= this.nside) { + y -= this.nside; + nbnum += 3; + } + let f = this.facearray[nbnum][face_num]; + if (f >= 0) { + let bits = this.swaparray[nbnum][face_num >>> 2]; + if ((bits & 1) > 0) { + x = Math.floor(this.nside - x - 1); + } + if ((bits & 2) > 0) { + y = Math.floor(this.nside - y - 1); + } + if ((bits & 4) > 0) { + let tint = x; + x = y; + y = tint; + } + result[i] = this.xyf2nest(BigInt(x), BigInt(y), BigInt(f)); } else { - for (let i = 0; i < 8; ++i) { - let x = ix + this.xoffset[i]; - let y = iy + this.yoffset[i]; - let nbnum = 4; - if (x < 0) { - x += this.nside; - nbnum -= 1; - } - else if (x >= this.nside) { - x -= this.nside; - nbnum += 1; - } - if (y < 0) { - y += this.nside; - nbnum -= 3; - } - else if (y >= this.nside) { - y -= this.nside; - nbnum += 3; - } - let f = this.facearray[nbnum][face_num]; - if (f >= 0) { - let bits = this.swaparray[nbnum][face_num >>> 2]; - if ((bits & 1) > 0) { - x = Math.floor(this.nside - x - 1); - } - if ((bits & 2) > 0) { - y = Math.floor(this.nside - y - 1); - } - if ((bits & 4) > 0) { - let tint = x; - x = y; - y = tint; - } - result[i] = this.xyf2nest(x, y, f); - } - else { - result[i] = -1; - } - } + result[i] = -1n; } - return result; + } } - ; - nside2order(nside) { - return ((nside & (nside - 1)) != 0) ? -1 : Math.log2(nside); + return result; + } + ; + nside2order(nside) { + return ((nside & (nside - 1)) != 0) ? -1 : Math.log2(nside); + } + ; + // x = (...x2 x1 x0)_2 <- in binary + // y = (...y2 y1 y0)_2 + // p = (...y2 x2 y1 x1 y0 x0)_2 + // returns x, y + bit_decombine(p) { + p = BigInt(p); + let x = 0n; + let y = 0n; + let shift = 0n; + while (p > 0n) { + x |= (p & 1n) << shift; + y |= ((p >> 1n) & 1n) << shift; + p >>= 2n; + shift += 1n; } - ; - nest2xyf(ipix) { - ipix = BigInt(ipix); - let pix = ipix & (this.npface - 1n); - let xyf = new Xyf(Number(this.compress_bits(pix)), Number(this.compress_bits(pix >> 1n)), Number(ipix >> (2n * BigInt(this.order)))); - return xyf; + return { x, y }; + } + + // f: base pixel index + // x: north east index in base pixel + // y: north west index in base pixel + nest2fxy(ipix) { + ipix = BigInt(ipix); + let nside2 = BigInt(this.nside) * BigInt(this.nside); + let f = ipix / nside2; // base pixel index + let k = ipix % nside2; // nested pixel index in base pixel + let { x, y } = this.bit_decombine(k); + return { f: f, x: x, y: y }; + } + + fxy2ring(f, x, y) { + let bnside = BigInt(this.nside); + f = BigInt(f); + x = BigInt(x); + y = BigInt(y); + + let f_row = f / 4n; // {0 .. 2} + let f1 = f_row + 2n; // {2 .. 4} + let v = x + y; + let i = f1 * bnside - v - 1n; + + if (i < bnside) { // north polar cap + let f_col = f % 4n; + let ipix = 2n * i * (i - 1n) + (i * f_col) + bnside - y - 1n; + return ipix; } - ; - xyf2nest(ix, iy, face_num) { - return (BigInt(face_num) << (2n * BigInt(this.order))) + + if (i < 3n * bnside) { // equatorial belt + let h = x - y; + let f2 = 2n * (f % 4n) - (f_row % 2n) + 1n; // {0 .. 7} + let k = (f2 * bnside + h + (8n * bnside)) % (8n * bnside); + let offset = 2n * bnside * (bnside - 1n); + let ipix = offset + (i - bnside) * 4n * bnside + (k >> 1n); + return ipix; + } else { // south polar cap + let i_i = 4n * bnside - i; + let i_f_col = 3n - (f % 4n); + let j = 4n * i_i - (i_i * i_f_col) - y; + let i_j = 4n * i_i - j + 1n; + let ipix = 12n * bnside * bnside - 2n * i_i * (i_i - 1n) - i_j; + return ipix; + } + } + + // Convert HEALPix NESTED index to RING index + nest2ring(ipix) { + let { f, x, y } = this.nest2fxy(ipix); + return this.fxy2ring(f, x, y); + } + + nest2xyf(ipix) { + ipix = BigInt(ipix); + let pix = ipix & (this.npface - 1n); + let xyf = new Xyf(Number(this.compress_bits(pix)), Number(this.compress_bits(pix >> 1n)), Number(ipix >> (2n * BigInt(this.order)))); + return xyf; + } + ; + xyf2nest(ix, iy, face_num) { + return (BigInt(face_num) << (2n * BigInt(this.order))) + this.spread_bits(BigInt(ix)) + (this.spread_bits(BigInt(iy)) << 1n); + } + ; + loc2pix(hploc) { + let z = hploc.z; + let phi = hploc.phi; + let za = Math.abs(z); + let tt = this.fmodulo((phi * this.inv_halfpi), 4.0); // in [0,4) + let pixNo; + if (za <= this.twothird) { // Equatorial region + let temp1 = this.nside * (0.5 + tt); + let temp2 = this.nside * (z * 0.75); + let jp = Math.floor(temp1 - temp2); // index of ascending edge line + let jm = Math.floor(temp1 + temp2); // index of descending edge line + let ifp = Math.floor(jp >>> this.order); // in {0,4} + let ifm = Math.floor(jm >>> this.order); + let face_num = Math.floor((ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8))); + let ix = Math.floor(jm & (this.nside - 1)); + let iy = Math.floor(this.nside - (jp & (this.nside - 1)) - 1); + pixNo = this.xyf2nest(ix, iy, face_num); } - ; - loc2pix(hploc) { - let z = hploc.z; - let phi = hploc.phi; - let za = Math.abs(z); - let tt = this.fmodulo((phi * this.inv_halfpi), 4.0); // in [0,4) - let pixNo; - if (za <= this.twothird) { // Equatorial region - let temp1 = this.nside * (0.5 + tt); - let temp2 = this.nside * (z * 0.75); - let jp = Math.floor(temp1 - temp2); // index of ascending edge line - let jm = Math.floor(temp1 + temp2); // index of descending edge line - let ifp = Math.floor(jp >>> this.order); // in {0,4} - let ifm = Math.floor(jm >>> this.order); - let face_num = Math.floor((ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8))); - let ix = Math.floor(jm & (this.nside - 1)); - let iy = Math.floor(this.nside - (jp & (this.nside - 1)) - 1); - pixNo = this.xyf2nest(ix, iy, face_num); - } - else { // polar region, za > 2/3 - let ntt = Math.min(3, Math.floor(tt)); - let tp = tt - ntt; - let tmp = ((za < 0.99) || (!hploc.have_sth)) ? - this.nside * Math.sqrt(3 * (1 - za)) : - this.nside * hploc.sth / Math.sqrt((1.0 + za) / 3.); - let jp = Math.floor(tp * tmp); // increasing edge line index - let jm = Math.floor((1.0 - tp) * tmp); // decreasing edge line index - if (jp >= this.nside) { - jp = this.nside - 1; // for points too close to the boundary - } - if (jm >= this.nside) { - jm = this.nside - 1; - } - if (z >= 0) { - pixNo = this.xyf2nest(Math.floor(this.nside - jm - 1), Math.floor(this.nside - jp - 1), ntt); - } - else { - pixNo = this.xyf2nest(Math.floor(jp), Math.floor(jm), ntt + 8); - } - } - return pixNo; + else { // polar region, za > 2/3 + let ntt = Math.min(3, Math.floor(tt)); + let tp = tt - ntt; + let tmp = ((za < 0.99) || (!hploc.have_sth)) ? + this.nside * Math.sqrt(3 * (1 - za)) : + this.nside * hploc.sth / Math.sqrt((1.0 + za) / 3.); + let jp = Math.floor(tp * tmp); // increasing edge line index + let jm = Math.floor((1.0 - tp) * tmp); // decreasing edge line index + if (jp >= this.nside) { + jp = this.nside - 1; // for points too close to the boundary + } + if (jm >= this.nside) { + jm = this.nside - 1; + } + if (z >= 0) { + pixNo = this.xyf2nest(Math.floor(this.nside - jm - 1), Math.floor(this.nside - jp - 1), ntt); + } + else { + pixNo = this.xyf2nest(Math.floor(jp), Math.floor(jm), ntt + 8); + } } - ; - /** Returns the normalized 3-vector corresponding to the center of the + return pixNo; + } + ; + /** Returns the normalized 3-vector corresponding to the center of the supplied pixel. @param pix long the requested pixel number. @return the pixel's center coordinates. */ - pix2vec(pix) { - return this.pix2loc(pix).toVec3(); - } - ; - /** Returns the Zphi corresponding to the center of the supplied pixel. + pix2vec(pix) { + return this.pix2loc(pix).toVec3(); + } + ; + /** Returns the Zphi corresponding to the center of the supplied pixel. @param pix the requested pixel number. @return the pixel's center coordinates. */ - pix2zphi(pix) { - return this.pix2loc(pix).toZphi(); - } - pix2ang(pix, mirror) { - return this.pix2loc(pix).toPointing(mirror); - } - /** + pix2zphi(pix) { + return this.pix2loc(pix).toZphi(); + } + pix2ang(pix, mirror) { + return this.pix2loc(pix).toPointing(mirror); + } + /** * @param pix long * @return Hploc */ - pix2loc(pix) { - let loc = new Hploc(undefined); - let xyf = this.nest2xyf(pix); - let jr = ((this.jrll[xyf.face]) << this.order) - xyf.ix - xyf.iy - 1; - let nr; - if (jr < this.nside) { - nr = jr; - let tmp = (nr * nr) * this.fact2; - loc.z = 1 - tmp; - if (loc.z > 0.99) { - loc.sth = Math.sqrt(tmp * (2. - tmp)); - loc.have_sth = true; - } - } - else if (jr > this.nl3) { - nr = this.nl4 - jr; - let tmp = (nr * nr) * this.fact2; - loc.z = tmp - 1; - if (loc.z < -0.99) { - loc.sth = Math.sqrt(tmp * (2. - tmp)); - loc.have_sth = true; - } - } - else { - nr = this.nside; - loc.z = (this.nl2 - jr) * this.fact1; - } - let tmp = (this.jpll[xyf.face]) * nr + xyf.ix - xyf.iy; - // assert(tmp<8*nr); // must not happen - if (tmp < 0) { - tmp += 8 * nr; - } - loc.phi = (nr == this.nside) ? 0.75 * Constants.halfpi * tmp * this.fact1 : (0.5 * Constants.halfpi * tmp) / nr; - // loc.setPhi((nr == this.nside) ? 0.75 * Constants.halfpi * tmp * this.fact1 : (0.5 * Constants.halfpi * tmp)/nr); - return loc; + pix2loc(pix) { + let loc = new Hploc(undefined); + let xyf = this.nest2xyf(pix); + let jr = (this.jrll[xyf.face] * Math.pow(2, this.order)) - xyf.ix - xyf.iy - 1; + let nr; + if (jr < this.nside) { + nr = jr; + let tmp = (nr * nr) * this.fact2; + loc.z = 1 - tmp; + if (loc.z > 0.99) { + loc.sth = Math.sqrt(tmp * (2. - tmp)); + loc.have_sth = true; + } } - ; - za2vec(z, a) { - const sin_theta = Math.sqrt(1 - z * z); - const X = sin_theta * Math.cos(a); - const Y = sin_theta * Math.sin(a); - return new Vec3(X, Y, z); + else if (jr > this.nl3) { + nr = this.nl4 - jr; + let tmp = (nr * nr) * this.fact2; + loc.z = tmp - 1; + if (loc.z < -0.99) { + loc.sth = Math.sqrt(tmp * (2. - tmp)); + loc.have_sth = true; + } } - ang2vec(theta, phi) { - const z = Math.cos(theta); - return this.za2vec(z, phi); + else { + nr = this.nside; + loc.z = (this.nl2 - jr) * this.fact1; } - vec2ang(v) { - const { z, a } = this.vec2za(v.getX(), v.getY(), v.getZ()); - return { theta: Math.acos(z), phi: a }; + let tmp = (this.jpll[xyf.face]) * nr + xyf.ix - xyf.iy; + // assert(tmp<8*nr); // must not happen + if (tmp < 0) { + tmp += 8 * nr; } - vec2za(X, Y, z) { - const r2 = X * X + Y * Y; - if (r2 == 0) - return { z: z < 0 ? -1 : 1, a: 0 }; - else { - const PI2 = Math.PI / 2; - const a = (Math.atan2(Y, X) + PI2) % PI2; - z /= Math.sqrt(z * z + r2); - return { z, a }; - } + loc.phi = (nr == this.nside) ? 0.75 * Constants.halfpi * tmp * this.fact1 : (0.5 * Constants.halfpi * tmp) / nr; + // loc.setPhi((nr == this.nside) ? 0.75 * Constants.halfpi * tmp * this.fact1 : (0.5 * Constants.halfpi * tmp)/nr); + return loc; + } + ; + za2vec(z, a) { + const sin_theta = Math.sqrt(1 - z * z); + const X = sin_theta * Math.cos(a); + const Y = sin_theta * Math.sin(a); + return new Vec3(X, Y, z); + } + ang2vec(theta, phi) { + const z = Math.cos(theta); + return this.za2vec(z, phi); + } + vec2ang(v) { + const { z, a } = this.vec2za(v.getX(), v.getY(), v.getZ()); + return { theta: Math.acos(z), phi: a }; + } + vec2za(X, Y, z) { + const r2 = X * X + Y * Y; + if (r2 == 0) + return { z: z < 0 ? -1 : 1, a: 0 }; + else { + const PI2 = Math.PI / 2; + const a = (Math.atan2(Y, X) + PI2) % PI2; + z /= Math.sqrt(z * z + r2); + return { z, a }; } - ang2pix(ptg, mirror) { - return this.loc2pix(new Hploc(ptg)); + } + ang2pix(ptg, mirror) { + return this.loc2pix(new Hploc(ptg)); + } + ; + fmodulo(v1, v2) { + if (v1 >= 0) { + return (v1 < v2) ? v1 : v1 % v2; } - ; - fmodulo(v1, v2) { - if (v1 >= 0) { - return (v1 < v2) ? v1 : v1 % v2; - } - var tmp = v1 % v2 + v2; - return (tmp === v2) ? 0.0 : tmp; + var tmp = v1 % v2 + v2; + return (tmp === v2) ? 0.0 : tmp; + } + ; + compress_bits(v) { + v = BigInt(v); + let compressed = 0n; + for (let i = 0n; i < 32n; i++) { + if ((v & (1n << (2n * i))) !== 0n) { + compressed |= (1n << i); + } } - ; - compress_bits(v) { - v = BigInt(v); - let raw = (v & 0x5555n) | ((v & 0x55550000n) >> 15n); - let compressed = BigInt(this.ctab[Number(raw & 0xffn)]) | (BigInt(this.ctab[Number(raw >> 8n)]) << 4n); - return compressed; + return compressed; + } + ; + spread_bits(v) { + v = BigInt(v); + let spread = 0n; + for (let i = 0n; i < 32n; i++) { + if ((v & (1n << i)) !== 0n) { + spread |= (1n << (2n * i)); + } } - ; - spread_bits(v) { - v = BigInt(v); - return BigInt(this.utab[Number(v & 0xffn)]) - | (BigInt(this.utab[Number((v >> 8n) & 0xffn)]) << 16n) - | (BigInt(this.utab[Number((v >> 16n) & 0xffn)]) << 32n) - | (BigInt(this.utab[Number((v >> 24n) & 0xffn)]) << 48n); - } - ; - /** + return spread; + } + ; + /** * Returns a range set of pixels that overlap with the convex polygon * defined by the {@code vertex} array. *

@@ -470,79 +547,79 @@ export class Healpix { * choice would be 4. * @return the requested set of pixel number ranges */ - queryPolygonInclusive(vertex, fact) { - let inclusive = (fact != 0); - let nv = vertex.length; - // let ncirc = inclusive ? nv+1 : nv; - if (!(nv >= 3)) { - console.log("not enough vertices in polygon"); - return; - } - let vv = new Array(); - for (let i = 0; i < nv; ++i) { - vv[i] = Vec3.pointing2Vec3(vertex[i]); - } - let normal = new Array(); - let flip = 0; - let index = 0; - let back = false; - while (index < vv.length) { - let first = vv[index]; - let medium = null; - let last = null; - if (index == vv.length - 1) { - last = vv[1]; - medium = vv[0]; - } - else if (index == vv.length - 2) { - last = vv[0]; - medium = vv[index + 1]; - } - else { - medium = vv[index + 1]; - last = vv[index + 2]; - } - normal[index] = first.cross(medium).norm(); - let hnd = normal[index].dot(last); - if (index == 0) { - flip = (hnd < 0.) ? -1 : 1; - let tmp = new Pointing(first); // TODO not used - back = false; - } - else { - let flipThnd = flip * hnd; - if (flipThnd < 0) { - let tmp = new Pointing(medium); - vv.splice(index + 1, 1); - normal.splice(index, 1); - back = true; - index -= 1; - continue; - } - else { - let tmp = new Pointing(first); - back = false; - } - } - normal[index].scale(flip); - index += 1; - } - nv = vv.length; - let ncirc = inclusive ? nv + 1 : nv; - let rad = new Array(ncirc); - rad = rad.fill(Constants.halfpi); - // rad = rad.fill(1.5707963267948966); - // let p = "1.5707963267948966"; - // rad = rad.fill(parseFloat(p)); - if (inclusive) { - let cf = new CircleFinder(vv); - normal[nv] = cf.getCenter(); - rad[nv] = Hploc.acos(cf.getCosrad()); - } - return this.queryMultiDisc(normal, rad, fact); + queryPolygonInclusive(vertex, fact) { + let inclusive = (fact != 0); + let nv = vertex.length; + // let ncirc = inclusive ? nv+1 : nv; + if (!(nv >= 3)) { + console.log("not enough vertices in polygon"); + return; } - ; - /** + let vv = new Array(); + for (let i = 0; i < nv; ++i) { + vv[i] = Vec3.pointing2Vec3(vertex[i]); + } + let normal = new Array(); + let flip = 0; + let index = 0; + let back = false; + while (index < vv.length) { + let first = vv[index]; + let medium = null; + let last = null; + if (index == vv.length - 1) { + last = vv[1]; + medium = vv[0]; + } + else if (index == vv.length - 2) { + last = vv[0]; + medium = vv[index + 1]; + } + else { + medium = vv[index + 1]; + last = vv[index + 2]; + } + normal[index] = first.cross(medium).norm(); + let hnd = normal[index].dot(last); + if (index == 0) { + flip = (hnd < 0.) ? -1 : 1; + let tmp = new Pointing(first); // TODO not used + back = false; + } + else { + let flipThnd = flip * hnd; + if (flipThnd < 0) { + let tmp = new Pointing(medium); + vv.splice(index + 1, 1); + normal.splice(index, 1); + back = true; + index -= 1; + continue; + } + else { + let tmp = new Pointing(first); + back = false; + } + } + normal[index].scale(flip); + index += 1; + } + nv = vv.length; + let ncirc = inclusive ? nv + 1 : nv; + let rad = new Array(ncirc); + rad = rad.fill(Constants.halfpi); + // rad = rad.fill(1.5707963267948966); + // let p = "1.5707963267948966"; + // rad = rad.fill(parseFloat(p)); + if (inclusive) { + let cf = new CircleFinder(vv); + normal[nv] = cf.getCenter(); + rad[nv] = Hploc.acos(cf.getCosrad()); + } + return this.queryMultiDisc(normal, rad, fact); + } + ; + /** * For NEST schema only * * @param normal: @@ -556,85 +633,85 @@ export class Healpix { * choice would be 4. * @return RangeSet the requested set of pixel number ranges */ - queryMultiDisc(norm, rad, fact) { - this.computeBn(); - let inclusive = (fact != 0); - let nv = norm.length; - // HealpixUtils.check(nv==rad.lengt0,"inconsistent input arrays"); - if (!(nv == rad.length)) { - console.error("inconsistent input arrays"); - return; - } - let res = new RangeSet(4 << 1); - // Removed code for Scheme.RING - let oplus = 0; - if (inclusive) { - if (!(Math.pow(2, this.order_max - this.order) >= fact)) { - console.error("invalid oversampling factor"); - } - if (!((fact & (fact - 1)) == 0)) { - console.error("oversampling factor must be a power of 2"); - } - oplus = this.ilog2(fact); - } - let omax = this.order + oplus; // the order up to which we test - // TODO: ignore all disks with radius>=pi - // let crlimit = new Float32Array[omax+1][nv][3]; - let crlimit = new Array(omax + 1); - let o; - let i; - for (o = 0; o <= omax; ++o) { // prepare data at the required orders - crlimit[o] = new Array(nv); - let dr = this.bn[o].maxPixrad(); // safety distance - for (i = 0; i < nv; ++i) { - crlimit[o][i] = new Float64Array(3); - crlimit[o][i][0] = (rad[i] + dr > Math.PI) ? -1 : Hploc.cos(rad[i] + dr); - crlimit[o][i][1] = (o == 0) ? Hploc.cos(rad[i]) : crlimit[0][i][1]; - crlimit[o][i][2] = (rad[i] - dr < 0.) ? 1. : Hploc.cos(rad[i] - dr); - } - } - let stk = new pstack(12 + 3 * omax); - for (let i = 0n; i < 12n; i++) { // insert the 12 base pixels in reverse - // order - stk.push(11n - i, 0); - } - while (stk.size() > 0) { // as long as there are pixels on the stack - // pop current pixel number and order from the stack - let pix = stk.ptop(); - let o = stk.otop(); - stk.pop(); - let pv = this.bn[o].pix2vec(pix); - let zone = 3; - for (let i = 0; (i < nv) && (zone > 0); ++i) { - let crad = pv.dot(norm[i]); - for (let iz = 0; iz < zone; ++iz) { - if (crad < crlimit[o][i][iz]) { - zone = iz; - } - } - } - if (zone > 0) { - this.check_pixel(o, omax, zone, res, pix, stk, inclusive); - } - } - return res; + queryMultiDisc(norm, rad, fact) { + this.computeBn(); + let inclusive = (fact != 0); + let nv = norm.length; + // HealpixUtils.check(nv==rad.lengt0,"inconsistent input arrays"); + if (!(nv == rad.length)) { + console.error("inconsistent input arrays"); + return; } - ; - /** Integer base 2 logarithm. + let res = new RangeSet(4 << 1); + // Removed code for Scheme.RING + let oplus = 0; + if (inclusive) { + if (!(Math.pow(2, this.order_max - this.order) >= fact)) { + console.error("invalid oversampling factor"); + } + if (!((fact & (fact - 1)) == 0)) { + console.error("oversampling factor must be a power of 2"); + } + oplus = this.ilog2(fact); + } + let omax = this.order + oplus; // the order up to which we test + // TODO: ignore all disks with radius>=pi + // let crlimit = new Float32Array[omax+1][nv][3]; + let crlimit = new Array(omax + 1); + let o; + let i; + for (o = 0; o <= omax; ++o) { // prepare data at the required orders + crlimit[o] = new Array(nv); + let dr = this.bn[o].maxPixrad(); // safety distance + for (i = 0; i < nv; ++i) { + crlimit[o][i] = new Float64Array(3); + crlimit[o][i][0] = (rad[i] + dr > Math.PI) ? -1 : Hploc.cos(rad[i] + dr); + crlimit[o][i][1] = (o == 0) ? Hploc.cos(rad[i]) : crlimit[0][i][1]; + crlimit[o][i][2] = (rad[i] - dr < 0.) ? 1. : Hploc.cos(rad[i] - dr); + } + } + let stk = new pstack(12 + 3 * omax); + for (let i = 0n; i < 12n; i++) { // insert the 12 base pixels in reverse + // order + stk.push(11n - i, 0); + } + while (stk.size() > 0) { // as long as there are pixels on the stack + // pop current pixel number and order from the stack + let pix = stk.ptop(); + let o = stk.otop(); + stk.pop(); + let pv = this.bn[o].pix2vec(pix); + let zone = 3; + for (let i = 0; (i < nv) && (zone > 0); ++i) { + let crad = pv.dot(norm[i]); + for (let iz = 0; iz < zone; ++iz) { + if (crad < crlimit[o][i][iz]) { + zone = iz; + } + } + } + if (zone > 0) { + this.check_pixel(o, omax, zone, res, pix, stk, inclusive); + } + } + return res; + } + ; + /** Integer base 2 logarithm. @param arg @return the largest integer {@code n} that fulfills {@code 2^n<=arg}. For negative arguments and zero, 0 is returned. */ - ilog2(arg) { - let max = Math.max(arg, 1); - return 31 - Math.clz32(max); - } - ; - /** Computes the cosine of the angular distance between two z, phi positions + ilog2(arg) { + let max = Math.max(arg, 1); + return 31 - Math.clz32(max); + } + ; + /** Computes the cosine of the angular distance between two z, phi positions on the unit sphere. */ - cosdist_zphi(z1, phi1, z2, phi2) { - return z1 * z2 + Hploc.cos(phi1 - phi2) * Math.sqrt((1.0 - z1 * z1) * (1.0 - z2 * z2)); - } - /** + cosdist_zphi(z1, phi1, z2, phi2) { + return z1 * z2 + Hploc.cos(phi1 - phi2) * Math.sqrt((1.0 - z1 * z1) * (1.0 - z2 * z2)); + } + /** * @param int o * @param int omax * @param int zone @@ -643,82 +720,82 @@ export class Healpix { * @param pstack stk * @param boolean inclusive */ - check_pixel(o, omax, zone, pixset, pix, stk, inclusive) { - if (zone == 0) - return; - if (o < this.order) { - if (zone >= 3) { // output all subpixels - let sdist = 2n * BigInt(this.order - o); // the "bit-shift distance" between map orders - pixset.append1(pix << sdist, ((pix + 1n) << sdist)); - } - else { // (zone>=1) - for (let i = 0n; i < 4n; ++i) { - stk.push(4n * pix + 3n - i, o + 1); // add children - } - } - } - else if (o > this.order) { // this implies that inclusive==true - if (zone >= 2) { // pixel center in shape - pixset.append(pix >> (2n * BigInt(o - this.order))); // output the parent pixel at order - stk.popToMark(); // unwind the stack - } - else { // (zone>=1): pixel center in safety range - if (o < omax) { // check sublevels - for (let i = 0n; i < 4n; ++i) { // add children in reverse order - stk.push(4n * pix + 3n - i, o + 1); // add children - } - } - else { // at resolution limit - pixset.append(pix >> (2n * BigInt(o - this.order))); // output the parent pixel at order - stk.popToMark(); // unwind the stack - } - } - } - else { // o==order - if (zone >= 2) { - pixset.append(pix); - } - else if (inclusive) { // and (zone>=1) - if (this.order < omax) { // check sublevels - stk.mark(); // remember current stack position - for (let i = 0n; i < 4n; ++i) { // add children in reverse order - stk.push(4n * pix + 3n - i, o + 1); // add children - } - } - else { // at resolution limit - pixset.append(pix); // output the pixel - } - } + check_pixel(o, omax, zone, pixset, pix, stk, inclusive) { + if (zone == 0) + return; + if (o < this.order) { + if (zone >= 3) { // output all subpixels + let sdist = 2n * BigInt(this.order - o); // the "bit-shift distance" between map orders + pixset.append1(pix << sdist, ((pix + 1n) << sdist)); + } + else { // (zone>=1) + for (let i = 0n; i < 4n; ++i) { + stk.push(4n * pix + 3n - i, o + 1); // add children } + } } - /** Returns the maximum angular distance between a pixel center and its + else if (o > this.order) { // this implies that inclusive==true + if (zone >= 2) { // pixel center in shape + pixset.append(pix >> (2n * BigInt(o - this.order))); // output the parent pixel at order + stk.popToMark(); // unwind the stack + } + else { // (zone>=1): pixel center in safety range + if (o < omax) { // check sublevels + for (let i = 0n; i < 4n; ++i) { // add children in reverse order + stk.push(4n * pix + 3n - i, o + 1); // add children + } + } + else { // at resolution limit + pixset.append(pix >> (2n * BigInt(o - this.order))); // output the parent pixel at order + stk.popToMark(); // unwind the stack + } + } + } + else { // o==order + if (zone >= 2) { + pixset.append(pix); + } + else if (inclusive) { // and (zone>=1) + if (this.order < omax) { // check sublevels + stk.mark(); // remember current stack position + for (let i = 0n; i < 4n; ++i) { // add children in reverse order + stk.push(4n * pix + 3n - i, o + 1); // add children + } + } + else { // at resolution limit + pixset.append(pix); // output the pixel + } + } + } + } + /** Returns the maximum angular distance between a pixel center and its corners. @return maximum angular distance between a pixel center and its corners. */ - maxPixrad() { - let zphia = new Zphi(2. / 3., Math.PI / this.nl4); - let xyz1 = this.convertZphi2xyz(zphia); - let va = new Vec3(xyz1[0], xyz1[1], xyz1[2]); - let t1 = 1. - 1. / this.nside; - t1 *= t1; - let zphib = new Zphi(1 - t1 / 3, 0); - let xyz2 = this.convertZphi2xyz(zphib); - let vb = new Vec3(xyz2[0], xyz2[1], xyz2[2]); - return va.angle(vb); - } - ; - /** + maxPixrad() { + let zphia = new Zphi(2. / 3., Math.PI / this.nl4); + let xyz1 = this.convertZphi2xyz(zphia); + let va = new Vec3(xyz1[0], xyz1[1], xyz1[2]); + let t1 = 1. - 1. / this.nside; + t1 *= t1; + let zphib = new Zphi(1 - t1 / 3, 0); + let xyz2 = this.convertZphi2xyz(zphib); + let vb = new Vec3(xyz2[0], xyz2[1], xyz2[2]); + return va.angle(vb); + } + ; + /** * this is a workaround replacing the Vec3(Zphi) constructor. */ - convertZphi2xyz(zphi) { - let sth = Math.sqrt((1.0 - zphi.z) * (1.0 + zphi.z)); - let x = sth * Hploc.cos(zphi.phi); - let y = sth * Hploc.sin(zphi.phi); - let z = zphi.z; - return [x, y, z]; - } - ; - /** Returns a range set of pixels which overlap with a given disk.

+ convertZphi2xyz(zphi) { + let sth = Math.sqrt((1.0 - zphi.z) * (1.0 + zphi.z)); + let x = sth * Hploc.cos(zphi.phi); + let y = sth * Hploc.sin(zphi.phi); + let z = zphi.z; + return [x, y, z]; + } + ; + /** Returns a range set of pixels which overlap with a given disk.

This method is more efficient in the RING scheme.

This method may return some pixels which don't overlap with the polygon at all. The higher {@code fact} is chosen, the fewer false @@ -729,53 +806,120 @@ export class Healpix { {@code fact*nside}. For NESTED ordering, {@code fact} must be a power of 2, else it can be any positive integer. A typical choice would be 4. @return the requested set of pixel number ranges */ - queryDiscInclusive(ptg, radius, fact) { - this.computeBn(); - let inclusive = (fact != 0); - let pixset = new RangeSet(); - if (radius >= Math.PI) { // disk covers the whole sphere - pixset.append1(0n, this.npix); - return pixset; - } - let oplus = 0; - if (inclusive) { - // HealpixUtils.check ((1L<=fact,"invalid oversampling factor"); - if (!((fact & (fact - 1)) == 0)) { - console.error("oversampling factor must be a power of 2"); - } - oplus = this.ilog2(fact); - } - let omax = Math.min(this.order_max, this.order + oplus); // the order up to which we test - let vptg = Vec3.pointing2Vec3(ptg); - let crpdr = new Array(omax + 1); - let crmdr = new Array(omax + 1); - let cosrad = Hploc.cos(radius); - let sinrad = Hploc.sin(radius); - for (let o = 0; o <= omax; o++) { // prepare data at the required orders - let dr = this.mpr[o]; // safety distance - let cdr = this.cmpr[o]; - let sdr = this.smpr[o]; - crpdr[o] = (radius + dr > Math.PI) ? -1. : cosrad * cdr - sinrad * sdr; - crmdr[o] = (radius - dr < 0.) ? 1. : cosrad * cdr + sinrad * sdr; - } - let stk = new pstack(12 + 3 * omax); - for (let i = 0n; i < 12n; i++) { // insert the 12 base pixels in reverse order - stk.push(11n - i, 0); - } - while (stk.size() > 0) { // as long as there are pixels on the stack - // pop current pixel number and order from the stack - let pix = stk.ptop(); - let curro = stk.otop(); - stk.pop(); - let pos = this.bn[curro].pix2zphi(pix); - // cosine of angular distance between pixel center and disk center - let cangdist = this.cosdist_zphi(vptg.z, ptg.phi, pos.z, pos.phi); - if (cangdist > crpdr[curro]) { - let zone = (cangdist < cosrad) ? 1 : ((cangdist <= crmdr[curro]) ? 2 : 3); - this.check_pixel(curro, omax, zone, pixset, pix, stk, inclusive); - } - } - return pixset; + queryDiscInclusive(ptg, radius, fact) { + this.computeBn(); + let inclusive = (fact != 0); + let pixset = new RangeSet(); + if (radius >= Math.PI) { // disk covers the whole sphere + pixset.append1(0n, this.npix); + return pixset; } + let oplus = 0; + if (inclusive) { + // HealpixUtils.check ((1L<=fact,"invalid oversampling factor"); + if (!((fact & (fact - 1)) == 0)) { + console.error("oversampling factor must be a power of 2"); + } + oplus = this.ilog2(fact); + } + let omax = Math.min(this.order_max, this.order + oplus); // the order up to which we test + let vptg = Vec3.pointing2Vec3(ptg); + let crpdr = new Array(omax + 1); + let crmdr = new Array(omax + 1); + let cosrad = Hploc.cos(radius); + let sinrad = Hploc.sin(radius); + for (let o = 0; o <= omax; o++) { // prepare data at the required orders + let dr = this.mpr[o]; // safety distance + let cdr = this.cmpr[o]; + let sdr = this.smpr[o]; + crpdr[o] = (radius + dr > Math.PI) ? -1. : cosrad * cdr - sinrad * sdr; + crmdr[o] = (radius - dr < 0.) ? 1. : cosrad * cdr + sinrad * sdr; + } + let stk = new pstack(12 + 3 * omax); + for (let i = 0n; i < 12n; i++) { // insert the 12 base pixels in reverse order + stk.push(11n - i, 0); + } + while (stk.size() > 0) { // as long as there are pixels on the stack + // pop current pixel number and order from the stack + let pix = stk.ptop(); + let curro = stk.otop(); + stk.pop(); + let pos = this.bn[curro].pix2zphi(pix); + // cosine of angular distance between pixel center and disk center + let cangdist = this.cosdist_zphi(vptg.z, ptg.phi, pos.z, pos.phi); + if (cangdist > crpdr[curro]) { + let zone = (cangdist < cosrad) ? 1 : ((cangdist <= crmdr[curro]) ? 2 : 3); + this.check_pixel(curro, omax, zone, pixset, pix, stk, inclusive); + } + } + return pixset; + } + + clipPolygonAgainstZPlane(points) { + if (points.length === 0) return []; + + const clipped = []; + + for (let i = 0; i < points.length; i++) { + const currPt = points[i]; + const prevPt = points[i === 0 ? points.length - 1 : i - 1]; // Wraparound + + const currInside = currPt.z >= 0; + const prevInside = prevPt.z >= 0; + + if (currInside !== prevInside) { + const t = prevPt.z / (prevPt.z - currPt.z); + const intersectPt = { + x: prevPt.x + t * (currPt.x - prevPt.x), + y: prevPt.y + t * (currPt.y - prevPt.y), + z: 0 + }; + + const len = Math.sqrt(intersectPt.x * intersectPt.x + intersectPt.y * intersectPt.y); + if (len > 0) { + intersectPt.x /= len; + intersectPt.y /= len; + } + + intersectPt.isExit = prevInside; + intersectPt.isEntry = currInside; + + clipped.push(intersectPt); + } + + if (currInside) { + clipped.push({ x: currPt.x, y: currPt.y, z: currPt.z }); + } + } + + const finalClipped = []; + for (let i = 0; i < clipped.length; i++) { + const pt = clipped[i]; + finalClipped.push({ x: pt.x, y: pt.y, z: pt.z }); + + if (pt.isExit) { + const nextIdx = (i + 1) % clipped.length; + const nextPt = clipped[nextIdx]; + + if (nextPt.isEntry) { + const theta1 = Math.atan2(pt.y, pt.x); + const theta2 = Math.atan2(nextPt.y, nextPt.x); + let dTheta = theta2 - theta1; + while (dTheta > Math.PI) dTheta -= 2 * Math.PI; + while (dTheta < -Math.PI) dTheta += 2 * Math.PI; + + const segments = Math.max(2, Math.ceil(Math.abs(dTheta) * 32 / Math.PI)); + + for (let step = 1; step < segments; step++) { + const t = step / segments; + const theta = theta1 + dTheta * t; + finalClipped.push({ x: Math.cos(theta), y: Math.sin(theta), z: 0 }); + } + } + } + } + + return finalClipped; + } } //# sourceMappingURL=Healpix.js.map \ No newline at end of file