Add spherical peak counting functions for HEALPix maps

[Copilot]

This PR implements new functionality for calculating multi-scale peak counts on spherical HEALPix maps, providing a spherical analog to the existing get_wtpeaks function for 2D images.

New Functions Added

get_peaks_sphere(healpix_map, threshold=None, ordered=True, mask=None, nside=None)

Core function for finding local maxima in HEALPix maps using proper spherical geometry:

• Uses healpy.get_all_neighbours() to find pixel neighbors on the sphere
• Handles masking and thresholding for selective peak detection
• Returns peak pixel indices and heights in descending order

get_wtpeaks_sphere(Map, nscales, nbins=None, Mask=None, min_snr=None, max_snr=None, noise_std=None, peak_threshold=None)

Main multi-scale peak counting function that:

• Uses CMRStarlet() for spherical wavelet transform (same as get_wtl1_sphere)
• Identifies peaks at each wavelet scale using spherical neighbor detection
• Returns binned peak count histograms plus detailed peak information
• Supports SNR analysis with noise normalization

Usage Example

import numpy as np
import healpy as hp
from pycs.astro.wl.hos_peaks_l1 import get_wtpeaks_sphere

# Create or load a HEALPix map
nside = 32
healpix_map = hp.read_map("cosmic_map.fits")

# Perform multi-scale peak analysis
bins, counts, pixels_list, heights_list = get_wtpeaks_sphere(
    healpix_map, 
    nscales=5, 
    nbins=20,
    noise_std=0.3  # For SNR analysis
)

# counts[i] contains peak count histogram for scale i
# pixels_list[i] contains peak pixel indices for scale i
# heights_list[i] contains peak amplitudes for scale i

Design Features

• Spherical geometry: Properly handles HEALPix pixel neighbors on the sphere, unlike 2D grid-based approaches
• API consistency: Similar parameter structure to existing functions like get_wtl1_sphere and get_wtpeaks
• Comprehensive analysis: Returns both statistical summaries (histograms) and detailed peak information
• Flexible parameters: Supports masking, thresholding, SNR analysis, and custom binning
• Robust implementation: Handles edge cases including boundary pixels, uniform maps, and missing neighbors

Applications

This functionality enables cosmologists to perform peak-based statistical analysis on spherical maps such as:

• Cosmic microwave background (CMB) temperature and polarization maps
• Weak lensing convergence maps
• Galaxy overdensity maps
• Other cosmological observables mapped to the celestial sphere

The new functions complement the existing get_wtl1_sphere() function by providing peak-based statistics rather than L1-norm based analysis, offering a different perspective on multi-scale structure in spherical data.

Testing

The implementation has been thoroughly tested with:

• Various HEALPix map sizes and configurations
• Different numbers of wavelet scales
• Masking and thresholding scenarios
• SNR-based analysis with noise normalization
• Edge cases including uniform maps and single peaks

All functions integrate seamlessly with the existing CMRStarlet spherical wavelet transform infrastructure.

---

💬 Share your feedback on Copilot coding agent for the chance to win a $200 gift card! Click here to start the survey.