Examples

Setup

You will need the following packages installed to replicate this tutorial

julia> ]add Distributions LACosmic Plots PSFModels

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))

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

imdata .= rand.(rng, Poisson.(imdata))

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

clean_image = copy(imdata)

# Add Nc fake cosmic rays
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
end

# Make a mask where the detected cosmic rays should be
end

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

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(
)
data.mask == mask