(Re)Introduction to Image Processing

Version 0.1

During Session 1 of the DSFP, Robert Lupton provided a problem that brilliantly introduced some of the basic challenges associated with measuring the flux of a point source. As such, we will revisit that problem as a review/introduction to the remainder of the week.

By AA Miller (CIERA/Northwestern & Adler)
[But please note that this is essentially a copy of Robert's lecture.]



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt

%matplotlib notebook

Problem 1) An (oversimplified) 1-D Model

For this introductory problem we are going to simulate a 1 dimensional detector (the more complex issues associated will real stars on 2D detectors will be covered tomorrow by Dora). We will generate stars as Gaussians $N(\mu, \sigma^2)$, with mean $\mu$ and variance $\sigma^2$.

As observed by LSST, all stars are point sources that reflect the point spread function (PSF), which is produced by a combination of the atmosphere, telescope, and detector. A standard measure of the PSF's width is the Full Width Half Maximum (FWHM).

There is also a smooth background of light from several sources that I previously mentioned (the atmosphere, the detector, etc). We will refer to this background simply as "The Sky".

Problem 1a

Write a function phi() to simulate a (noise-free) 1D Gaussian PSF. The function should take mu and fwhm as arguments, and evaluate the PSF along a user-supplied array x.

Hint - for a Gaussian $N(0, \sigma^2)$, the FWHM is $2\sqrt{2\ln(2)}\,\sigma \approx 2.3548\sigma$.



In [ ]:

    
def phi(x, mu, fwhm):
    # complete

Problem 1b

Plot the noise-free PSF for a star with $\mu = 10$ and $\mathrm{FWHM} = 3$. What is the flux of this star?



In [ ]:

    
x = # complete

plt.plot( # complete
    
print("The flux of the star is: {:.3f}".format( # complete

Problem 1c

Add Sky noise (a constant in this case) to your model. Define the sky as S, with total stellar flux F.

Plot the model for S = 100 and F = 500.



In [ ]:

    
plt.plot( # complete

Problem 2) Add Noise

We will add noise to this simulation assuming that photon counting contributes the only source of uncertainty (this assumption is far from sufficient in real life). Within each pixel, $n$ photons are detected with an uncertainty that follows a Poisson distribution, which has the property that the mean $\mu$ is equal to the variance $\mu$. If $n \gg 1$ then $P(\mu) \approx N(\mu, \mu)$ [you can safely assume we will be in this regime for the remainder of this problem].

Problem 2a

Calculate the noisy flux for the simulated star in Problem 1c.

Hint - you may find the function np.random.normal() helpful.



In [ ]:

    
# complete
noisy_flux = # complete

Problem 2b

Overplot the noisy signal, with the associated uncertainties, on top of the noise-free signal.



In [ ]:

    
plt.plot( # complete
plt.errorbar( # complete

Problem 3) Flux Measurement

We will now attempt to measure the flux from a simulated star.

Problem 3a

Write a function simulate() to simulate the noisy flux measurements of a star with centroid mu, FWHM fwhm, sky background S, and flux F.

Hint - it may be helpful to plot the output of your function.



In [ ]:

    
def simulate(# complete
    # complete
    # complete

Problem 3b

Using an aperture with radius of 5 pixels centered on the source, measure the flux from a star centered at mu = 0, with fwhm = 5, S = 100, and F = 1000.

Hint - assume you can perfectly measure the background, and subtract this prior to the measurement.



In [ ]:

    
# complete
sim_star = simulate( # complete

ap_flux = # complete

print("The star has flux = {:.3f}".format( # complete

Problem 3c

Write a Monte Carlo simulator to estimate the mean and standard deviation of the flux from the simulated star.

Food for thought - what do you notice if you run your simulator many times?



In [ ]:

    
sim_fluxes = # complete
for # complete

print("The mean flux = {:.3f} with variance = {:.3f}".format( # complete

Problem 4) PSF Flux measurement

In this problem we are going to use our knowledge of the PSF to estimate the flux of the star. We will compare these measurements to the aperture flux measurements above.

Problem 4a

Create the psf model, psf, which is equivalent to a noise-free star with fwhm = 5.



In [ ]:

    
psf = # complete

Problem 4b

Using the same parameters as problem 3, simulate a star and measure it's PSF flux.



In [ ]:

    
sim_star = simulate( # complete

psf_flux = # complete

print("The PSF flux is {:.3f}".format( # complete

Problem 4c

As before, write a Monte Carlo simulator to estimate the PSF flux of the star. How do your results compare to above?



In [ ]:

    
sim_fluxes = # complete
for # complete

print("The mean flux = {:.3f} with variance = {:.3f}".format( # complete

Challenge Problem

Problem C1

Simulate several stars with flux F and measure their flux to determine the "detection limit" relative to a sky brightness S.



In [ ]:

Problem C2

Simulate multiple stars and determine the minimum separation for which multiple stars can be detected as a function of FWHM. Is this only a function of FHWM?



In [ ]: