Add zscale for AAVSO output images

[rzellem]

using this function, stolen from ChatGPT, here:

import numpy as np
import matplotlib.pyplot as plt

def zscale(image, nsamples=1000, contrast=0.25):
"""
Replicates IRAF/DS9 zscale algorithm for image normalization.

Parameters
----------
image : 2D numpy array
    Image array (e.g., from a FITS file).
nsamples : int
    Number of random pixels to sample for scaling (default 1000).
contrast : float
    Contrast parameter (default 0.25, like DS9).

Returns
-------
z1, z2 : float
    Lower and upper limits for display scaling.
"""

# Flatten and sample pixels
samples = image.flatten()
if len(samples) > nsamples:
    samples = np.random.choice(samples, nsamples, replace=False)
samples = np.sort(samples)

n = len(samples)
if n < 10:
    raise ValueError("Not enough pixels for zscale")

# Median
median = np.median(samples)

# Fit line: index vs. value
x = np.arange(n)
fit = np.polyfit(x, samples, 1)
slope = fit[0]

# Adjust slope by contrast
slope /= contrast

# Compute z1, z2
z1 = median - (n/2)*slope
z2 = median + (n/2)*slope

# Clip to min/max
z1 = max(z1, samples[0])
z2 = min(z2, samples[-1])

return z1, z2

----------------- DEMONSTRATION -----------------

if name == "main":
# Make a fake "astronomy-like" image
np.random.seed(42)
background = np.random.normal(1000, 50, (500, 500)) # sky background
stars = np.zeros((500, 500))
for _ in range(50):
x, y = np.random.randint(0, 500, 2)
stars[x, y] = np.random.randint(2000, 5000)
image = background + stars

# Compute zscale
z1, z2 = zscale(image)

# Compare raw vs. zscale-normalized
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].imshow(image, origin='lower', cmap='gray')
axs[0].set_title("Raw Image (unscaled)")
axs[1].imshow(image, origin='lower', cmap='gray', vmin=z1, vmax=z2)
axs[1].set_title(f"Zscale Normalized\nz1={z1:.1f}, z2={z2:.1f}")
plt.show()