diff --git a/pkg/healpix/healpix.go b/pkg/healpix/healpix.go index b8c621b..1e76ff6 100644 --- a/pkg/healpix/healpix.go +++ b/pkg/healpix/healpix.go @@ -162,6 +162,20 @@ func (h *HealPIX) GetFaceXY(pixel int) (face int, x int, y int) { } /*****************************************************************************************************************/ + +func (h *HealPIX) GetPixelIndexFromFaceXY(face int, x int, y int) int { + switch h.Scheme { + case RING: + return getRingPixelIndexFromFaceXY(face, x, y, h.NSide) + case NESTED: + return getNestedPixelIndexFromFaceXY(face, x, y, h.NSide) + default: + return getRingPixelIndexFromFaceXY(face, x, y, h.NSide) + } +} + +/*****************************************************************************************************************/ + // GetPixelArea returns the area of each pixel in the HEALPix projection, in degrees. func (h *HealPIX) GetPixelArea() float64 { // Get the number of pixels for the given NSide: @@ -645,3 +659,29 @@ func getNestedFaceXY(nside, index int) (face, x, y int) { } /*****************************************************************************************************************/ + +// getRingPixelIndexFromFaceXY converts the given (face, x, y) coordinates to a RING pixel index. +// It does so by first computing the corresponding NESTED pixel index, converting that to spherical +// coordinates, and then converting those spherical coordinates to the RING index. +func getRingPixelIndexFromFaceXY(face, x, y, nside int) int { + // First, get the NESTED pixel index from (face, x, y) + nestedIndex := getNestedPixelIndexFromFaceXY(face, x, y, nside) + // Convert the nested index to spherical coordinates (θ, φ) + theta, phi := convertNestedIndexToSpherical(nside, nestedIndex) + // Convert spherical coordinates to a RING pixel index + return convertSphericalToRingIndex(nside, theta, phi) +} + +/*****************************************************************************************************************/ + +// getNestedPixelIndexFromFaceXY encodes the (face, x, y) into a NESTED pixel index. +// It uses the lookup table `utab` to interleave the bits of x and y. +func getNestedPixelIndexFromFaceXY(face, x, y, nside int) int { + return face*nside*nside + + (utab[x&0xff] | + (utab[(x>>8)&0xff] << 16) | + (utab[y&0xff] << 1) | + (utab[(y>>8)&0xff] << 17)) +} + +/*****************************************************************************************************************/ diff --git a/pkg/healpix/healpix_test.go b/pkg/healpix/healpix_test.go index 851112f..9329430 100644 --- a/pkg/healpix/healpix_test.go +++ b/pkg/healpix/healpix_test.go @@ -1212,3 +1212,23 @@ func TestGetFaceXY(t *testing.T) { } /*****************************************************************************************************************/ + +func TestGetPixelIndexFromFaceXY(t *testing.T) { + nside := 2 + + healpix := NewHealPIX(nside, RING) + + face := 4 + x := 1 + y := 1 + + pixelIndex := healpix.GetPixelIndexFromFaceXY(face, x, y) + + expectedPixelIndex := 12 + + if pixelIndex != expectedPixelIndex { + t.Errorf("Expected Pixel Index=%d, Got Pixel Index=%d", expectedPixelIndex, pixelIndex) + } +} + +/*****************************************************************************************************************/