-
Notifications
You must be signed in to change notification settings - Fork 10
/
masking.jl
67 lines (55 loc) · 2.16 KB
/
masking.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
function make_mask(
rng::AbstractRNG, Nside, θpix;
edge_padding_deg = 2,
edge_rounding_deg = 1,
apodization_deg = 1,
ptsrc_radius_arcmin = 7, # roughly similar to Story et al. 2015
num_ptsrcs = round(Int, prod(Nside .* (1,1)) * (θpix/60)^2 * 120/100) # SPT-like density of sources
)
deg2npix(x) = round(Int, x/θpix*60)
arcmin2npix(x) = round(Int, x/θpix)
ptsrc = num_ptsrcs==0 ? 1 : .!bleed(sim_ptsrcs(rng, Nside, num_ptsrcs), arcmin2npix(ptsrc_radius_arcmin))
boundary = boundarymask(Nside, deg2npix(edge_padding_deg))
mask_array = if apodization_deg in (false, 0)
boundary .& ptsrc
else
apod_ptsrc = num_ptsrcs==0 ? 1 : cos_apod(ptsrc, arcmin2npix(ptsrc_radius_arcmin));
cos_apod(boundary, deg2npix(apodization_deg), deg2npix(edge_rounding_deg)) .* apod_ptsrc
end
FlatMap(Float32.(mask_array), θpix=θpix)
end
make_mask(args...; kwargs...) = make_mask(Random.default_rng(), args...; kwargs...)
# all padding/smoothing/etc... quantities below are in units of numbers of pixels
function boundarymask(Nside, pad)
m = fill(true, Nside .* (1,1,))
m[1:pad,:] .= false
m[:,1:pad] .= false
m[end-pad+1:end,:] .= false
m[:,end-pad+1:end] .= false
m
end
function bleed(img, w)
Ny,Nx = size(img)
nearest = getfield.(@ondemand(ImageMorphology.feature_transform)(img),:I)
[norm(nearest[i,j] .- [i,j]) < w for i=1:Ny,j=1:Nx]
end
function cos_apod(img, w, smooth_distance=false)
Ny,Nx = size(img)
nearest = getfield.(@ondemand(ImageMorphology.feature_transform)(.!img),:I)
distance = [norm(nearest[i,j] .- [i,j]) for i=1:Ny,j=1:Nx]
if smooth_distance!=false
distance = @ondemand(ImageFiltering.imfilter)(distance, @ondemand(ImageFiltering.Kernel.gaussian)(smooth_distance))
end
@. (1-cos(min(distance,w)/w*pi))/2
end
function round_edges(img, w)
.!(@ondemand(ImageFiltering.imfilter)(img, @ondemand(ImageFiltering.Kernel.gaussian)(w)) .< 0.5)
end
function sim_ptsrcs(rng,Nside,nsources)
Ny, Nx = Nside .* (1,1)
m = fill(false,Ny,Nx);
for i=1:nsources
m[rand(rng,1:Ny),rand(rng,1:Nx)] = true
end
m
end