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

    # 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 = [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(
    imshow(data.mask, title="true cosmics", cbar=false),
    imshow(mask, title="detected cosmics", cbar=false),
    size=(700, 400)
)

data.mask == mask

true