Examples

Setup

You will need the following packages installed to replicate this tutorial

julia> ] # pressing ']' should drop you into pkg-mode
pkg> add Distributions LACosmic Plots PSFModels

Removing bad pixels with LACosmic.jl

First, let's create some fake data with Gaussian sources

using Distributions
using PSFModels: gaussian
using Random

function make_data(rng, N; N_sources=20, N_cosmics=20)
    imdata = fill(200.0, (N, N))

    # Add some fake sources
    for _ in 1:N_sources
        x = rand(rng, Uniform(1, N + 1))
        y = rand(rng, Uniform(1, N + 1))
        brightness = rand(rng, Uniform(1000, 30000)) / (2π * 3.5^2)
        model = gaussian(;x, y, fwhm=3.5, amp=brightness)
        imdata .+= model.(axes(imdata)...)
    end

    # Add the poisson noise
    imdata .= rand.(rng, Poisson.(imdata))

    # Add readnoise
    imdata .+= rand(rng, Normal(0, 10), (N, N))

    clean_image = copy(imdata)

    # Add Nc fake cosmic rays
    crmask = falses((N, N))
    for i in 1:N_cosmics
        cr_x = round(Int, rand(rng, Uniform(6, N - 5)))
        cr_y = round(Int, rand(rng, Uniform(6, N - 5)))
        cr_brightnesses = rand(rng, Uniform(1000, 30000))
        imdata[cr_y, cr_x] += cr_brightnesses
        crmask[cr_y, cr_x] = true
    end

    # Make a mask where the detected cosmic rays should be
    return (image=imdata, clean_image, mask=crmask)
end

rng = MersenneTwister(808)
data = make_data(rng, 201)

(image = [170.2087856720863 187.49469665521792 … 221.19896556405303 202.87099187261302; 211.55829508818854 177.30442832081616 … 196.66652116315373 217.1615384412527; … ; 207.7234068656278 239.5194657153742 … 178.88040560478606 203.36915666154326; 201.32532419058822 209.79986780607024 … 213.03558525423801 183.38472950015483], clean_image = [170.2087856720863 187.49469665521792 … 221.19896556405303 202.87099187261302; 211.55829508818854 177.30442832081616 … 196.66652116315373 217.1615384412527; … ; 207.7234068656278 239.5194657153742 … 178.88040560478606 203.36915666154326; 201.32532419058822 209.79986780607024 … 213.03558525423801 183.38472950015483], mask = Bool[0 0 … 0 0; 0 0 … 0 0; … ; 0 0 … 0 0; 0 0 … 0 0])

let's inspect it

using Plots

function imshow(image; kwargs...)
    axy, axx = axes(image)
    heatmap(axy, axx, image;
            aspect_ratio=1,
            ticks=false,
            xlim=extrema(axx),
            ylim=extrema(axy),
            kwargs...)
end

plot(
    imshow(log10.(data.clean_image), title="original image"),
    imshow(log10.(data.image), title="image w/cosmics"),
    size=(775, 350)
)

now we can clean it using lacosmic

using LACosmic

clean_image, mask = lacosmic(data.image, sigma_clip=6, contrast=5, neighbor_thresh=1)

plot(
    imshow(log10.(data.clean_image), title="original image"),
    imshow(log10.(clean_image), title="cleaned image"),
    size=(775, 350)
)

plot(
    imshow(data.mask, title="true cosmics", cbar=false),
    imshow(mask, title="detected cosmics", cbar=false),
    size=(700, 400)
)

data.mask == mask

true