# Examples

## Setup

You will need the following packages installed to replicate this tutorial

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

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)
(image = [188.58170768150052 214.09362006840362 … 204.15202016042792 174.58109503250085; 202.09914267128877 183.6972497338309 … 203.6666879254898 204.72346481491394; … ; 193.09513217045617 166.82428246193842 … 163.70511533689998 183.12325388575027; 228.6225661704681 189.5528435387213 … 202.46161451180063 199.53550046618867], clean_image = [188.58170768150052 214.09362006840362 … 204.15202016042792 174.58109503250085; 202.09914267128877 183.6972497338309 … 203.6666879254898 204.72346481491394; … ; 193.09513217045617 166.82428246193842 … 163.70511533689998 183.12325388575027; 228.6225661704681 189.5528435387213 … 202.46161451180063 199.53550046618867], 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(
)
data.mask == mask
true