import numpy as np
import mrcfile
import math
from statistics import mean
import sys
from numpy.fft import fftn, ifftn, ifftshift
import argparse
import os
from write_mrc import write_mrc
from scipy.spatial import cKDTree
from loc_spiral import open_map_mask, cart_2_sph, spiral_filter_vect, spiral_filter

#   locBSharpen - method to local sharpening by B-factor correction

#   Method to improve the interpretability of cryo-EM maps improving the 
#   contrast and interpretability at high-resolutions specially for maps
#   affected by heterogeneous resolutions and SNRs

#   INPUT PARAMETERS:
#   map: cryo-EM map to be processed. This maps should be the average of
#   the half maps obtained after refinment or the final map without any
#   filtering, postprocessing (masking, b-factor correction).

#   mask: Tight binary map that differenciates between the macromolecule (1)
#   and the background (0).

#   pixel_size: sampling rate of the cryo-EM map in Angstroms.

#   minRes: minimun resolution of the sweeping resolution range in Angstroms.
#   Usually a value of 25-35 Angstroms is good.

#   maxRes: maximum resolution of the sweeping resolution range. Provide here 
#   the obtained FSC global resolution value in Angstroms.

#   threshold: noise signficance in the comparision with noise. Excellent  
#   values are 0.9 or 0.95.

#   bw: bandwith of the bandpass filter applied in Fourier Space. For high
#   quality maps with global resolution (<4 Angstroms) and showing homogeneous 
#   distribution of resolutions a value of 3.5-5 is good. For maps presenting 
#   heterogeneous distributions of resolutions, or worse global resolution, 
#   values between 6-8 provide excellents results.

#   OUTPUT PARAMETERS:
#   M: Improved output map
#   W: Sum of all amplitude maps during the resolution sweeping that can be
#   consider as an quality map representing the amount of signal

#   If this method was useful to you, please cite the following paper:
#   <NAME>, <NAME>, <NAME>, <NAME>, <NAME>,
#   <NAME>, <NAME>, <NAME>, <NAME>, "Local computational methods to 
#   improve the interpretability and analysis of cryo-EM maps" BioRxiv, 
#   doi: https://doi.org/10.1101/2020.05.11.088013 (2020)

#   <NAME>, 01/06/20
#   <EMAIL>
#   Copyright 2020, Universidad Complutense de Madrid
#   $ Revision: 1.0.0.0
#   $ Date: 01/06/20

# Copyright 2020 <NAME> @UCM

# Redistribution and use in source and binary forms, with or without modification, 
# are permitted provided that the following conditions are met:

# 1. Redistributions of source code must retain the above copyright notice, 
# this list of conditions and the following disclaimer.

# 2. Redistributions in binary form must reproduce the above copyright notice, 
# this list of conditions and the following disclaimer in the documentation 
# and/or other materials provided with the distribution.

# 3. Neither the name of the copyright holder nor the names of its contributors 
# may be used to endorse or promote products derived from this software without 
# specific prior written permission.

# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 
# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


# The eco option is an aproximation for the spiral phase transform that
# works fine in all the cases tested and saves a lot of computational
# resources. If case you use a computer with limited RAM memory (<32Gb and
# box sizes smaller than 400x400x400px) use the eco option. If not, it is
# more appropiate to use eco = false.


def loc_b_sharpen(map, mask, pixel_size, min_res, max_res, threshold, bandwidth, eco = True):

    if not np.any(mask == 0):
        print('Your mask contains no zero points. Please enter a threshold for binarization')
        binary_thresh = float(input('Threshold: '))
        bin_mask = mask.copy()
        bin_mask[bin_mask <= binary_thresh] = 0
        mask = bin_mask

    M = np.multiply(map, 0)
    W = np.multiply(map, 0)

    [n, n1, n2] = np.shape(map)
    max_freq = math.floor(pixel_size * n / max_res)
    min_freq = math.floor(pixel_size * n / min_res)

    # We create necessary staff for the bank of filters once:
    [NR, NC, NZ] = np.shape(map)
    [u, v, w] = np.meshgrid(range(1, NC + 1), range(1, NR + 1), range(1, NZ + 1))

    # Temporal Variables
    ## Rich - calculating the mean of all entries in a grid is just shape/2
    u0 = math.floor(n/2)
    v0 = math.floor(n/2)
    w0 = math.floor(n/2)

    u = u - u0
    v = v - v0
    w = w - w0

    # Circular mask to be used in the quantile evaluation
    mask_c = np.power((np.power(u, 2) + np.power(v, 2) + np.power(w, 2)), 0.5)
    mask_c = mask_c < u0

    [theta, phi, radial] = cart_2_sph(u, v, w)
    
    # Gaussian DC filter with sigma = 0.05
    Hr = 1-np.exp(-(np.power(u, 2) + np.power(v, 2) + np.power(w, 2))/(2*0.05**2))

    if eco:
        Hs = spiral_filter([u, v, w])
    else:
        [Hsx, Hsy, Hsz] = spiral_filter_vect([u, v, w])

    C = fftn(map)
    num = abs(max_freq - min_freq)
    freqs = np.linspace(min_freq, max_freq, num)

    # Sharpening starts at 15--10 A
    sharpening_freq = math.floor(pixel_size * n / 15)

    # In matlab, M' is the complex conjugate of M, but sharpening freq should be real,
    # so I think the transpose should work fine.
    #
    # We need to make a KDTree first, then look up the points
    # From https://stackoverflow.com/questions/66494042/dsearchn-equivalent-in-python

    kdt = cKDTree(np.matrix(freqs).T)
    inx_sharpening_freq = kdt.query(sharpening_freq)[1]

    # Amplitude map at the sharpening resolution (15A) that is going to be
    # impose for higher resolutions

    aMap = np.multiply(0, map)

    for i in range(0, len(freqs)):
        if i%10 == 0:
            print(f'On iteration {i} out of {len(freqs)}')
        # actual filter

        # Annular filter with radius=R and sigma=S
        H = np.exp(-np.divide(
            np.power((radial - freqs[i]), 2),
            2*bandwidth**2))

        H = np.multiply(Hr, H)
        H = H/sum(H.flatten())
        CH = np.multiply(C, ifftshift(H))
        ch = np.real(ifftn(CH))

        if eco:
            CH = np.multiply(CH, ifftshift(Hs))
            s = abs(np.conj(ifftn(CH)))
        else:
            # using vectorized spiral phase transform
            CHx = ifftn(np.multiply(CH, ifftshift(Hsx)))
            CHx = np.multiply(CHx, np.conj(CHx))

            CHy = ifftn(np.multiply(CH, ifftshift(Hsy)))
            CHy = np.multiply(CHy, np.conj(CHy))

            CHz = ifftn(np.multiply(CH, ifftshift(Hsz)))
            CHz = np.multiply(CHz, np.conj(CHz))

            s = np.power((CHx + CHy + CHz), 0.5)

        # normalized igram
        ## Rich - original code doesn't coerce to real here, but
        ## it seems like the imaginary part is 0
        if not np.all(np.isreal(s)) and not np.all(np.isreal(ch)):
            print('s or ch are non-real')
            sys.exit(2)

        cn = np.cos(np.arctan2(np.real(s), np.real(ch)))

        # modulation
        m = abs(ch + 1j*s)

        # I'm assuming what ~mask does here
        q = np.quantile(m[np.multiply(mask == 0, mask_c) > 0.5], threshold)
        Cref = np.divide(m, (m + q))

        m = m >= q
        m = m.astype(np.double)

        if i < inx_sharpening_freq:
            M = M + np.multiply(Cref, np.multiply(cn, m))
            W = W + np.multiply(Cref, m)
        else:
            if i == inx_sharpening_freq:
                aMap = m
            M = M + np.multiply(Cref, np.multiply(cn, aMap))
            W = W + np.multiply(Cref, aMap)

    M = M/max_freq
    W = W/max_freq

    M = M - mean(M[np.multiply(mask == 0, mask_c) > 0.5])
    M = np.divide(M, np.std(M[np.multiply(mask == 0, mask_c) > 0.5]))

    W = W - mean(W[np.multiply(mask == 0, mask_c) > 0.5])
    W = np.divide(W, np.std(W[np.multiply(mask == 0, mask_c) > 0.5]))

    # Impose positivity
    M[M<0] = 0

    return (M, W)

parser = argparse.ArgumentParser()
parser.add_argument(
    'map',
    help = 'cryoEM map to be processed. Should not be pre-processed.',
)
parser.add_argument(
    'mask',
    help = 'Tight binary mask that differenciates between macromolecule and background',
)
parser.add_argument(
    '-p', '--pixel-size',
    help = 'Pixel size in Angstroms of map',
    type = float
)
parser.add_argument(
    '-l', '--min-res',
    help = 'Minimum resolution. Try a value from 25--35 A',
    type = float
)
parser.add_argument(
    '-m', '--max-res',
    help = 'Maximum resolution provided by the GSFSC in Angstroms',
    type = float
)
parser.add_argument(
    '-t', '--threshold',
    help = 'Noise significance in comparison with noise. Try 0.9--0.95',
    type = float
)
parser.add_argument(
    '-b', '--bandwidth',
    help = 'Bandwidth of bandpass filter applied in Fourier Space. For high-quality maps, try 3.5--5. For more heterogeneous maps, try 6--8',
    type = float
)
parser.add_argument(
    '-e', '--eco',
    help = 'Use a faster approximation for the spiral filter',
    action = 'store_true'
)
parser.add_argument(
    '-o', '--output',
    help = 'Where to put sharpened map',
    default = 'sharpened.mrc'
)

def main(args):
    map_path = os.path.abspath(args.map)
    mask_path = os.path.abspath(args.mask)
    output = os.path.abspath(args.output)
    map, mask = open_map_mask(map_path, mask_path)

    sharp_map, W = loc_b_sharpen(
        map, mask, args.pixel_size,
        args.min_res, args.max_res,
        args.threshold, args.bandwidth, 
        args.eco)

    if output[-4:]!= '.mrc':
        output = output + '.mrc'

    write_mrc(sharp_map, args.pixel_size, output)

if __name__ == '__main__':
    args = parser.parse_args()
    main(args)
import numpy as np
import mrcfile
import math
from statistics import mean
import sys
from numpy.fft import fftn, ifftn, ifftshift
import argparse
import os
from write_mrc import write_mrc
from scipy.spatial import cKDTree
from loc_spiral import open_map_mask, cart_2_sph, spiral_filter_vect, spiral_filter

#   locBSharpen - method to local sharpening by B-factor correction

#   Method to improve the interpretability of cryo-EM maps improving the 
#   contrast and interpretability at high-resolutions specially for maps
#   affected by heterogeneous resolutions and SNRs

#   INPUT PARAMETERS:
#   map: cryo-EM map to be processed. This maps should be the average of
#   the half maps obtained after refinment or the final map without any
#   filtering, postprocessing (masking, b-factor correction).

#   mask: Tight binary map that differenciates between the macromolecule (1)
#   and the background (0).

#   pixel_size: sampling rate of the cryo-EM map in Angstroms.

#   minRes: minimun resolution of the sweeping resolution range in Angstroms.
#   Usually a value of 25-35 Angstroms is good.

#   maxRes: maximum resolution of the sweeping resolution range. Provide here 
#   the obtained FSC global resolution value in Angstroms.

#   threshold: noise signficance in the comparision with noise. Excellent  
#   values are 0.9 or 0.95.

#   bw: bandwith of the bandpass filter applied in Fourier Space. For high
#   quality maps with global resolution (<4 Angstroms) and showing homogeneous 
#   distribution of resolutions a value of 3.5-5 is good. For maps presenting 
#   heterogeneous distributions of resolutions, or worse global resolution, 
#   values between 6-8 provide excellents results.

#   OUTPUT PARAMETERS:
#   M: Improved output map
#   W: Sum of all amplitude maps during the resolution sweeping that can be
#   consider as an quality map representing the amount of signal

#   If this method was useful to you, please cite the following paper:
#   <NAME>, <NAME>, <NAME>, <NAME>, <NAME>,
#   <NAME>, <NAME>, <NAME>, <NAME>, "Local computational methods to 
#   improve the interpretability and analysis of cryo-EM maps" BioRxiv, 
#   doi: https://doi.org/10.1101/2020.05.11.088013 (2020)

#   <NAME>, 01/06/20
#   <EMAIL>
#   Copyright 2020, Universidad Complutense de Madrid
#   $ Revision: 1.0.0.0
#   $ Date: 01/06/20

# Copyright 2020 <NAME> @UCM

# Redistribution and use in source and binary forms, with or without modification, 
# are permitted provided that the following conditions are met:

# 1. Redistributions of source code must retain the above copyright notice, 
# this list of conditions and the following disclaimer.

# 2. Redistributions in binary form must reproduce the above copyright notice, 
# this list of conditions and the following disclaimer in the documentation 
# and/or other materials provided with the distribution.

# 3. Neither the name of the copyright holder nor the names of its contributors 
# may be used to endorse or promote products derived from this software without 
# specific prior written permission.

# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 
# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


# The eco option is an aproximation for the spiral phase transform that
# works fine in all the cases tested and saves a lot of computational
# resources. If case you use a computer with limited RAM memory (<32Gb and
# box sizes smaller than 400x400x400px) use the eco option. If not, it is
# more appropiate to use eco = false.


def loc_b_sharpen(map, mask, pixel_size, min_res, max_res, threshold, bandwidth, eco = True):

    if not np.any(mask == 0):
        print('Your mask contains no zero points. Please enter a threshold for binarization')
        binary_thresh = float(input('Threshold: '))
        bin_mask = mask.copy()
        bin_mask[bin_mask <= binary_thresh] = 0
        mask = bin_mask

    M = np.multiply(map, 0)
    W = np.multiply(map, 0)

    [n, n1, n2] = np.shape(map)
    max_freq = math.floor(pixel_size * n / max_res)
    min_freq = math.floor(pixel_size * n / min_res)

    # We create necessary staff for the bank of filters once:
    [NR, NC, NZ] = np.shape(map)
    [u, v, w] = np.meshgrid(range(1, NC + 1), range(1, NR + 1), range(1, NZ + 1))

    # Temporal Variables
    ## Rich - calculating the mean of all entries in a grid is just shape/2
    u0 = math.floor(n/2)
    v0 = math.floor(n/2)
    w0 = math.floor(n/2)

    u = u - u0
    v = v - v0
    w = w - w0

    # Circular mask to be used in the quantile evaluation
    mask_c = np.power((np.power(u, 2) + np.power(v, 2) + np.power(w, 2)), 0.5)
    mask_c = mask_c < u0

    [theta, phi, radial] = cart_2_sph(u, v, w)
    
    # Gaussian DC filter with sigma = 0.05
    Hr = 1-np.exp(-(np.power(u, 2) + np.power(v, 2) + np.power(w, 2))/(2*0.05**2))

    if eco:
        Hs = spiral_filter([u, v, w])
    else:
        [Hsx, Hsy, Hsz] = spiral_filter_vect([u, v, w])

    C = fftn(map)
    num = abs(max_freq - min_freq)
    freqs = np.linspace(min_freq, max_freq, num)

    # Sharpening starts at 15--10 A
    sharpening_freq = math.floor(pixel_size * n / 15)

    # In matlab, M' is the complex conjugate of M, but sharpening freq should be real,
    # so I think the transpose should work fine.
    #
    # We need to make a KDTree first, then look up the points
    # From https://stackoverflow.com/questions/66494042/dsearchn-equivalent-in-python

    kdt = cKDTree(np.matrix(freqs).T)
    inx_sharpening_freq = kdt.query(sharpening_freq)[1]

    # Amplitude map at the sharpening resolution (15A) that is going to be
    # impose for higher resolutions

    aMap = np.multiply(0, map)

    for i in range(0, len(freqs)):
        if i%10 == 0:
            print(f'On iteration {i} out of {len(freqs)}')
        # actual filter

        # Annular filter with radius=R and sigma=S
        H = np.exp(-np.divide(
            np.power((radial - freqs[i]), 2),
            2*bandwidth**2))

        H = np.multiply(Hr, H)
        H = H/sum(H.flatten())
        CH = np.multiply(C, ifftshift(H))
        ch = np.real(ifftn(CH))

        if eco:
            CH = np.multiply(CH, ifftshift(Hs))
            s = abs(np.conj(ifftn(CH)))
        else:
            # using vectorized spiral phase transform
            CHx = ifftn(np.multiply(CH, ifftshift(Hsx)))
            CHx = np.multiply(CHx, np.conj(CHx))

            CHy = ifftn(np.multiply(CH, ifftshift(Hsy)))
            CHy = np.multiply(CHy, np.conj(CHy))

            CHz = ifftn(np.multiply(CH, ifftshift(Hsz)))
            CHz = np.multiply(CHz, np.conj(CHz))

            s = np.power((CHx + CHy + CHz), 0.5)

        # normalized igram
        ## Rich - original code doesn't coerce to real here, but
        ## it seems like the imaginary part is 0
        if not np.all(np.isreal(s)) and not np.all(np.isreal(ch)):
            print('s or ch are non-real')
            sys.exit(2)

        cn = np.cos(np.arctan2(np.real(s), np.real(ch)))

        # modulation
        m = abs(ch + 1j*s)

        # I'm assuming what ~mask does here
        q = np.quantile(m[np.multiply(mask == 0, mask_c) > 0.5], threshold)
        Cref = np.divide(m, (m + q))

        m = m >= q
        m = m.astype(np.double)

        if i < inx_sharpening_freq:
            M = M + np.multiply(Cref, np.multiply(cn, m))
            W = W + np.multiply(Cref, m)
        else:
            if i == inx_sharpening_freq:
                aMap = m
            M = M + np.multiply(Cref, np.multiply(cn, aMap))
            W = W + np.multiply(Cref, aMap)

    M = M/max_freq
    W = W/max_freq

    M = M - mean(M[np.multiply(mask == 0, mask_c) > 0.5])
    M = np.divide(M, np.std(M[np.multiply(mask == 0, mask_c) > 0.5]))

    W = W - mean(W[np.multiply(mask == 0, mask_c) > 0.5])
    W = np.divide(W, np.std(W[np.multiply(mask == 0, mask_c) > 0.5]))

    # Impose positivity
    M[M<0] = 0

    return (M, W)

parser = argparse.ArgumentParser()
parser.add_argument(
    'map',
    help = 'cryoEM map to be processed. Should not be pre-processed.',
)
parser.add_argument(
    'mask',
    help = 'Tight binary mask that differenciates between macromolecule and background',
)
parser.add_argument(
    '-p', '--pixel-size',
    help = 'Pixel size in Angstroms of map',
    type = float
)
parser.add_argument(
    '-l', '--min-res',
    help = 'Minimum resolution. Try a value from 25--35 A',
    type = float
)
parser.add_argument(
    '-m', '--max-res',
    help = 'Maximum resolution provided by the GSFSC in Angstroms',
    type = float
)
parser.add_argument(
    '-t', '--threshold',
    help = 'Noise significance in comparison with noise. Try 0.9--0.95',
    type = float
)
parser.add_argument(
    '-b', '--bandwidth',
    help = 'Bandwidth of bandpass filter applied in Fourier Space. For high-quality maps, try 3.5--5. For more heterogeneous maps, try 6--8',
    type = float
)
parser.add_argument(
    '-e', '--eco',
    help = 'Use a faster approximation for the spiral filter',
    action = 'store_true'
)
parser.add_argument(
    '-o', '--output',
    help = 'Where to put sharpened map',
    default = 'sharpened.mrc'
)

def main(args):
    map_path = os.path.abspath(args.map)
    mask_path = os.path.abspath(args.mask)
    output = os.path.abspath(args.output)
    map, mask = open_map_mask(map_path, mask_path)

    sharp_map, W = loc_b_sharpen(
        map, mask, args.pixel_size,
        args.min_res, args.max_res,
        args.threshold, args.bandwidth, 
        args.eco)

    if output[-4:]!= '.mrc':
        output = output + '.mrc'

    write_mrc(sharp_map, args.pixel_size, output)

if __name__ == '__main__':
    args = parser.parse_args()
    main(args)
#   locBSharpen - method to local sharpening by B-factor correction #   Method to improve the interpretability of cryo-EM maps improving the #   contrast and interpretability at high-resolutions specially for maps #   affected by heterogeneous resolutions and SNRs #   INPUT PARAMETERS: #   map: cryo-EM map to be processed. This maps should be the average of #   the half maps obtained after refinment or the final map without any #   filtering, postprocessing (masking, b-factor correction). #   mask: Tight binary map that differenciates between the macromolecule (1) #   and the background (0). #   pixel_size: sampling rate of the cryo-EM map in Angstroms. #   minRes: minimun resolution of the sweeping resolution range in Angstroms. #   Usually a value of 25-35 Angstroms is good. #   maxRes: maximum resolution of the sweeping resolution range. Provide here #   the obtained FSC global resolution value in Angstroms. #   threshold: noise signficance in the comparision with noise. Excellent #   values are 0.9 or 0.95. #   bw: bandwith of the bandpass filter applied in Fourier Space. For high #   quality maps with global resolution (<4 Angstroms) and showing homogeneous #   distribution of resolutions a value of 3.5-5 is good. For maps presenting #   heterogeneous distributions of resolutions, or worse global resolution, #   values between 6-8 provide excellents results. #   OUTPUT PARAMETERS: #   M: Improved output map #   W: Sum of all amplitude maps during the resolution sweeping that can be #   consider as an quality map representing the amount of signal #   If this method was useful to you, please cite the following paper: #   <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, #   <NAME>, <NAME>, <NAME>, <NAME>, "Local computational methods to #   improve the interpretability and analysis of cryo-EM maps" BioRxiv, #   doi: https://doi.org/10.1101/2020.05.11.088013 (2020) #   <NAME>, 01/06/20 #   <EMAIL> #   Copyright 2020, Universidad Complutense de Madrid #   $ Revision: 1.0.0.0 #   $ Date: 01/06/20 # Copyright 2020 <NAME> @UCM # Redistribution and use in source and binary forms, with or without modification, # are permitted provided that the following conditions are met: # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. # 3. Neither the name of the copyright holder nor the names of its contributors # may be used to endorse or promote products derived from this software without # specific prior written permission. # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR # SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE # USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. # The eco option is an aproximation for the spiral phase transform that # works fine in all the cases tested and saves a lot of computational # resources. If case you use a computer with limited RAM memory (<32Gb and # box sizes smaller than 400x400x400px) use the eco option. If not, it is # more appropiate to use eco = false. # We create necessary staff for the bank of filters once: # Temporal Variables ## Rich - calculating the mean of all entries in a grid is just shape/2 # Circular mask to be used in the quantile evaluation # Gaussian DC filter with sigma = 0.05 # Sharpening starts at 15--10 A # In matlab, M' is the complex conjugate of M, but sharpening freq should be real, # so I think the transpose should work fine. # # We need to make a KDTree first, then look up the points # From https://stackoverflow.com/questions/66494042/dsearchn-equivalent-in-python # Amplitude map at the sharpening resolution (15A) that is going to be # impose for higher resolutions # actual filter # Annular filter with radius=R and sigma=S # using vectorized spiral phase transform # normalized igram ## Rich - original code doesn't coerce to real here, but ## it seems like the imaginary part is 0 # modulation # I'm assuming what ~mask does here # Impose positivity