#!/usr/bin/env python3 """Lucky-imaging / stacking demonstration with a measurable resolution ceiling. Ground truth is a synthetic phantom, not a photograph: a Jupiter-like banded field beside a line-pair resolution ladder at known angular spacings. Nothing is imported that could smuggle real detail in. Each frame is the phantom convolved with (diffraction PSF for the aperture) x (a randomly drawn seeing PSF), displaced by random tip-tilt, then given Poisson photon noise. "Lucky imaging" keeps the sharpest few per cent, registers them by cross-correlation, and averages. Contrast recovered at each rung is measured as the Fourier amplitude at that rung's own spatial frequency, which is a modulation-transfer measurement rather than a max-minus- min guess. The result to look for: the stack recovers detail toward 1.22*lambda/D and stops there, whatever the frame count. """ import numpy as np from scipy import ndimage RNG = np.random.default_rng(20260827) LAMBDA = 550e-9 PIX = 0.04 # arcsec per pixel H, W = 200, 300 # 8 x 12 arcsec SEEING_MEDIAN = 2.5 # arcsec FWHM NFRAMES, KEEP = 400, 0.05 RL_ITERS = 40 LADDER = [1.6, 0.8, 0.4, 0.2] # arcsec, line-pair period APERTURES = [(100, '100 mm'), (200, '200 mm / 8 in'), (356, '356 mm / 14 in')] PHOTONS = 1600.0 LX0, LX1 = 194, 294 # ladder column range RUNG_H, RUNG_GAP, RUNG_Y0 = 38, 47, 12 def diffraction_fwhm(d_mm): return np.degrees(1.22 * LAMBDA / (d_mm / 1000.0)) * 3600.0 def otf(fwhm_arcsec): sig = max(fwhm_arcsec / 2.3548 / PIX, 1e-6) fy = np.fft.fftfreq(H)[:, None] fx = np.fft.fftfreq(W)[None, :] return np.exp(-2 * (np.pi * sig) ** 2 * (fy ** 2 + fx ** 2)) def rung_rows(i): y0 = RUNG_Y0 + i * RUNG_GAP return y0, y0 + RUNG_H def phantom(): y, x = np.mgrid[:H, :W].astype(float) img = np.full((H, W), 0.50) disc = x < LX0 - 8 for cy, hw, amp in [(38, 9, -0.20), (72, 15, 0.15), (118, 12, -0.24), (162, 8, 0.12)]: img += np.where(disc, amp * np.exp(-((y - cy) ** 2) / (2 * hw ** 2)), 0.0) img += np.where(disc, 0.28 * np.exp(-(((x - 92) / 27.0) ** 2 + ((y - 120) / 14.0) ** 2)), 0.0) for i, sp in enumerate(LADDER): # separated rungs, own rows per = sp / PIX y0, y1 = rung_rows(i) col = np.arange(LX0, LX1) bar = ((col - LX0) % per) < (per / 2.0) img[y0:y1, LX0:LX1] += 0.28 * bar[None, :] return np.clip(img, 0.02, 1.0) def richardson_lucy(img, psf, iters): """Standard RL deconvolution, done in the Fourier domain.""" P = np.fft.fft2(np.fft.ifftshift(psf)) Pc = np.conj(P) est = np.full_like(img, img.mean()) obs = np.clip(img, 1e-6, None) for _ in range(iters): conv = np.real(np.fft.ifft2(np.fft.fft2(est) * P)) est = est * np.real(np.fft.ifft2(np.fft.fft2(obs / np.clip(conv, 1e-6, None)) * Pc)) est = np.clip(est, 0, None) return est def register(ref_f, frame): cc = np.real(np.fft.ifft2(ref_f * np.conj(np.fft.fft2(frame)))) dy, dx = np.unravel_index(np.argmax(cc), cc.shape) if dy > H // 2: dy -= H if dx > W // 2: dx -= W return ndimage.shift(frame, (dy, dx), order=1, mode='nearest') def simulate(truth, d_mm, nframes=NFRAMES): T = np.fft.fft2(truth) D = otf(diffraction_fwhm(d_mm)) frames, sharp, sees = [], [], [] for _ in range(nframes): see = max(0.55, RNG.lognormal(np.log(SEEING_MEDIAN), 0.55)) blur = np.real(np.fft.ifft2(T * D * otf(see))) blur = ndimage.shift(blur, RNG.normal(0, see / PIX / 2.5, 2), order=1, mode='nearest') noisy = RNG.poisson(np.clip(blur, 0, None) * PHOTONS) / PHOTONS frames.append(noisy) sharp.append(np.abs(np.diff(ndimage.gaussian_filter(noisy, 1.0), axis=1)).mean()) sees.append(see) frames = np.array(frames); sees = np.array(sees) best = frames[np.argsort(sharp)[-max(2, int(nframes * KEEP)):]] idx = np.argsort(sharp)[-max(2, int(nframes * KEEP)):] ref_f = np.fft.fft2(best[-1]) stack = np.mean([best[-1]] + [register(ref_f, f) for f in best[:-1]], axis=0) # effective PSF still in the stack: diffraction x the mean retained seeing. # A real imager estimates this by photographing a star; here we use it directly, # which is the best case the pipeline can hope for rather than a typical one. eff = np.real(np.fft.ifft2(D * otf(float(np.mean(sees[idx]))))) eff = np.fft.fftshift(eff); eff = np.clip(eff, 0, None); eff /= eff.sum() return frames[nframes // 2], stack, richardson_lucy(stack, eff, RL_ITERS) def modulation(img): """Fourier amplitude at each rung's own frequency, normalised to its mean.""" out = [] for i, sp in enumerate(LADDER): y0, y1 = rung_rows(i) line = img[y0 + 4:y1 - 4, LX0:LX1].mean(axis=0) n = line.size line = line - line.mean() k = n * PIX / sp # cycles across the strip j = int(round(k)) if j < 1 or j >= n // 2: out.append(0.0); continue amp = np.abs(np.fft.rfft(line * np.hanning(n)))[j] * 4.0 / n out.append(float(amp / img[y0 + 4:y1 - 4, LX0:LX1].mean())) return out if __name__ == '__main__': t = phantom() np.save('/tmp/_truth.npy', t) print(f'phantom {W}x{H} px at {PIX}"/px; ladder periods {LADDER} arcsec') print(f'seeing median {SEEING_MEDIAN}" FWHM; {NFRAMES} frames/aperture, best {KEEP:.0%} kept\n') hdr = ' '.join(f'{s:>6}"' for s in LADDER) print(f'{"aperture":>16s} {"1.22L/D":>9s} {hdr}') print(f'{"(ground truth)":>16s} {"":>9s} ' + ' '.join(f'{v:7.3f}' for v in modulation(t))) for d, name in APERTURES: one, st, dec = simulate(t, d) np.save(f'/tmp/_one_{d}.npy', one); np.save(f'/tmp/_stack_{d}.npy', st) np.save(f'/tmp/_dec_{d}.npy', dec) print(f'{name:>16s} {diffraction_fwhm(d):8.3f}" single ' + ' '.join(f'{v:7.3f}' for v in modulation(one))) print(f'{"":>16s} {"":>9s} stack ' + ' '.join(f'{v:7.3f}' for v in modulation(st))) print(f'{"":>16s} {"":>9s} +decon ' + ' '.join(f'{v:7.3f}' for v in modulation(dec)) + ''.join(' <-- finer than 1.22L/D' if s_ < diffraction_fwhm(d) and i == next( (j for j, q in enumerate(LADDER) if q < diffraction_fwhm(d)), -1) else '' for i, s_ in enumerate(LADDER)))