Module cswt_tr_mod

module cswt_tr_mod

        ! Uses
    use s2_types_mod
    use s2_sky_mod
    use s2_vect_mod
    use cswt_error_mod
    use cswt_swav_mod

        ! Types
    public type cswt_tr

        ! Variables
    integer, public, parameter :: CSWT_TR_DIL_UNIT_RAD = 1
    integer, public, parameter :: CSWT_TR_DIL_UNIT_ARCMIN = 2
    character (len=*), private, parameter :: CSWT_TR_DIL_UNIT_STR_RAD = 'radian'
    character (len=*), private, parameter :: CSWT_TR_DIL_UNIT_STR_ARCMIN = 'arcmin'
    integer, public, parameter :: CSWT_TR_METHOD_FAST_FFT_REAL = 1
    integer, public, parameter :: CSWT_TR_METHOD_FAST_FFT = 2
    integer, public, parameter :: CSWT_TR_METHOD_FAST_DFT = 3
    integer, public, parameter :: CSWT_TR_METHOD_FAST_ISOTROPIC = 4
    integer, public, parameter :: CSWT_TR_METHOD_DIRECT = 5
    integer, public, parameter :: CSWT_TR_KERNEL_TYPE_SQR = 1
    integer, public, parameter :: CSWT_TR_KERNEL_TYPE_CIRC = 2
    integer, public, parameter :: CSWT_TR_MASK_GEN_CSWT = 1
    integer, public, parameter :: CSWT_TR_MASK_GEN_MORPH = 2
    integer, public, parameter :: CSWT_TR_THRES_NSIGMA_MODE_ABS = 1
    integer, public, parameter :: CSWT_TR_THRES_NSIGMA_MODE_ABOVE = 2
    integer, public, parameter :: CSWT_TR_THRES_NSIGMA_MODE_BELOW = 3

        ! Interfaces
    public interface cswt_tr_init
    public interface cswt_tr_const_multiply
    public interface cswt_tr_const_subtract
    public interface cswt_tr_wcoeff_sigma

        ! Subroutines and functions
    private function cswt_tr_init_data (swav_mother, dilation, lmax, mmax, n_gamma) result (tr)
    private function cswt_tr_init_file_dil (swav_mother, filename_dil, lmax, mmax, n_gamma) result (tr)
    private function cswt_tr_init_file (filename_swav, filename_dil, lmax, mmax, n_gamma, wav_name) result (tr)
    private function cswt_tr_init_read (filename) result (tr)
    private function cswt_tr_init_computed (dilation, wcoeff, lmax, mmax) result (tr)
    private function cswt_tr_init_copy (orig) result (copy)
    private function cswt_tr_init_merg (tr) result (tr_merg)
    public subroutine cswt_tr_free (tr)
    public subroutine cswt_tr_analysis (tr, data, method, message_status, time_it, filename_isotropic)
    private subroutine cswt_tr_anadirect (tr, data, message)
    private subroutine cswt_tr_anafast (tr, data, method, message, filename_isotropic)
    private subroutine cswt_tr_anafast_isotropic (tr, i_dil, sky, swav, message, filename)
    private subroutine cswt_tr_anafast_spcv (tr, i_dil, sky1, sky2, method_in, message)
    private subroutine cswt_tr_anafast_fft3d (x, backward)
    private function cswt_tr_anafast_imdft (tmmm, alpha, beta, gamma) result (t)
    public subroutine cswt_tr_norm (tr, norm, weight_in)
    public subroutine cswt_tr_multiply (tr1, tr2, weight_in)
    private subroutine cswt_tr_const_multiply_array (tr, mfactor)
    private subroutine cswt_tr_const_multiply_val (tr, mfactor)
    private subroutine cswt_tr_const_subtract_array (tr, val)
    private subroutine cswt_tr_const_subtract_val (tr, val)
    private subroutine cswt_tr_wcoeff_sigma_each (tr, n_effective_samples, mean, sigma)
    private subroutine cswt_tr_wcoeff_sigma_all (tr, n_effective_samples, mean, sigma)
    public subroutine cswt_tr_wcoeff_thres_nsigma (tr, nsigma, mode, n_effective_samples_in)
    public subroutine cswt_tr_localmax_ab (tr, i_dil, n_regions, max_val, max_loc, max_siz, filename_connected)
    private subroutine cswt_tr_connected_ab (tr, i_dil, i_gamma)
    private subroutine cswt_tr_equivclass (tab, n_equiv_pair, equiv)
    private subroutine cswt_tr_reallocate_mat (mat, x_axis)
    private function cswt_tr_present_vect_int (vect, val) result (pres)
    private function cswt_tr_present_mat_int (mat, val) result (pres)
    public subroutine cswt_tr_mask_nonzero (tr_mask, ncoeff_eff, ncoeff_tot)
    public subroutine cswt_tr_mask_nonzero_weight (tr_mask, prop)
    public subroutine cswt_tr_mask_invert (tr_mask)
    public subroutine cswt_tr_mask_apply (tr, tr_mask, display_in)
    public subroutine cswt_tr_mask_thres_nsigma (tr, nsigma, n_effective_samples, mode)
    public subroutine cswt_tr_mask_copy (tr_mask, sky_mask)
    public subroutine cswt_tr_mask_gen (tr_mask, sky_mask, mode, run_stage1_in, run_stage2_in, orig_mask_only_in, thres1_proportion_in, thres2_proportion_in, thres1_min_in, std_in)
    private subroutine cswt_tr_morph_erode_ab (tr, i_dil, i_gamma, kernel_size, kernel_type)
    private subroutine cswt_tr_morph_dilate_ab (tr, i_dil, i_gamma, kernel_size, kernel_type)
    public subroutine cswt_tr_smooth_gaussian (tr, kernel_size1, kernel_size2, std)
    private subroutine cswt_tr_smooth_gaussian_ab (tr, kernel_size1, kernel_size2, std, i_dil, i_gamma)
    public function cswt_tr_extract_wcoeff_sky (tr, i_dil, i_gamma, nside, interp, pix_scheme_in) result (sky)
    public subroutine cswt_tr_io_fits_write_wcoeff_sky (filename, tr, i_dil, i_gamma, nside, interp, pix_scheme_in, comment)
    public subroutine cswt_tr_io_fits_write_wcoeff (filename, tr, comment)
    private subroutine cswt_tr_io_fits_read_wcoeff (filename, tr)
    private subroutine cswt_tr_io_fits_error_check (status, halt)
    private subroutine cswt_tr_io_fits_exists (filename, status, exists)
    private subroutine cswt_tr_io_fits_del (filename, status)
    public subroutine cswt_tr_io_txt_dilation_write (filename, tr, unit)
    public subroutine cswt_tr_io_txt_dilation_read (filename, dilation)
    public function cswt_tr_get_init (tr) result (init)
    public function cswt_tr_get_swav_mother (tr) result (swav_mother)
    public subroutine cswt_tr_get_wcoeff (tr, wcoeff)
    public subroutine cswt_tr_get_dilation (tr, dilation)
    public function cswt_tr_get_lmax (tr) result (lmax)
    public function cswt_tr_get_mmax (tr) result (mmax)
    public function cswt_tr_get_n_alpha (tr) result (n_alpha)
    public function cswt_tr_get_n_beta (tr) result (n_beta)
    public function cswt_tr_get_n_gamma (tr) result (n_gamma)
    public function cswt_tr_get_n_dilation (tr) result (n_dilation)
    public function cswt_tr_get_wcoeff_status (tr) result (wcoeff_status)
    public function cswt_tr_get_swav_mother_status (tr) result (swav_mother_status)

end module cswt_tr_mod

Provides functionality to support and compute the directional spherical wavelet transform. In addition the mother wavelet, the dilation values and all wavelet parameters are stored herein. The resulting wavelet coefficients for each dilation are sampled on an ecp (equi-sampled) alpha-beta-gamma Euler angle grid. Functionality is also provided to convert the alpha-beta dimensions to Healpix sky representations and vice versa.

Author: J. D. McEwen (mcewen@mrao.cam.ac.uk)

Version: 0.1 - November 2004

Description of Types

cswt_tr

public type cswt_tr
    private
    logical :: init = .false.
    type (cswt_swav) :: swav_mother
    real (kind=s2_sp), allocatable, dimension (:,:,:,:) :: wcoeff
    real (kind=s2_sp), allocatable, dimension (:,:) :: dilation
    integer :: lmax = 0
    integer :: mmax = 0
    integer :: n_alpha = 0
    integer :: n_beta = 0
    integer :: n_gamma = 0
    integer :: n_dilation = 0
    logical :: wcoeff_status = .false.
    logical :: swav_mother_status = .false.
end type cswt_tr

Description of Variables

CSWT_TR_DIL_UNIT_RAD

integer, public, parameter :: CSWT_TR_DIL_UNIT_RAD = 1

Radian dilation unit specifier.

CSWT_TR_DIL_UNIT_ARCMIN

integer, public, parameter :: CSWT_TR_DIL_UNIT_ARCMIN = 2

Arcminute dilation unit specifier.

CSWT_TR_DIL_UNIT_STR_RAD

character (len=*), private, parameter :: CSWT_TR_DIL_UNIT_STR_RAD = 'radian'

CSWT_TR_DIL_UNIT_STR_ARCMIN

character (len=*), private, parameter :: CSWT_TR_DIL_UNIT_STR_ARCMIN = 'arcmin'

CSWT_TR_METHOD_FAST_FFT_REAL

integer, public, parameter :: CSWT_TR_METHOD_FAST_FFT_REAL = 1

Fast FFT analysis method specifier.

CSWT_TR_METHOD_FAST_FFT

integer, public, parameter :: CSWT_TR_METHOD_FAST_FFT = 2

Fast FFT analysis method specifier.

CSWT_TR_METHOD_FAST_DFT

integer, public, parameter :: CSWT_TR_METHOD_FAST_DFT = 3

Fast DFT analysis method specifier.

CSWT_TR_METHOD_FAST_ISOTROPIC

integer, public, parameter :: CSWT_TR_METHOD_FAST_ISOTROPIC = 4

Fast isotropic analysis method specifier.

CSWT_TR_METHOD_DIRECT

integer, public, parameter :: CSWT_TR_METHOD_DIRECT = 5

Direct analysis method specifier.

CSWT_TR_KERNEL_TYPE_SQR

integer, public, parameter :: CSWT_TR_KERNEL_TYPE_SQR = 1

Square kernel type specifier.

CSWT_TR_KERNEL_TYPE_CIRC

integer, public, parameter :: CSWT_TR_KERNEL_TYPE_CIRC = 2

Circular kernel type specifier.

CSWT_TR_MASK_GEN_CSWT

integer, public, parameter :: CSWT_TR_MASK_GEN_CSWT = 1

Generate mask from wavelet coefficients specifier.

CSWT_TR_MASK_GEN_MORPH

integer, public, parameter :: CSWT_TR_MASK_GEN_MORPH = 2

Generate mask by morphological operations specifier.

CSWT_TR_THRES_NSIGMA_MODE_ABS

integer, public, parameter :: CSWT_TR_THRES_NSIGMA_MODE_ABS = 1

Generate thresholded nsigma mask by considering absolte value.

CSWT_TR_THRES_NSIGMA_MODE_ABOVE

integer, public, parameter :: CSWT_TR_THRES_NSIGMA_MODE_ABOVE = 2

Generate thresholded nsigma mask by leaving values above threshold.

CSWT_TR_THRES_NSIGMA_MODE_BELOW

integer, public, parameter :: CSWT_TR_THRES_NSIGMA_MODE_BELOW = 3

Generate thresholded nsigma mask by leaving values below threshold.

Description of Interfaces

cswt_tr_init

public interface cswt_tr_init
    module procedure cswt_tr_init_data
    module procedure cswt_tr_init_file_dil
    module procedure cswt_tr_init_file
    module procedure cswt_tr_init_read
    module procedure cswt_tr_init_copy
    module procedure cswt_tr_init_merg
end interface cswt_tr_init

cswt_tr_const_multiply

public interface cswt_tr_const_multiply
    module procedure cswt_tr_const_multiply_array
    module procedure cswt_tr_const_multiply_val
end interface cswt_tr_const_multiply

cswt_tr_const_subtract

public interface cswt_tr_const_subtract
    module procedure cswt_tr_const_subtract_array
    module procedure cswt_tr_const_subtract_val
end interface cswt_tr_const_subtract

cswt_tr_wcoeff_sigma

public interface cswt_tr_wcoeff_sigma
    module procedure cswt_tr_wcoeff_sigma_each
    module procedure cswt_tr_wcoeff_sigma_all
end interface cswt_tr_wcoeff_sigma

Description of Subroutines and Functions

cswt_tr_init_data

private function cswt_tr_init_data (swav_mother, dilation, lmax, mmax, n_gamma) result (tr)
    type (cswt_swav), intent(in) :: swav_mother
    real (kind=s2_sp), intent(in), dimension (:,:) :: dilation
    integer, intent(in) :: lmax
    integer, intent(in) :: mmax
    integer, intent(in) :: n_gamma
    type (cswt_tr) :: tr
    ! Calls: cswt_error
end function cswt_tr_init_data

Initialise a tr structure from a mother wavelet, dilations and resolution parameters.

Notes:

Space for wcoeff is allocated when tr is initialised (before actually calculate coefficients).
Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.
n_gamma must be odd for the fast cswt.

Variables:

swav_mother: Mother wavelet.
dilation: Dilation values.
lmax: Alm lmax.
mmax: Alm mmax.
n_gamma: Number of gamma Euler orientation to consider in the directional transform.
tr: Returned tr structure initialised with a mother wavelet, dilations and resolution parameters.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_init_file_dil

private function cswt_tr_init_file_dil (swav_mother, filename_dil, lmax, mmax, n_gamma) result (tr)
    type (cswt_swav), intent(in) :: swav_mother
    character (len=*), intent(in) :: filename_dil
    integer, intent(in) :: lmax
    integer, intent(in) :: mmax
    integer, intent(in) :: n_gamma
    type (cswt_tr) :: tr
    ! Calls: cswt_error, cswt_tr_io_txt_dilation_read
end function cswt_tr_init_file_dil

Initialise a tr strucutre from a mother wavelet, dilations read from an input dilation text file and resolution parameters.

Notes:

Once dilations are read cswt_tr_init_data is called.
Space for wcoeff is allocated when tr is initialised (before actually calculate coefficients).
Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.
n_gamma must be odd for the fast cswt.

Variables:

swav_mother: Mother wavelet.
filename_dil: Name of input dilation file.
lmax: Alm lmax.
mmax: Alm mmax.
n_gamma: Number of gamma Euler orientation to consider in the directional transform.
tr: Returned tr structure initialised with a mother wavelet, dilations and resolution parameters.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_init_file

private function cswt_tr_init_file (filename_swav, filename_dil, lmax, mmax, n_gamma, wav_name) result (tr)
    character (len=*), intent(in) :: filename_swav
    character (len=*), intent(in) :: filename_dil
    integer, intent(in) :: lmax
    integer, intent(in) :: mmax
    integer, intent(in) :: n_gamma
    character (len=*), optional, intent(in) :: wav_name
    type (cswt_tr) :: tr
    ! Calls: cswt_error, cswt_swav_free
end function cswt_tr_init_file

Initialise a tr strucutre from a mother wavelet read from an input fits file, dilations read from an input dilation text file and resolution parameters.

Notes:

Once the wavelet is read cswt_tr_init_dil is called.
Space for wcoeff is allocated when tr is initialised (before actually calculate coefficients).
Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.
n_gamma must be odd for the fast cswt.

Variables:

filename_swav: Name of input fits sky file containing the mother wavelet.
filename_dil: Name of input dilation file.
lmax: Alm lmax.
mmax: Alm mmax.
n_gamma: Number of gamma Euler orientation to consider in the directional transform.
[wav_name]: Name of mother wavelet type.
tr: Returned tr structure initialised with a mother wavelet, dilations and resolution parameters.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_init_read

private function cswt_tr_init_read (filename) result (tr)
    character (len=*), intent(in) :: filename
    type (cswt_tr) :: tr
    ! Calls: cswt_error, cswt_tr_io_fits_read_wcoeff
end function cswt_tr_init_read

Wrapper to initalise a tr data structure from a file. The tr structure is read and initialised by the routine cswt_tr_io_fits_read_wcoeff.

Notes:

If a tr structure is initialised in this manner then no wavelet is stored and further analysis calls may not be made.
Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.

Variables:

filename: Name of wavelet coefficient fits file containing the tr data to be read.
tr: Returned tr structure initialised with the data contained in the input fits file.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_init_computed

private function cswt_tr_init_computed (dilation, wcoeff, lmax, mmax) result (tr)
    real (kind=s2_sp), intent(in), dimension (:,:) :: dilation
    real (kind=s2_sp), intent(in), dimension (:,:,:,:) :: wcoeff
    integer, intent(in) :: lmax
    integer, intent(in) :: mmax
    type (cswt_tr) :: tr
    ! Calls: cswt_error
end function cswt_tr_init_computed

Initialise a tr structure from precomputed wavelet coefficients (if, say, read from file).

Notes:

If a tr structure is initialised in this manner then no wavelet is stored and further analysis calls may not be made.
Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.

Variables:

dilation: Array of dilation values.
wcoeff: Array of wavelet coefficients.
lmax: Alm lmax.
mmax: Alm mmax.
tr: Returned tr strucutre initialised with the wavelet coefficients and dilations passed as input.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_init_copy

private function cswt_tr_init_copy (orig) result (copy)
    type (cswt_tr), intent(in) :: orig
    type (cswt_tr) :: copy
    ! Calls: cswt_error
end function cswt_tr_init_copy

Initialise a tr structure as a copy of another tr.

Variables:

orig: Original tr to be copied.
copy: Returned tr strucutre initialised with a copy of the original.

Notes:

Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_init_merg

private function cswt_tr_init_merg (tr) result (tr_merg)
    type (cswt_tr), intent(in), dimension (:) :: tr
    type (cswt_tr) :: tr_merg
    ! Calls: cswt_error
end function cswt_tr_init_merg

Initialise a tr structure by merging a number of tr structures that have only a single dilation.

Variables:

tr: Array of tr objects to be merged.
tr_merg: Returned merged tr strucutre initialised from the merged tr data.

Notes:

Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_free

public subroutine cswt_tr_free (tr)
    type (cswt_tr), intent(inout) :: tr
    ! Calls: cswt_error, cswt_swav_free
end subroutine cswt_tr_free

Free all data associated with an initialised tr structure and reset all other attributes.

Variables:

tr: The tr structure to be freed.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_analysis

public subroutine cswt_tr_analysis (tr, data, method, message_status, time_it, filename_isotropic)
    type (cswt_tr), intent(inout) :: tr
    type (s2_sky), intent(inout) :: data
    integer, optional, intent(in) :: method
    logical, optional, intent(in) :: message_status
    logical, optional, intent(in) :: time_it
    character (len=*), optional, intent(in) :: filename_isotropic
    ! Calls: cswt_error, cswt_tr_anadirect, cswt_tr_anafast, s2_sky_free, system_clock
end subroutine cswt_tr_analysis

Perform directional continuous spherical wavelet transform to compute wavelet coefficients. The actual method employed to compute wavelet coefficients depends on the method specifier.

Notes:

Data must have intent inout so can write alms to it.

Variables:

tr: Spherical wavelet structure to compute coefficients of. The computed coefficients are written to the tr%wcoeff array.
data: The data sky to take the wavelet transform of.
[method]: Optional method specifier to specify the method used to compute the spherical wavelet coefficients.
[message_status]: Logical to indicate whether messages are to be displyed during the analysis routine.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_anadirect

private subroutine cswt_tr_anadirect (tr, data, message)
    type (cswt_tr), intent(inout) :: tr
    type (s2_sky), intent(inout) :: data
    logical, optional, intent(in) :: message
    ! Calls: cswt_error, cswt_swav_dilate, cswt_swav_free, s2_sky_free, s2_sky_get_map, s2_sky_rotate
end subroutine cswt_tr_anadirect

Perform the direct directional spherical wavelet transform.

Variables:

tr: Spherical wavelet structure to compute coefficients of. The computed coefficients are written to the tr%wcoeff array.
data: The data sky to take the wavelet transform of.
[message_status]: Logical to indicate whether messages are to be displyed during the analysis routine.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_anafast

private subroutine cswt_tr_anafast (tr, data, method, message, filename_isotropic)
    type (cswt_tr), intent(inout) :: tr
    type (s2_sky), intent(inout) :: data
    integer, optional, intent(in) :: method
    logical, optional, intent(in) :: message
    character (len=*), optional, intent(in) :: filename_isotropic
    ! Calls: cswt_error, cswt_swav_compute_alm, cswt_swav_dilate, cswt_swav_free, cswt_tr_anafast_isotropic, cswt_tr_anafast_spcv, s2_sky_compute_alm, s2_sky_free
end subroutine cswt_tr_anafast

Perform the fast directional spherical wavelet transform.

Notes:

For case where alms are already defined (and map isn't) then set dilation to unitary; no dilation will actually be performed and analysis can continue.

Variables:

tr: Spherical wavelet structure to compute coefficients of. The computed coefficients are written to the tr%wcoeff array.
data: The data sky to take the wavelet transform of.
[method]: Optional method specifier to specify the method used to compute the spherical wavelet coefficients.
[message_status]: Logical to indicate whether messages are to be displyed during the analysis routine.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_anafast_isotropic

private subroutine cswt_tr_anafast_isotropic (tr, i_dil, sky, swav, message, filename)
    type (cswt_tr), intent(inout) :: tr
    integer, intent(in) :: i_dil
    type (s2_sky), intent(in) :: sky
    type (s2_sky), intent(in) :: swav
    logical, optional, intent(in) :: message
    character (len=*), optional, intent(in) :: filename
    ! Calls: cswt_error, s2_sky_compute_map, s2_sky_extract_ab, s2_sky_free, s2_sky_get_alm, s2_sky_write_map_file
end subroutine cswt_tr_anafast_isotropic

Perform the fast isotropic spherical convolution for a particular dilation (noting the the convolution is simply given by a product of spherical harmonic coefficients).

Notes:

Sizes are checked before here so no checking of size consistency is performed herein. (All skies must have same lmax and mmax as values passed to function.)
tr%wcoeff array must be already allocated and must be indexed from 0:N-1 for Euler angle dimensions.

Variables:

tr: Spherical wavelet structure to compute coefficients of. The computed coefficients are written to the tr%wcoeff array.
i_dil: Dilation index to store computed coefficient in.
sky: Data to convolve.
swav: Dilated wavelet to convolve.
[message_status]: Logical to indicate whether messages are to be displyed during the analysis routine.
[filename]: Name of fits file to save healpix pixelisation of wavelet coefficients in (if present, else not saved)

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_anafast_spcv

private subroutine cswt_tr_anafast_spcv (tr, i_dil, sky1, sky2, method_in, message)
    type (cswt_tr), intent(inout) :: tr
    integer, intent(in) :: i_dil
    type (s2_sky), intent(in) :: sky1
    type (s2_sky), intent(in) :: sky2
    integer, optional, intent(in) :: method_in
    logical, optional, intent(in) :: message
    ! Calls: cswt_error, cswt_tr_anafast_fft3d, rfftw3d_f77_create_plan, rfftwnd_f77_destroy_plan, rfftwnd_f77_one_complex_to_real, s2_dl_beta_operator, s2_sky_get_alm
end subroutine cswt_tr_anafast_spcv

Perform the fast spherical convolution for a particular dilation using the fast algorithm of Wandlet and Gorski.

Notes:

Sizes are checked before here so no checking of size consistency is performed herein. (All skies must have same lmax and mmax as values passed to function.)
tr%wcoeff array must be already allocated and must be indexed from 0:N-1 for Euler angle dimensions.

Variables:

tr: Spherical wavelet structure to compute coefficients of. The computed coefficients are written to the tr%wcoeff array.
i_dil: Dilation index to store computed coefficient in.
sky1: Data to convolve.
sky2: Dilated wavelet to convolve.
[method_in]: Optional method specifier to specify the method used to compute the spherical wavelet coefficients.
[message_status]: Logical to indicate whether messages are to be displyed during the analysis routine.

Author: D. J. Mortlock and J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_anafast_fft3d

private subroutine cswt_tr_anafast_fft3d (x, backward)
    complex (kind=s2_dpc), intent(inout), dimension (:,:,:) :: x
    logical, intent(in) :: backward
    ! Calls: complex_fft
end subroutine cswt_tr_anafast_fft3d

Perform 3D FFT using healpix fft routines (fftw if linked).

Notes:

Since the Fourier transform is separable the 3D transform is calculated by taking suitable 1D transforms over each dimension of the data.
No scaling performed - even on reverse transform.

Variables:

x: The 3D input data on input. The 3D FFT of the inputted data on exit.
backward: Logical specifying whether the forward FFT (.false.) or reverse inverse FFT (.true.) is performed.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_anafast_imdft

private function cswt_tr_anafast_imdft (tmmm, alpha, beta, gamma) result (t)
    complex (kind=s2_dpc), intent(in), dimension (:, :, :) :: tmmm
    real (kind=s2_dp), intent(in) :: alpha
    real (kind=s2_dp), intent(in) :: beta
    real (kind=s2_dp), intent(in) :: gamma
    real (kind=s2_sp) :: t
end function cswt_tr_anafast_imdft

Perform the modified inverse Fourier transform required to convert T_mmm to T_abg directly.

Notes:

t may be computed for any continuous values of alpha, beta and gamma.

Variables:

tmmm: Input data representing the modified 3d Fourier transform of required t_abg. (Now allocated array rather than pointer so indexed from 1 when accessed herein.)
alpha: Alpha Euler angle to calculate t for.
beta: Beta Euler angle to calculate t for.
gamma: Gamma Euler angle to calculate t for.
t: The modified inverse Fourier transform of tmmm for the specified angles.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_norm

public subroutine cswt_tr_norm (tr, norm, weight_in)
    type (cswt_tr), intent(in) :: tr
    real (kind=s2_sp), intent(out), dimension (:) :: norm
    logical, optional, intent(in) :: weight_in
    ! Calls: cswt_error, cswt_tr_free
end subroutine cswt_tr_norm

Compute the norm of the wavelet coefficients stored in the tr strucutre for each dilation. The norm is defined as the sum, over all Euler angles, of the squared coefficient values, i.e. norm(a) = 1/N_s * \sum_i ( k_i (W_i)^2 ) where W_i are the wavelet coefficients, k_i the weights (if applied) and N_s is the number of samples per dilation, i.e. N_s = N_alpha * N_beta * N_gamma. The weights may be optionally applied. Note that N_s=1 now

Notes:

Take note that the norm computed is scaled by one over the number of samples per dilation. -- Not any longer: must do scaling separately, since in some cases want to scale by number of effective coefficients, not total number.

Variables:

tr: The tr structure containing the wavelet coefficients to compute the dilation of.
norm: Computed norm values for each dilation.
[weight_in]: Logical to specified whether pixel weights are to be applied. If not specified then WEIGHT_DEFAULT is assumed.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_multiply

public subroutine cswt_tr_multiply (tr1, tr2, weight_in)
    type (cswt_tr), intent(inout) :: tr1
    type (cswt_tr), intent(in) :: tr2
    logical, optional, intent(in) :: weight_in
    ! Calls: cswt_error
end subroutine cswt_tr_multiply

Multiply two wavelet coefficients contained in tr structure arrays. Weight the multiplication for each pixel by the size of the pixel on the sky if desired. The product coefficients are written to the tr1 structure wavelet coefficient array on exit.

Variables:

tr1: First tr structure containing the coefficients to multiply. Contains product of coefficients on exit.
tr2: Second tr structure containing the coefficients to multiply. Unchanged on exit.
[weight_in]: Logical to specified whether pixel weights are to be applied. If not specified then WEIGHT_DEFAULT is assumed.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_const_multiply_array

private subroutine cswt_tr_const_multiply_array (tr, mfactor)
    type (cswt_tr), intent(inout) :: tr
    real (kind=s2_sp), intent(in), dimension (:) :: mfactor
    ! Calls: cswt_error
end subroutine cswt_tr_const_multiply_array

Multiply the wavelet coefficients of the tr structure for a given dilation by the value of mfactor for the corresponding dilation.

Variables:

tr: tr structure containing the coefficients to multiply.
mfactor: Array containing multiplication factor for each dilation.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_const_multiply_val

private subroutine cswt_tr_const_multiply_val (tr, mfactor)
    type (cswt_tr), intent(inout) :: tr
    real (kind=s2_sp), intent(in) :: mfactor
    ! Calls: cswt_error
end subroutine cswt_tr_const_multiply_val

Multiply the wavelet coefficients of the tr structure by a constant value (for all dilations).

Variables:

tr: tr structure containing the coefficients to multiply.
mfactor: Constant value to multiply all coefficients by (for all dilations).

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_const_subtract_array

private subroutine cswt_tr_const_subtract_array (tr, val)
    type (cswt_tr), intent(inout) :: tr
    real (kind=s2_sp), intent(in), dimension (:) :: val
    ! Calls: cswt_error
end subroutine cswt_tr_const_subtract_array

Subtract the wavelet coefficients of the tr structure for a given dilation by the value of val for the corresponding dilation.

Variables:

tr: tr structure containing the coefficients to subtract.
val: Array containing subtraction values for each dilation.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_const_subtract_val

private subroutine cswt_tr_const_subtract_val (tr, val)
    type (cswt_tr), intent(inout) :: tr
    real (kind=s2_sp), intent(in) :: val
    ! Calls: cswt_error
end subroutine cswt_tr_const_subtract_val

Subtract the wavelet coefficients of the tr structure by a constant value (for all dilations).

Variables:

tr: tr structure containing the coefficients to subtract.
val: Constant value to subtract from all coefficients (for all dilations).

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_wcoeff_sigma_each

private subroutine cswt_tr_wcoeff_sigma_each (tr, n_effective_samples, mean, sigma)
    type (cswt_tr), intent(in) :: tr
    integer, intent(in), dimension (:) :: n_effective_samples
    real (kind=s2_sp), intent(out), dimension (:) :: mean
    real (kind=s2_sp), intent(out), dimension (:) :: sigma
    ! Calls: cswt_error
end subroutine cswt_tr_wcoeff_sigma_each

Compute the mean and standard deviation of the wavelet coefficients of the tr strucutre for each dilation.

Notes:

If tr wavelet coefficients masked, sigma and mean should be computed only over the remaining pixels. If these pixels are zero, then they won't compute to the calculation of mean and sigma, however the number of samples to divide by should be the effective number of pixels. Thus divide by this value, for each dilation, that is passed as an input.

Variables:

tr: tr structure containing the coefficients to calculate the mean and standard deviation of.
n_effective_samples: Number of effective samples per dilation for previously applied mask, if not masked just give full number of Euler angle samples (for each dilation - which of course will be an array of the same number (n_alpha*n_beta*n_gamma)).
mean: Computed mean array for each dilation.
sigma: Computed standard deviation array for each dilation.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_wcoeff_sigma_all

private subroutine cswt_tr_wcoeff_sigma_all (tr, n_effective_samples, mean, sigma)
    type (cswt_tr), intent(in) :: tr
    integer, intent(in) :: n_effective_samples
    real (kind=s2_sp), intent(out) :: mean
    real (kind=s2_sp), intent(out) :: sigma
    ! Calls: cswt_error
end subroutine cswt_tr_wcoeff_sigma_all

Compute the mean and standard deviation of the wavelet coefficients of the tr strucutre over all dilations.

Notes:

If tr wavelet coefficients masked, sigma and mean should be computed only over the remaining pixels. If these pixels are zero, then they won't compute to the calculation of mean and sigma, however the number of samples to divide by should be the effective number of pixels. Thus divide by this value that is passed as an input.

Variables:

tr: tr structure containing the coefficients to calculate the mean and standard deviation of.
n_effective_samples: Number of effective samples, if not masked just give full number of samples: size(tr%wcoeff).
mean: Computed mean over all dilations.
sigma: Computed standard deviation over all dilations.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_wcoeff_thres_nsigma

public subroutine cswt_tr_wcoeff_thres_nsigma (tr, nsigma, mode, n_effective_samples_in)
    type (cswt_tr), intent(inout) :: tr
    real (kind=s2_sp) :: nsigma
    integer, intent(in) :: mode
    integer, optional, intent(in) :: n_effective_samples_in
    ! Calls: cswt_error, cswt_tr_wcoeff_sigma_all
end subroutine cswt_tr_wcoeff_thres_nsigma

Threshold wavelet coefficients. Coefficients that remain keep their original value and are not set to 1.0. Differs to cswt_tr_mask_thres_nsigma routine since the same nsigma is used over all dilations.

Notes:

If tr wavelet coefficients masked, sigma and mean should be computed only over the remaining pixels. If these pixels are zero, then they won't compute to the calculation of mean and sigma, however the number of samples to divide by should be the effective number of pixels. Thus divide by this value that is passed as an input.

Variables:

tr: tr structure containing the coefficients to be thresholded. On exit, tr structure containing thresholded coefficients.
n_effective_samples: Number of effective samples, if not masked just give full number of samples: size(tr%wcoeff).
mode: Mode of thresholding to perform. If mode=CSWT_TR_THRES_NSIGMA_MODE_ABS then the absolute value of coefficients are considered and those coefficients below (in absolute value) nsigma * sigma are set to zero. If mode=CSWT_TR_THRES_NSIGMA_MODE_ABOVE then those coefficients below nsigma * sigma are set to zero. If mode=CSWT_TR_THRES_NSIGMA_MODE_BELOW then those coefficients above -nsigma * sigma (note minus sign) are set to zero.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_localmax_ab

public subroutine cswt_tr_localmax_ab (tr, i_dil, n_regions, max_val, max_loc, max_siz, filename_connected)
    type (cswt_tr), intent(in) :: tr
    integer, intent(in) :: i_dil
    integer, allocatable, intent(out), dimension (:) :: n_regions
    real (kind=s2_sp), allocatable, intent(out), dimension (:,:) :: max_val
    integer, allocatable, intent(out), dimension (:,:,:) :: max_loc
    integer, allocatable, intent(out), dimension (:,:) :: max_siz
    character (len=*), optional, intent(in) :: filename_connected
    ! Calls: cswt_error, cswt_tr_connected_ab, cswt_tr_free, cswt_tr_io_fits_write_wcoeff
end subroutine cswt_tr_localmax_ab

Find local maxima in a single ab slice of the wavelet coefficients (performed for each orientation separately). The input coefficients should already be thresholded so that only local non-zero regions remain for which the maximum of each is found.

Notes:

Note the tr structure passed should be thresholded so only need to find local maximum in each connected component region.
Data arrays of maximum local values and locations are allocated herein, since no prior knowledge of number of objects contained in map. These arrays must be freed by the calling routine.
Both the dilation (i_dil) and orientation (i_gamma) are indexed from 1:N.
The algorithm works by finding the connected components in thresholded wavelet coefficient map. The maximum of each connected component is taken as the local maximum.

Variables:

tr: Thresholded tr structure containing non-zero local regions to find maxima of.
i_dil: Dilation index of tr%wcoeff map to consider.
i_gamma: Gamma Euler angle index of tr%wcoeff to consider (indexed from 1:N).
max_val: List of local maximum values found.
max_loc: List of locations of local maximum values found.
[filename_connected]: If present save cswt file of connected components.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_connected_ab

private subroutine cswt_tr_connected_ab (tr, i_dil, i_gamma)
    type (cswt_tr), intent(inout) :: tr
    integer, intent(in) :: i_dil
    integer, intent(in) :: i_gamma
    ! Calls: cswt_error, cswt_tr_equivclass, cswt_tr_reallocate_mat
end subroutine cswt_tr_connected_ab

Find wavelet coefficient connected components (defined by non-zero pixels) and overwrite with integer label for each connected component.

Notes:

Both the dilation (i_dil) and orientation (i_gamma) are indexed from 1:N.
Connected components defined by non-zero pixels (non-zero values are irrelevant).

Variables:

tr: Thresholded tr structure to find connected components of. On exit, coefficients contain integer values (although real precision) labelling connected components.
i_dil: Dilation index of tr%wcoeff map to consider.
i_gamma: Gamma Euler angle index of tr%wcoeff to consider (indexed from 1:N).

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_equivclass

private subroutine cswt_tr_equivclass (tab, n_equiv_pair, equiv)
    integer, intent(in), dimension (:,:) :: tab
    integer, intent(in) :: n_equiv_pair
    integer, allocatable, intent(out), dimension (:,:) :: equiv
    ! Calls: cswt_error, cswt_tr_reallocate_mat
end subroutine cswt_tr_equivclass

Resolve a table of equivalence pairs into a table of equivalence classes.

Notes:

Note memory for equic is allocated here and must be freed be calling routine.

Variables:

tab: Original table of equivalence pairs. Table may be larger than actual number of pairs - unused cells filled with zeros.
n_equiv_pair: Number of equivalence pairs in tab (number of non-zero rows of table).
equiv: Table of equivalence classes. Table may be larger than required - unused cells filled with zeros.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_reallocate_mat

private subroutine cswt_tr_reallocate_mat (mat, x_axis)
    integer, allocatable, dimension (:,:) :: mat
    logical, intent(in) :: x_axis
    ! Calls: cswt_error
end subroutine cswt_tr_reallocate_mat

Dynamically reallocate a 2D array to increase its size, either in the x or y direction.

Notes:

Original mat must be allocatable.
Current implementation is only for integer mat.

Variables:

mat: Original matrix to increase the size of. On exit, mat contains the same data but has a larger size.
x_axis: Logical indication whether to increase the size of mat in the x (down) direction. If false then size is increased in the y (across) direction.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_present_vect_int

private function cswt_tr_present_vect_int (vect, val) result (pres)
    integer, intent(in), dimension (:) :: vect
    integer, intent(in) :: val
    logical :: pres
end function cswt_tr_present_vect_int

Check if a value is present in a vector.

Notes:

Could simply subtract value from vect, then take product of all elements. Result will be zero if value is present. Errors with this approach however (error with product intrinsic function?), so search directly.

Variables:

vect: Vector to search for value in.
val: Test value to look for in vector.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_present_mat_int

private function cswt_tr_present_mat_int (mat, val) result (pres)
    integer, intent(in), dimension (:,:) :: mat
    integer, intent(in) :: val
    logical :: pres
end function cswt_tr_present_mat_int

Check if a value is present in a 2D matrix.

Notes:

Could simply subtract value from mat, then take product of all elements. Result will be zero if value is present. Errors with this approach however (error with product intrinsic function?), so search directly.

Variables:

mat: matrix to search for value in.
val: Test value to look for in vector.

Author: J. D. McEwen

Version: 0.1 - April 2005

cswt_tr_mask_nonzero

public subroutine cswt_tr_mask_nonzero (tr_mask, ncoeff_eff, ncoeff_tot)
    type (cswt_tr), intent(in) :: tr_mask
    integer, intent(out), dimension (:,:) :: ncoeff_eff
    integer, intent(out) :: ncoeff_tot
    ! Calls: cswt_error
end subroutine cswt_tr_mask_nonzero

Compute the number of non-zero values (i.e. effective remaining coefficients) in a given mask. (Computed for each dilation and orientation, i.e. for each alpha-beta `sky' of the mask.)

Variables:

tr_mask: tr structure containing the mask in the tr_mask%wcoeff array.
ncoeff_eff: Number of effective coefficients in the mask for each dilation and orientation, i.e. for each alpha-beta `sky'.
ncoeff_tot: Total number of coefficients in a single alpha-beta mask.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_mask_nonzero_weight

public subroutine cswt_tr_mask_nonzero_weight (tr_mask, prop)
    type (cswt_tr), intent(in) :: tr_mask
    real (kind=s2_sp), intent(out), dimension (:,:) :: prop
    ! Calls: cswt_error
end subroutine cswt_tr_mask_nonzero_weight

Compute the proportion of weighted coefficients in a given mask. (Computed for each dilation and orientation, i.e. for each alpha-beta `sky' of the mask.)

Notes:

The proportion is computed from . \sum_{all i s.t. mask_i > thres} w_i . ------------------------------------ . \sum_{all j} w_j for the map corresponding to each dilation and orientation.

Variables:

tr_mask: tr structure containing the mask in the tr_mask%wcoeff array.
prop: Proportion of effective weighted coefficients for each dilation and orientation, i.e. for each alpha-beta `sky'.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_mask_invert

public subroutine cswt_tr_mask_invert (tr_mask)
    type (cswt_tr), intent(inout) :: tr_mask
    ! Calls: cswt_error
end subroutine cswt_tr_mask_invert

Invert a coefficient mask by converting ones to zeros, and zeros to ones.

Variables:

tr_mask: Mask to invert. On output contains ones in the place of original zeros, and zeros in the place of original ones.

Author: J. D. McEwen

Version: 0.1 - January 2005

cswt_tr_mask_apply

public subroutine cswt_tr_mask_apply (tr, tr_mask, display_in)
    type (cswt_tr), intent(inout) :: tr
    type (cswt_tr), intent(in) :: tr_mask
    logical, optional, intent(in) :: display_in
    ! Calls: cswt_error
end subroutine cswt_tr_mask_apply

Apply an extended coefficient mask by multiplying the wavelet coefficients of the specified tr structure with the mask (contained in the wavelet coefficients of the tr_mask structure).

Notes:

Important: Only specify the display flag if the output may is only used for display purposes. If subsequent processing is applied erroneous results amy occur.
Obviously the sizes of the tr and tr_mask structure must be the same.

Variables:

tr: Wavelet transform structure to be masked.
tr_mask: Mask to apply. Mask is stored in the tr_mask%wcoeff array.
[display]: Logical to specified whether a display mask is to be applied. If that is the case then the pixels in the tr structure corresponding to the zeros with be replaced with a large negative value so that they appear grey in output images.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_mask_thres_nsigma

public subroutine cswt_tr_mask_thres_nsigma (tr, nsigma, n_effective_samples, mode)
    type (cswt_tr), intent(inout) :: tr
    real (kind=s2_sp), intent(in) :: nsigma
    integer, intent(in), dimension (:) :: n_effective_samples
    integer, intent(in) :: mode
    ! Calls: cswt_error, cswt_tr_wcoeff_sigma
end subroutine cswt_tr_mask_thres_nsigma

Produce a mask from the wavelet coefficients contained in the tr structure by thresholding the coefficients to binary 1,0 values. A range of thresholding 'directions' may be performed (see mode below), nevertheless each is realitive to the nsigma * sigma(a) (where sigma(a) is computed herein).

Notes:

Note that on output the binary {0,1} mask is stored in the tr wavelet coefficient array.
If tr wavelet coefficients masked, sigma and mean should be computed only over the remaining pixels. If these pixels are zero, then they won't compute to the calculation of mean and sigma, however the number of samples to divide by should be the effective number of pixels. Thus divide by this value, for each dilation, that is passed as an input.

Variables:

tr: Strucutre containing wavelet coefficients to be thresholded to construct the mask. On exit, the binary mask is contained in the tr structure wavelet coefficient array.
nsigma: Number of sigma threshold is set at.
n_effective_samples: Number of effective samples per dilation for previously applied mask, if not masked just give full number of Euler angle samples (for each dilation - which of course will be an array of the same number (n_alpha*n_beta*n_gamma)).
mode: Mode of thresholding to perform. If mode=CSWT_TR_THRES_NSIGMA_MODE_ABS then the absolute value of coefficients are considered and those coefficients above (in absolute value) nsigma * sigma(a) are set to one, and the remainder to zero. If mode=CSWT_TR_THRES_NSIGMA_MODE_ABOVE then those coefficients above nsigma * sigma(a) are set to one, and the remainder to zero. If mode=CSWT_TR_THRES_NSIGMA_MODE_BELOW then those coefficients below -nsigma * sigma(a) (note minus sign) are set to one, and the remainder to zero.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_mask_copy

public subroutine cswt_tr_mask_copy (tr_mask, sky_mask)
    type (cswt_tr), intent(inout) :: tr_mask
    type (s2_sky), intent(in), dimension (:) :: sky_mask
    ! Calls: cswt_error, s2_sky_extract_ab
end subroutine cswt_tr_mask_copy

Copy masks defined on the sky (Healpix format) for each dilation to a cswt tr structure mask defined in the spherical wavelet coefficient domain.

Notes:

The wavelet coefficients are not actually required in the algorithm, however the wavelet transform structure of the correct size must be initialised. To overcome this problem the tr_mask structure containing the wavelet coefficient masks is often used.

Variables:

tr_mask: Spherical wavelet transform of sky mask. Overwritten with extended coefficient mask on output. (Note that this strucutre doesn't actually have to have the wavelet coefficients of the original sky map mask, however it is easier to just supply that.)
sky_mask: Array of sky masks for each dilation that are to be convert to the wavelet coefficient domain.

Author: J. D. McEwen

Version: 0.1 - December 2004

cswt_tr_mask_gen

public subroutine cswt_tr_mask_gen (tr_mask, sky_mask, mode, run_stage1_in, run_stage2_in, orig_mask_only_in, thres1_proportion_in, thres2_proportion_in, thres1_min_in, std_in)
    type (cswt_tr), intent(inout) :: tr_mask
    type (s2_sky), intent(in) :: sky_mask
    integer, intent(in) :: mode
    logical, optional, intent(in) :: run_stage1_in
    logical, optional, intent(in) :: run_stage2_in
    logical, optional, intent(in) :: orig_mask_only_in
    real (kind=s2_sp), optional, intent(in) :: thres1_proportion_in
    real (kind=s2_sp), optional, intent(in) :: thres2_proportion_in
    real (kind=s2_sp), optional, intent(in) :: thres1_min_in
    real (kind=s2_sp), optional, intent(in), dimension (:) :: std_in
    ! Calls: cswt_error, cswt_tr_morph_dilate_ab, cswt_tr_morph_erode_ab, cswt_tr_smooth_gaussian_ab, s2_sky_extract_ab
end subroutine cswt_tr_mask_gen

Generate an extended coefficient exclusion mask from an original sky mask and the spherical wavelet transform of the mask.

Notes:

This routine may perform two different methods for generating the extended mask: (a) from spherical wavelet coefficients or (b) by using morphological operations. These are subsequently described.
(a) The extended coefficient exclusion mask is generated by taking the transform of the original sky mask. The only non-zero wavelet coefficients are those that are distorted by the mask boundary. These distorted coefficients are detected by thresholding. An intermediate stage of smoothing (by convolving with a Gaussian kernel) is performed to smooth out whispy regions that otherwise remain when the wavelet support is exactly half contained in the masked region and half not. A second stage of thresholding is performed to remove these regions once blurred.
(b) The extended coefficient exclusion mask is generated by performing morphological operations on the original mask. The wavelet coefficients are not requried, however the wavelet transform structure of the correct size must be initialised. To overcome this problem the tr_mask structure is actually used (since it is already available) and simply overwritten with the new extended mask. The method itself involves first removing the point source mask regions by opening (morphological dialtion followed by erosion) the original mask. The glactic plane exclusion region then remains. This is extended simply by a morphological erosion operation where the erosion kernel size is given by the 2*effective size of the wavelet on the sky. A circular kernel shape is used.
Additional important note: Morphological mask generation is only valid for the Mexhat, Morlet and Butterfly wavelets (or other wavelets that are defined to have the same effective size on the sky). Modifications must be made for wavelets with a different effective size on the sky.

Variables:

tr_mask: Spherical wavelet transform of sky mask. Overwritten with extended coefficient mask on output.
sky_mask: Original sky mask to be extended.
mode: Integer to specifier what method of generating extended coefficient mask to use.
[run_stage1]: Logical to specify whether to perform first stage of thresholding.
[run_stage2]: Logical to specify whether to perform second stage of thresholding.
[orig_mask_only_in]: Logical to specify if to just convert the original sky mask to the wavelet domain without extending it.
[thres1_proportion_in]: Threshold proportion for the first stage of dynamic thresholding.
[thres2_proportion_in]: Threshold proportion for the second stage of dynamic thresholding.
[thres1_min_in]: Minimum threshold for the first stage of dynamic thresholding.
[std_in]: Standard deviation of Gaussian smoothing for each dilation index. If not specified then calculated by default method from dilation values.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_morph_erode_ab

private subroutine cswt_tr_morph_erode_ab (tr, i_dil, i_gamma, kernel_size, kernel_type)
    type (cswt_tr), intent(inout) :: tr
    integer, intent(in) :: i_dil
    integer, intent(in) :: i_gamma
    integer, intent(in) :: kernel_size
    integer, intent(in) :: kernel_type
    ! Calls: cswt_error
end subroutine cswt_tr_morph_erode_ab

Morphologically erode an ecp alpha-beta map with either a square or circular kernel.

Notes:

The input tr wavelet coefficient alpha-beta map must be a binary map (the ones in the map are eroded, i.e. zero component dilated).
Kernel size variable is positive size (so total size of kernel is 2*kernel_size + 1).
Both the dilation (i_dil) and orientation (i_gamma) are indexed from 1:N.
Overwrites the tr alpha and beta components with the morphologically dilated values.

Variables:

tr: tr structure containing wavelet coefficient that are morphologically eroded.
i_dil: Dilation index of tr%wcoeff map to erode.
i_gamma: Gamma Euler angle index of tr%wcoeff to erode.
kernel_size: Size of 2D kernel in alpha and beta Euler angle dimension (must be odd).
kernel_type: Specifier to indicate shape of kernel to apply.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_morph_dilate_ab

private subroutine cswt_tr_morph_dilate_ab (tr, i_dil, i_gamma, kernel_size, kernel_type)
    type (cswt_tr), intent(inout) :: tr
    integer, intent(in) :: i_dil
    integer, intent(in) :: i_gamma
    integer, intent(in) :: kernel_size
    integer, intent(in) :: kernel_type
    ! Calls: cswt_error
end subroutine cswt_tr_morph_dilate_ab

Morphologically dilate an ecp alpha-beta map with either a square or circular kernel.

Notes:

The input tr wavelet coefficient alpha-beta map must be a binary map (the ones in the map are dilated).
Kernel size variable is positive size (so total size of kernel is 2*kernel_size + 1).
Both the dilation (i_dil) and orientation (i_gamma) are indexed from 1:N.
Overwrites the tr alpha and beta components with the morphologically dilated values.

Variables:

tr: tr structure containing wavelet coefficient that are morphologically dilated.
i_dil: Dilation index of tr%wcoeff map to dilate.
i_gamma: Gamma Euler angle index of tr%wcoeff to dilate.
kernel_size: Size of 2D kernel in alpha and beta Euler angle dimension (must be odd).
kernel_type: Specifier to indicate shape of kernel to apply.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_smooth_gaussian

public subroutine cswt_tr_smooth_gaussian (tr, kernel_size1, kernel_size2, std)
    type (cswt_tr), intent(inout) :: tr
    integer, intent(in) :: kernel_size1
    integer, intent(in) :: kernel_size2
    real (kind=s2_sp), intent(in), dimension (:,:) :: std
    ! Calls: cswt_error, cswt_tr_smooth_gaussian_ab
end subroutine cswt_tr_smooth_gaussian

Smooth all tr%wcoeff maps over alpha-beta for each dilation and gamma Euler orientation.

Notes:

Overwrites the tr wavelet coefficients with the smoothed values.

Variables:

tr: tr structure containing wavelet coefficient that are smoothed.
kernel_size1: Size of 2D kernel in alpha Euler angle dimension (must be odd).
kernel_size2: Size of 2D kernel in beta Euler angle dimension (must be odd).
std: Standard deviation of Gaussian used to smooth. Array indexed over i_dil and i_gamma containing the standard deviation for each alpha-beta map.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_smooth_gaussian_ab

private subroutine cswt_tr_smooth_gaussian_ab (tr, kernel_size1, kernel_size2, std, i_dil, i_gamma)
    type (cswt_tr), intent(inout) :: tr
    integer, intent(in) :: kernel_size1
    integer, intent(in) :: kernel_size2
    real (kind=s2_sp), intent(in) :: std
    integer, intent(in) :: i_dil
    integer, intent(in) :: i_gamma
    ! Calls: cswt_error
end subroutine cswt_tr_smooth_gaussian_ab

Smooth an ecp alpha-beta map by convolving it with a Gaussian kernel.

Notes:

Both the dilation (i_dil) and orientation (i_gamma) are indexed from 1:N.
Overwrites the tr alpha and beta components with the smoothed values.

Variables:

tr: tr structure containing wavelet coefficient that are smoothed.
kernel_size1: Size of 2D kernel in alpha Euler angle dimension (must be odd).
kernel_size2: Size of 2D kernel in beta Euler angle dimension (must be odd).
std: Standard deviation of Gaussian used to smooth.
i_dil: Dilation index of tr%wcoeff map to smooth.
i_gamma: Gamma Euler angle index of tr%wcoeff to smooth.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_extract_wcoeff_sky

public function cswt_tr_extract_wcoeff_sky (tr, i_dil, i_gamma, nside, interp, pix_scheme_in) result (sky)
    type (cswt_tr), intent(in) :: tr
    integer, intent(in) :: i_dil
    integer, intent(in) :: i_gamma
    integer, intent(in) :: nside
    logical, intent(in) :: interp
    integer, optional, intent(in) :: pix_scheme_in
    type (s2_sky) :: sky
    ! Calls: cswt_error
end function cswt_tr_extract_wcoeff_sky

Extract a sky representation of the ecp (equi-sampled) spherical wavelet coefficients for a given dilation and orientation.

Notes:

Initialises a new sky form the wavelet coefficients of the tr structure.
The returned sky is subsequently independed of the tr structure and should be freed by the calling routine at some point.
Both the dilation (i_dil) and orientation (i_gamma) are indexed from 1:N.

Variables:

tr: The tr structure containing the computed spherical wavelet coefficients that are to be written to a sky representation.
i_dil: The dilation index of the spherical wavelet coefficients to write the sky representation of.
i_gamma: The orientation index of the spherical wavelet coefficients to write the sky representation of.
nside: The Healpix nside resolution parameter of the constructed sky.
interp: Logical to specify whether linear interpolation is to be performed.
[pix_scheme_in]: Pixelisation scheme of initialised sky. If not specified then the default pixelisation scheme is used.
sky: Sky representation of the wavelet coefficients for the specified dilation and orientation.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_fits_write_wcoeff_sky

public subroutine cswt_tr_io_fits_write_wcoeff_sky (filename, tr, i_dil, i_gamma, nside, interp, pix_scheme_in, comment)
    character (len=*), intent(in) :: filename
    type (cswt_tr), intent(in) :: tr
    integer, intent(in) :: i_dil
    integer, intent(in) :: i_gamma
    integer, intent(in) :: nside
    logical, intent(in) :: interp
    integer, optional, intent(in) :: pix_scheme_in
    character (len=*), optional, intent(in) :: comment
    ! Calls: cswt_error, s2_sky_free, s2_sky_write_map_file
end subroutine cswt_tr_io_fits_write_wcoeff_sky

Extract a sky representation of the ecp (equi-sampled) spherical wavelet coefficients for a given dilation and orientation and write them directly to a fits sky file.

Notes:

Both the dilation (i_dil) and orientation (i_gamma) are indexed from 1:N.

Variables:

filename: Name of the output fits file to write the sky representation of the wavelet coefficients to.
tr: The tr structure containing the computed spherical wavelet coefficients that are to be written fits file sky representation.
i_dil: The dilation index of the spherical wavelet coefficients to write the sky representation of.
i_gamma: The orientation index of the spherical wavelet coefficients to write the sky representation of.
nside: The Healpix nside resolution parameter of the constructed sky.
interp: Logical to specify whether linear interpolation is to be performed.
[pix_scheme_in]: Pixelisation scheme of initialised sky. If not specified then the default pixelisation scheme is used.
[comment]: Optional comment to add to the header of the output fits file.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_fits_write_wcoeff

public subroutine cswt_tr_io_fits_write_wcoeff (filename, tr, comment)
    character (len=*), intent(in) :: filename
    type (cswt_tr), intent(in) :: tr
    character (len=*), optional, intent(in) :: comment
    ! Calls: cswt_error, cswt_tr_io_fits_error_check, cswt_tr_io_fits_exists, ftclos, ftfiou, ftgiou, ftibin, ftiimg, ftinit, ftp3de, ftpcle, ftpcom, ftpdat, ftphpr, ftpkye, ftpkyj, ftpkys, s2_sky_admiss_dil, s2_sky_free, s2_sky_get_param
end subroutine cswt_tr_io_fits_write_wcoeff

Write a tr structure to an output fits file. Both the wavelet coefficients and the dilation values are written (and other additional wavelet transform parameters).

Variables:

filename: Name of the output fits file to write the tr structure data to.
tr: The tr structure containing the data to be written to the output fits file.
[comment]: Optional comment string to be added to the output fits file header if present.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_fits_read_wcoeff

private subroutine cswt_tr_io_fits_read_wcoeff (filename, tr)
    character (len=*), intent(in) :: filename
    type (cswt_tr), intent(out) :: tr
    ! Calls: cswt_error, cswt_tr_io_fits_error_check, cswt_tr_io_fits_exists, ftclos, ftfiou, ftg3de, ftgcve, ftgiou, ftgkye, ftgkyj, ftmahd, ftopen, ftthdu
end subroutine cswt_tr_io_fits_read_wcoeff

Read a fits wavelet coefficient file and allocate a new tr structure with the data read.

Notes:

The tr data structure is allocated herein by calling cswt_tr_init_computed
No wavelet specific attributes are read from headers here since the wavelet itself is not contained in the file. In some cases the file may also be written from a previously read file (after manipulation) and wavelet attributes may not be present.

Variables:

filename: Name of the fits input file containing the wavelet transform data to be read.
tr: The initialised tr structure containing the data read from the input fits file.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_fits_error_check

private subroutine cswt_tr_io_fits_error_check (status, halt)
    integer, intent(inout) :: status
    logical, intent(in) :: halt
    ! Calls: ftgerr, ftgmsg
end subroutine cswt_tr_io_fits_error_check

Check if a fits error has occured and print error message. Halt program execution if halt flag is set.

Variables:

status: Fits integer status code.
halt: Logical to indicate whether to halt program execution if an error is detected.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_fits_exists

private subroutine cswt_tr_io_fits_exists (filename, status, exists)
    character (len=*), intent(in) :: filename
    integer, intent(inout) :: status
    logical, intent(out) :: exists
    ! Calls: cswt_tr_io_fits_error_check, ftclos, ftfiou, ftgiou, ftopen
end subroutine cswt_tr_io_fits_exists

Check if a fits file exists.

Variables:

filename: Name of fits file to check existence of.
status: Fits integer status code.
exists: Logical indicating whether the fits file already exists.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_fits_del

private subroutine cswt_tr_io_fits_del (filename, status)
    character (len=*), intent(in) :: filename
    integer, intent(inout) :: status
    ! Calls: ftcmsg, ftdelt, ftfiou, ftgiou, ftopen
end subroutine cswt_tr_io_fits_del

Delete a fits file.

Variables:

filename: Name of fits file to detele.
status: Fits integer status code.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_txt_dilation_write

public subroutine cswt_tr_io_txt_dilation_write (filename, tr, unit)
    character (len=*), intent(in) :: filename
    type (cswt_tr), intent(in) :: tr
    integer, intent(in) :: unit
    ! Calls: cswt_error
end subroutine cswt_tr_io_txt_dilation_write

Write tr structure dilation values to an output text file.

Notes:

Units stored in tr%dilation is always in radians. However units read and written to file may be in either radians or arcmins.

Variables:

filename: Name of file to write dilations to.
tr: Tr struture containing the dilations to be written.
unit: Integer specifying the units to write the dilation values in (either in arcmins or radians depending of the global CSWT_TR_DIL_UNIT_* specifier passed).

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_io_txt_dilation_read

public subroutine cswt_tr_io_txt_dilation_read (filename, dilation)
    character (len=*), intent(in) :: filename
    real (kind=s2_sp), intent(out), allocatable, dimension (:,:) :: dilation
    ! Calls: cswt_error
end subroutine cswt_tr_io_txt_dilation_read

Read dilation values from an input dilation text file.

Notes:

Units stored in tr%dilation is always in radians. However units read and written to file may be in either radians or arcmins.
The dilation array is allocated herein and must be freed by the calling routine.

Variables:

filename: Name of the input dilation file to be read.
dilation: Array of dilation values read from the input file.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_init

public function cswt_tr_get_init (tr) result (init)
    type (cswt_tr), intent(in) :: tr
    logical :: init
end function cswt_tr_get_init

Get init variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
init: Copy of object init variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_swav_mother

public function cswt_tr_get_swav_mother (tr) result (swav_mother)
    type (cswt_tr), intent(in) :: tr
    type (cswt_swav) :: swav_mother
    ! Calls: cswt_error
end function cswt_tr_get_swav_mother

Get copy of swav_mother variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
swav_mother: Copy of spherical mother wavelet extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_wcoeff

public subroutine cswt_tr_get_wcoeff (tr, wcoeff)
    type (cswt_tr), intent(in) :: tr
    real (kind=s2_sp), intent(out), dimension (1:,0:,0:,0:) :: wcoeff
    ! Calls: cswt_error
end subroutine cswt_tr_get_wcoeff

Get copy of the wavelet coefficients.

Notes:

Note that the wcoeff data array is indexed from one for the dilation dimensions, but from zero for the Euler angle dimensions.
The size of the input wcoeff array is not checked and must adhere to the dimensions of the tr%wcoeff array (and be index similarly as described by the previous note).

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
swav_mother: Copy of spherical mother wavelet extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_dilation

public subroutine cswt_tr_get_dilation (tr, dilation)
    type (cswt_tr), intent(in) :: tr
    real (kind=s2_sp), intent(out), dimension (:,:) :: dilation
    ! Calls: cswt_error
end subroutine cswt_tr_get_dilation

Get copy of dilation variable from the passed tr.

Notes:

The size of the input dilation array is not checked and must adhere to the dimensions of the tr%dilation array.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
dilation: Copy of object dilation variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_lmax

public function cswt_tr_get_lmax (tr) result (lmax)
    type (cswt_tr), intent(in) :: tr
    integer :: lmax
    ! Calls: cswt_error
end function cswt_tr_get_lmax

Get copy of lmax variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
lmax: Copy of object lmax variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_mmax

public function cswt_tr_get_mmax (tr) result (mmax)
    type (cswt_tr), intent(in) :: tr
    integer :: mmax
    ! Calls: cswt_error
end function cswt_tr_get_mmax

Get copy of mmax variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
mmax: Copy of object mmax variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_n_alpha

public function cswt_tr_get_n_alpha (tr) result (n_alpha)
    type (cswt_tr), intent(in) :: tr
    integer :: n_alpha
    ! Calls: cswt_error
end function cswt_tr_get_n_alpha

Get copy of n_alpha variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
n_alpha: Copy of object n_alpha variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_n_beta

public function cswt_tr_get_n_beta (tr) result (n_beta)
    type (cswt_tr), intent(in) :: tr
    integer :: n_beta
    ! Calls: cswt_error
end function cswt_tr_get_n_beta

Get copy of n_beta variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
n_beta: Copy of object n_beta variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_n_gamma

public function cswt_tr_get_n_gamma (tr) result (n_gamma)
    type (cswt_tr), intent(in) :: tr
    integer :: n_gamma
    ! Calls: cswt_error
end function cswt_tr_get_n_gamma

Get copy of n_gamma variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
n_gamma: Copy of object n_gamma variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_n_dilation

public function cswt_tr_get_n_dilation (tr) result (n_dilation)
    type (cswt_tr), intent(in) :: tr
    integer :: n_dilation
    ! Calls: cswt_error
end function cswt_tr_get_n_dilation

Get copy of n_dilation variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
n_dilaiton: Copy of object n_dilation variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_wcoeff_status

public function cswt_tr_get_wcoeff_status (tr) result (wcoeff_status)
    type (cswt_tr), intent(in) :: tr
    logical :: wcoeff_status
    ! Calls: cswt_error
end function cswt_tr_get_wcoeff_status

Get copy of wcoeff_status variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
wcoeff_status: Copy of object wcoeff_status variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004

cswt_tr_get_swav_mother_status

public function cswt_tr_get_swav_mother_status (tr) result (swav_mother_status)
    type (cswt_tr), intent(in) :: tr
    logical :: swav_mother_status
    ! Calls: cswt_error
end function cswt_tr_get_swav_mother_status

Get copy of swav_mother_status variable from the passed tr.

Variables:

tr: Spherical wavelet transform object to get copy of variable of.
swav_mother_status: Copy of object swav_mother_status variable extracted.

Author: J. D. McEwen

Version: 0.1 - November 2004