Metadata-Version: 2.1
Name: alpha-transform
Version: 0.0.1
Summary: Adaptive alpha shearlet transform for manifold-valued data
Home-page: https://github.com/dedale-fet/alpha-transform
Author: Felix Voigtlaender
Author-email: felix@voigtlaender.xyz
License: MIT
Description: # Adaptive transform for manifold-valued data
        
        This package contains all software developed as part of Workpackage 2.1 (Adaptive transform for manifold-valued data)
        of the DEDALE project. For an in-depth description of this workpackage, we refer to the associated technical report.
        
        The software consists of three main parts:
        
        1. An implementation of the (bandlimited) **α-shearlet transform** (in [AlphaTransform.py](../AlphaTransform.py),
           in three versions:
           * A *fully sampled* (non-decimated), translation invariant, fast, but memory-*consuming* implementation
           * A fully sampled, translation invariant, slightly slower, but memory-*efficient* implementation
           * A *subsampled* (decimated), *not* translation invariant, but fast and memory-efficient implementation.
        2. Implementations (in [AdaptiveAlpha.py](../AdaptiveAlpha.py)) of three criteria that can be used to adaptively choose
           the value of α, namely:
           * The [*asymptotic approximation rate*](#aar) (**AAR**),
           * the [*mean approximation error*](#mae) (**MAE**),
           * the [*thresholding denoising performance*](#tdp) (**TDP**).
        3. A [chart-based implementation](#sphere-valued) (in [SphereTransform.py](../SphereTransform.py)) of the α-shearlet transform for
           functions defined on the sphere.
        
        
        In the following, we provide brief explanations and hands-on experiments for all of these aspects. The following
        table of content can be used for easy navigation:
        
        ## Table of Contents
        1. [The α-shearlet transform: A crash course](#shearlet-trafo-crash-course)
        
        2. [Adaptively choosing α](#adaptivity)
        
             a) [Asymptotic approximation rate (AAR)](#aar)
        
             b) [Mean approximation error (MAE)](#mae)
        
             c) [Thresholding denoising performance (TDP)](#tdp)
        
        3. [The α-shearlet transform for functions defined on the sphere](#sphere-valued)
        
        
        ## <a name="shearlet-trafo-crash-course"></a>The α-shearlet transform: A crash course
        
        The following demonstrates a very simple use case of the α-shearlet transform: We compute the transform of an example image,
        threshold the coefficients, reconstruct and compute the error. The code is longer than strictly necessary, since
        along the way give a demonstration of the general usage of the transform.
        
        ```python
        >>> # Importing necessary packages
        >>> from AlphaTransform import AlphaShearletTransform as AST
        >>> import numpy as np
        >>> import matplotlib.pyplot as plt
        >>> from scipy import misc
        
        >>> im = misc.face(gray=True)
        >>> im.shape
        (768, 1024)
        
        >>> # Setting up the transform
        >>> trafo = AST(im.shape[1], im.shape[0], [0.5]*3) # 1
        Precomputing shearlet system: 100%|███████████████████████████████████████| 52/52 [00:05<00:00,  8.83it/s]
        
        >>> # Computing and understanding the α-shearlet coefficients
        >>> coeff = trafo.transform(im) # 2
        >>> coeff.shape # 3
        (53, 768, 1024)
        >>> trafo.indices # 4
        [-1,
         (0, -1, 'r'), (0, 0, 'r'), (0, 1, 'r'),
         (0, 1, 't'), (0, 0, 't'), (0, -1, 't'),
         (0, -1, 'l'), (0, 0, 'l'), (0, 1, 'l'),
         (0, 1, 'b'), (0, 0, 'b'), (0, -1, 'b'),
         (1, -2, 'r'), (1, -1, 'r'), (1, 0, 'r'), (1, 1, 'r'), (1, 2, 'r'),
         (1, 2, 't'), (1, 1, 't'), (1, 0, 't'), (1, -1, 't'), (1, -2, 't'),
         (1, -2, 'l'), (1, -1, 'l'), (1, 0, 'l'), (1, 1, 'l'), (1, 2, 'l'),
         (1, 2, 'b'), (1, 1, 'b'), (1, 0, 'b'), (1, -1, 'b'), (1, -2, 'b'),
         (2, -2, 'r'), (2, -1, 'r'), (2, 0, 'r'), (2, 1, 'r'), (2, 2, 'r'),
         (2, 2, 't'), (2, 1, 't'), (2, 0, 't'), (2, -1, 't'), (2, -2, 't'),
         (2, -2, 'l'), (2, -1, 'l'), (2, 0, 'l'), (2, 1, 'l'), (2, 2, 'l'),
         (2, 2, 'b'), (2, 1, 'b'), (2, 0, 'b'), (2, -1, 'b'), (2, -2, 'b')]
        
        >>> # Thresholding the coefficients and reconstructing
        >>> np.max(np.abs(coeff)) # 5
        2041.1017181588547
        >>> np.sum(np.abs(coeff) > 200) / coeff.size # 6
        0.020679905729473761
        >>> thresh_coeff = coeff * (np.abs(coeff) > 200) # 7
        >>> recon = trafo.inverse_transform(thresh_coeff, real=True) # 8
        >>> np.linalg.norm(im - recon) / np.linalg.norm(im) # 9
        0.13912540983541383
        
        >>> plt.imshow(recon, cmap='gray')
        <matplotlib.image.AxesImage object at 0x2b0f568c25c0>
        >>> plt.show()
        ```
        
        ### Importing necessary packages
        
        The first few (unnumbered) lines import relevant packages and load (a gray scale version of) the following
        test image, which has a resolution of 1024 x 768 pixels:
        <p align="center">
        <img src="http://www.scipy-lectures.org/_images/face.png" width="512"/>
        </p>
        
        
        ### Setting up the transform
        
        Then, in line 1, we create an instance `trafo` of the class `AlphaShearletTransform` (which is
        called `AST` above for brevity).
        During construction of this object, the necessary shearlet filters are precomputed. This may take some time
        (5 seconds in the example above), but speeds up later computations.
        
        The three parameters passed to the constructor of `trafo` require some explanation:
        * The first parameter is the *width* of the images which can be analyzed using the `trafo` object.
        * Similarly, the second parameter determines the *height*.
        * The third parameter is of the form `[alpha] * N`, where `N` determines the number of scales of the transform and `alpha`
          determines the value of α.
          
          The reason for this notation is that in principle, one can choose a
          different value of α on each scale. Since `[0.5] * 3 = [0.5, 0.5, 0.5]`, one can use this notation to obtain a
          system using a single value of α across all scales.
        
        All in all, we see that line 1 creates a shearlet system (i.e., α = 0.5) with 3 scale (plus a low-pass)
        for images of dimension 1024 x 768.
        
        ### Computing and understanding the α-shearlet coefficients
        
        Line 2 shows that the α-shearlet coefficients of the image `im` can be readily computed using the `transform` method of `trafo`.
        As seen in line 3, this results in an array of size `53`, where each of the elements of the array is an array (an image) of
        size 1024 x 768, i.e., of the same size as the input image.
        
        To help understand the meaning of each of these **coefficient images** `coeff[i]`, for `i = 0, ..., 52`, the output of line 4
        is helpful: Associated to each coefficient image `coeff[i]`, there is an **index** `trafo.indices[i]` which encodes the *meaning*
        of the coefficient image, i.e., the shearlet used to compute it.
        
        The special index `-1` stands for the **low-pass part**. All other indices are of the form `(j, k, c)`, where
        * `j` encodes the **scale**. In the present case, `j` ranges from `0` to `2`, since we have 3 scales.
        * `k` encodes the **amount of shearing**, ranging from -⌈2<sup>j(1 - α)</sup>⌉ to ⌈2<sup>j(1 - α)</sup>⌉ on scale ``j``.
        * `c` encodes the **cone** to which the shearlet belongs (in the Fourier domain). Precisely, we have the following correspondence
          between the value of `c` and the corresponding frequency cones:
        
          | value of `c`   | `'r'` | `'t'` | `'l'` | `'b'`  |
          | -------------- |:-----:|:-----:|:-----:|:------:|
          | Frequency cone | right |  top  | left  | bottom |
        
        
        Note that if we divide the frequency plane into four cones, such that each shearlet has a real-valued Fourier transform
        which is supported in one of these cones, then the shearlets themselves (i.e., in space) *can not* be real-valued.
        Hence, if real-valued shearlets are desired, one can pass the constructor of the class `AlphaShearletTransform` the
        additional argument `real=True`. In this case, the frequency plane is split into a horizontal (encoded by `'h'`) and
        a vertical (encoded by `'v'`) cone, as is indicated by the following example:
        ```python
        >>> trafo_real = AST(im.shape[1], im.shape[0], [0.5]*3, real=True)
        Precomputing shearlet system: 100%|██████████████████████████████████████| 26/26 [00:03<00:00,  5.62it/s]
        >>> trafo_real.indices
        [-1,
         (0, -1, 'h'), (0, 0, 'h'), (0, 1, 'h'),
         (0, 1, 'v'), (0, 0, 'v'), (0, -1, 'v'),
         (1, -2, 'h'), (1, -1, 'h'), (1, 0, 'h'), (1, 1, 'h'), (1, 2, 'h'),
         (1, 2, 'v'), (1, 1, 'v'), (1, 0, 'v'), (1, -1, 'v'), (1, -2, 'v'),
         (2, -2, 'h'), (2, -1, 'h'), (2, 0, 'h'), (2, 1, 'h'), (2, 2, 'h'),
         (2, 2, 'v'), (2, 1, 'v'), (2, 0, 'v'), (2, -1, 'v'), (2, -2, 'v')]
        ```
        
        ### Thresholding the coefficients & reconstructing
        
        When called without further parameters, the method `trafo.transform` computes a *normalized* transform,
        so that effectively all shearlets are normalized to have L² norm 1. With this normalization, line 5 shows
        that the largest coefficient has size about 2041. We now (arbitrarily) pick a threshold of 200 and see (in line 6)
        that only about 2% of the coefficients are larger than this threshold. Next, we set all coefficients which are smaller
        (in absolute value) than 200 to zero and save the resulting thresholded coefficients as `thresh_coeff`, in line 7.
        
        In line 8, we then use the method `inverse_transform` of the `trafo` object to compute the inverse transform.
        Since we know that the original image was real-valued, we pass the additional argument `real=True`.
        This has the same effect as reconstructing without this additional argument and then taking the real part.
        Line 9 shows that the relative error is about 13.9%.
        
        Finally, the last two lines display the reconstructed image.
        
        ### Changes for the subsampled transform
        
        Above, we showed how our implementation of the *fully sampled* α-shearlet transform can be used to compute the α-shearlet
        transform of an image and reconstruct (with thresholded coefficients). For the *subsampled* transform, this can be
        done very similarly; the main difference one has to keep in mind is that for the fully sampled transform,
        one obtains an array of "coefficient images" which are all of the same size. In contrast, due to the subsampling,
        the "coefficient images" for the subsampled transform are all of different sizes:
        
        ```python
        >>> from AlphaTransform import AlphaShearletTransform as AST
        >>> import matplotlib.pyplot as plt
        >>> import numpy as np
        >>> from scipy import misc
        
        >>> im = misc.face(gray=True)
        >>> trafo = AST(im.shape[1], im.shape[0], [0.5]*3, subsampled=True) # 1
        Precomputing shearlets: 100%|████████████████████████████████████████████| 52/52 [00:00<00:00, 69.87it/s]
        >>> coeff = trafo.transform(im) # 2
        
        >>> type(coeff).__name__ # 3
        'list'
        >>> [c.shape for c in coeff] # 4
        [(129, 129),
         (364, 161), (364, 161), (364, 161),
         (97, 257), (97, 257), (97, 257),
         (364, 161), (364, 161), (364, 161),
         (97, 257), (97, 257), (97, 257),
         (514, 321), (514, 321), (514, 321), (514, 321), (514, 321),
         (193, 364), (193, 364), (193, 364), (193, 364), (193, 364),
         (514, 321), (514, 321), (514, 321), (514, 321), (514, 321),
         (193, 364), (193, 364), (193, 364), (193, 364), (193, 364),
         (727, 641), (727, 641), (727, 641), (727, 641), (727, 641),
         (385, 513), (385, 513), (385, 513), (385, 513), (385, 513),
         (727, 641), (727, 641), (727, 641), (727, 641), (727, 641),
         (385, 513), (385, 513), (385, 513), (385, 513), (385, 513)]
        
        >>> np.max([np.max(np.abs(c)) for c in coeff]) # 5
        2031.0471969998314
        >>> np.sum([np.sum(np.abs(c) > 200) for c in coeff]) # 6
        22357
        >>> np.sum([np.sum(np.abs(c) > 200) for c in coeff]) / np.sum([c.size for c in coeff]) # 7
        0.0023520267754635542
        >>> thresh_coeff = [c * (np.abs(c) > 200) for c in coeff] # 8
        >>> recon = trafo.inverse_transform(thresh_coeff, real=True) # 9
        >>> np.linalg.norm(im - recon) / np.linalg.norm(im)
        0.13945789596375507
        ```
        
        Up to the first marked line, everything is identical to the code for the fully sampled transform.
        The only difference in line 1 is the additional argument `subsampled=True` to obtain a subsampled transform.
        Then, in line 2, the transform of the image `im` is computed just as for the fully sampled case.
        
        The main difference to the fully sampled transform becomes visible in lines 3 and 4:
        In contrast to the fully sampled transform, where the coefficients are a 3-dimensional numpy *array*,
        the subsampled transform yields a *list* of 2-dimensional numpy arrays, with **varying shapes**.
        This shape is constant with respect to the shear `k` as long as the scale `j` and the cone `c` are kept fixed,
        but varies strongly with `j` and `c`. In face, for a *quadratic* image, the shape would only depend on the scale `j`.
        
        Since we have a list of numpy arrays instead of a single numpy array, all operations on the coefficients are more
        cumbersome to write down (using list comprehensions), but are identical in spirit to the case of the fully sampled
        transform, cf. lines 5-8.
        
        The actual reconstruction (in line 9) is exactly identical to the fully sampled case. It is interesting to note that
        only about 0.24% of the coefficients - and thus much less than the 2% for the fully sampled transform - are larger than
        the threshold. Nevertheless, the relative error is essentially the same.
        
        ## <a name="adaptivity"></a>Adaptively choosing α
        
        In the following, we show for each of the three optimality criteria ([AAR](#aar), [MAE](#mae) and [TDP](#tdp)) how our
        implementation can be used to determine the optimal value of α for a given set of images.
        
        ### <a name="aar"></a>Asymptotic approximation rate (AAR)
        
        The following code uses a grid search to determine the value of α which yields the best *asymptotic approximation rate*
        (as described in the technical report) for the given set of images:
        
        <a name="aar-session"></a>
        ```python
        >>> from AdaptiveAlpha import optimize_AAR
        
        >>> shearlet_args = {'real' : True, 'verbose' : False} # 1
        >>> images = ['./Review/log_dust_test.npy'] # 2
        >>> num_scales = 4 # 3
        >>> num_alphas = 3 # 4
        
        >>> optimize_AAR(images, num_scales, 1 / (num_alphas - 1), shearlet_args=shearlet_args)
        First step: Determine the maximum relevant value...
        alpha loop: 100%|████████████████████████████████████████████████████████| 3/3 [00:10<00:00,  3.48s/it]
        image loop: 100%|████████████████████████████████████████████████████████| 1/1 [00:01<00:00,  1.60s/it]
        Maximum relevant value: 0.04408709120918954
        
        Second step: Computing the approximation errors...
        alpha loop: 100%|████████████████████████████████████████████████████████| 3/3 [01:40<00:00, 30.99s/it]
        Image loop: 100%|████████████████████████████████████████████████████████| 1/1 [00:46<00:00, 46.89s/it]
        Thresh. loop: 100%|██████████████████████████████████████████████████████| 50/50 [00:45<00:00,  1.10it/s]
        
        Third step: Computing the approximation rates...
        Common breakpoints: [0, 4, 34, 50]
        Last common linear part: [34, 50)
        
        last slopes:
           alpha = 1.00 : -0.161932 + 0.012791
           alpha = 0.50 : -0.161932 + 0.000108
         * alpha = 0.00 : -0.161932 - 0.012900
        
        Optimal value: alpha = 0.00
        ```
        In addition to the output shown above, executing this code will display the following plot:
        <p align="center">
        <a name="aar-plot"></a>
        <img src="./Review/AAR_dust_example_plot.png" width="900"/>
        </p>
        
        
        We now briefly explain the above [code](#aar-session) and output:
        In the short program above, we first import the function `optimize_AAR` which will do the actual work.
        Then, we define the parameters to be passed to this function:
        
        1. In line 1, we determine the properties of the α-shearlet systems that will be used:
            * `'real': True` ensures that real-valued shearlets are used.
            * `'verbose': False` suppresses some output, e.g. the progress bar for precomputing the shearlet system.
        
           Another possible option would be `'subsampled': True` if one wants to use the subsampled transform. Note though
           that this is incompatible with the `'real': True` option.
        
        2. In line 2, we determine the set of images that is to be used. To ensure fast computations, we only take
           a single image for this example. Specifically, the used image is the logarithm of one of the 12 faces of cosmic
           dust data provided by CEA, as depicted in the following figure:
        
           <p align="center">
           <a name="log-dust-image"></a>
           <img src="./Review/log_dust_graphic.png" width="400"/>
           </p>
        
        3. The variable `num_scales` determines how many scales the α-shearlet systems should use.
        
        4. Since we are using a grid search (i.e., we are only considering finitely many values of α), the
           variable `num_alphas` is used to determine how many different values of α should be distinguished.
           These are then uniformly spread in the interval [0,1]. Again, to ensure fast computations, we only consider
           three different values, namely α=0, α=0.5 and α=1.
        
        Finally, we invoke the function `optimize_AAR` with the chosen parameters. As described in the technical report,
        this function does the following:
        
        1. It determines a range [0, c<sub>0</sub>] such that for c≥c<sub>0</sub>, all α-shearlet transforms yield the
           *same error* when all coefficients of absolute value ≤c are set to zero ("thresholded").
        
        2. It computes the reconstruction errors for the different values of α after thresholding the coefficients with
           a threshold of c<sub>0</sub>·b<sup>k</sup> for k=0,...,K-1.
           
           The default value is K=50. This can be adjusted by passing e.g. the argument `num_x_values=40` as an
           additional argument to `optimize_AAR`. Likewise, the base `b` (with default value `b=1.25`) can be adjusted
           by passing e.g. `base=1.3` as a further argument.
        
        3. It determines a partition of {0, ..., K-1} into at most 4 intervals such that on each of these intervals,
           each of the (logarithmic) error curves is almost linear.
           
           In the example run above, the end points of the
           resulting partition are given by `Common breakpoints: [0, 4, 34, 50]`. In this case, the resulting
           partition has only three intervals instead of 4, since on each of these intervals, the best linear approximation
           is already sufficiently good.
        
        4. For each value of α, the function then determines the *slopes* of the (logarithmic) error curve on the
           last of these intervals and compares these slopes.
        
           The optimal value of α (in this case α=0) is the one with the smallest slope, i.e., with the highest
           decay of the error. To allow for a visual comparison, `optimize_AAR` also displays a [plot](#aar-plot) of the
           (logarithmic) error curves, including the partition into the almost linear parts.
        
        ### <a name="mae"></a>Mean approximation error (MAE)
        
        The following code uses a grid search to determine the value of α which yields the best *mean approximation error*
        (as described in the technical report) for the given set of images:
        
        <a name="mae-session"></a>
        ```python
        >>> from AdaptiveAlpha import optimize_MAE
        >>> shearlet_args = {'real' : True, 'verbose' : False} # 1
        >>> images = ['./Review/log_dust_test.npy'] # 2
        >>> num_scales = 4
        >>> num_alphas = 3
        >>> optimize_MAE(images, num_scales, 1 / (num_alphas - 1), shearlet_args=shearlet_args)
        First step: Determine the maximum relevant value...
        alpha loop: 100%|████████████████████████████████████████████████████████| 3/3 [00:11<00:00,  3.61s/it]
        image loop: 100%|████████████████████████████████████████████████████████| 1/1 [00:01<00:00,  1.61s/it]
        Maximum relevant value: 0.006053772894213893
        
        Second step: Computing the approximation errors...
        alpha loop: 100%|████████████████████████████████████████████████████████| 3/3 [01:40<00:00, 31.13s/it]
        image loop: 100%|████████████████████████████████████████████████████████| 1/1 [00:48<00:00, 48.25s/it]
        thresholding loop: 100%|█████████████████████████████████████████████████| 50/50 [00:46<00:00,  1.08it/s]
        
        Final step: Computing optimal value of alpha...
        mean errors:
           alpha = 1.00 : +0.024105 + 0.000258
           alpha = 0.50 : +0.024105 + 0.000015
         * alpha = 0.00 : +0.024105 - 0.000273
        
        Optimal value: alpha = 0.00
        ```
        
        In addition to the output shown above, executing this code will display the following plot:
        <p align="center">
        <a name="mae-plot"></a>
        <img src="./Review/MAE_example_plot.png" width="800"/>
        </p>
        
        We now briefly explain the above [code](#mae-session) and output:
        In the short program above, we first import the function `optimize_MAE` which will do the actual work.
        Then, we define the parameters to be passed to this function:
        
        1. In line 1, we determine the properties of the α-shearlet systems that will be used:
            * `'real': True` ensures that real-valued shearlets are used.
            * `'verbose': False` suppresses some output, e.g. the progress bar for precomputing the shearlet system.
        
           Another possible option would be `'subsampled': True` if one wants to use the subsampled transform. Note though
           that this is incompatible with the `'real': True` option.
        
        2. In line 2, we determine the set of images that is to be used. To ensure fast computations, we only take
           a single image for this example. Specifically, the used image is the logarithm of one of the 12 faces of cosmic
           dust data provided by CEA, as depicted in the [figure](#log-dust-image) above.
        
        3. The variable `num_scales` determines how many scales the α-shearlet systems should use.
        
        4. Since we are using a grid search (i.e., we are only considering finitely many values of α), the
           variable `num_alphas` is used to determine how many different values of α should be distinguished.
           These are then uniformly spread in the interval [0,1]. Again, to ensure fast computations, we only consider
           three different values, namely α=0, α=0.5 and α=1.
        
        Finally, we invoke the function `optimize_MAE` with the chosen parameters. As described in the technical report,
        this function does the following:
        
        1. It determines a range [0, c<sub>0</sub>] such that for c≥c<sub>0</sub>, all α-shearlet transforms yield the
           *same error* when all coefficients of absolute value ≤c are set to zero ("thresholded").
        
        2. It computes the reconstruction errors for the different values of α after thresholding the coefficients with
           a threshold of c = c<sub>0</sub>·i / K for i = 0, ..., K-1.
        
           The default value is K=50. This can be adjusted by passing e.g. the argument `num_x_values=40` as an additional
           argument to `optimize_MAE`.
        
        3. For each value of α, the function then determines the *mean* of all these approximation errors. The optimal
           value of α (in this case α=0) is the one with the smallest mean approximation error.
        
           To allow for a visual comparison, `optimize_MAE` also displays a [plot](#mae-plot) of the error curves.
        
        
        
        ### <a name="tdp"></a>Thresholding Denoising Performance (TDP)
        
        In the following, we show how the denoising performance of an α-shearlet system can be used as an optimality
        criterion for adaptively choosing the correct value of α.
        
        Since the [(logarithmic) dust data](#log-dust-image) used for the previous experiments does not allow for an
        easy visual comparison between the original image and the different denoised versions, we decided to instead
        use the following 684 x 684 cartoon image (taken from [SMBC](http://smbc-comics.com)) as a toy example:
        <p align="center">
        <a name=cartoon-example></a>
        <img src="./Review/cartoon_example.png" width="400" height="400" />
        </p>
        
        
        The following code uses a grid search over different values of α to determine the value with the
        optimal denoising performance:
        <a name="tdp-session"></a>
        ```python
        >>> from AdaptiveAlpha import optimize_denoising
        
        >>> image_paths=['./Review/cartoon_example.png']       # 1
        >>> num_alphas = 3                                     # 2
        >>> num_scales = 5                                     # 3
        >>> num_noise_levels = 5                               # 4
        >>> shearlet_args = {'real' : True, 'verbose' : False} # 5
        
        >>> optimize_denoising(image_paths,
                               num_scales,
                               1 / (num_alphas - 1),
                               num_noise_levels,
                               shearlet_args=shearlet_args)
        image loop: 100%|██████████████████████████████████████████████████████████| 1/1 [01:33<00:00, 93.16s/it]
        alpha loop: 100%|██████████████████████████████████████████████████████████| 3/3 [01:33<00:00, 28.44s/it]
        noise loop: 100%|██████████████████████████████████████████████████████████| 5/5 [00:46<00:00,  9.38s/it]
        
        Averaged error over all images and all noise levels:
        alpha = 1.00: 0.0961
        alpha = 0.50: 0.0900
        alpha = 0.00: 0.0948
        
        Optimal value on whole set: alpha = 0.50
        ```
        
        In addition to the output shown above, executing the sample code also displays the following plot:
        <p align="center">
        <a name="tdp-plot"></a>
        <img src="./Review/TDP_example_plot.png" width="900"/>
        </p>
        
        We now briefly explain the above [code](#tdp-session) and output:
        First, we import from `AdaptiveAlpha.py` the function `optimize_denoising` which will do the actual work.
        We then set the parameters to this function. In the present case, we want to
        * analyze the [cartoon image](#cartoon-example) shown above,
        * use α-shearlet transforms with *real-valued* shearlets (cf. line 5, the second part of that line suppresses some output),
        * use α-shearlet transforms with 5 scales (line 4),
        * use 5 different noise levels λ (line 4) which are uniformly spaced in [0.02, 0.4],
        * compare three different values of α, which are uniformly spaced in [0,1], i.e, α=1, α=0.5 and α=0.
        
        In more realistic experiment, one would of course use a larger set of test images and consider more different values of
        α, possibly also with a larger number of different noise levels. But here, we are mainly interested in a quick
        execution time, so we keep everything small.
        
        Finally, we invoke the function `optimize_denoising` which - briefly summarized - does the following:
        1. It normalizes each of the K x K input images to have L² norm equal to 1.
           Then, for each image, each value of α  and each noise level λ in [0.02,0.4], a distorted image is calculated
           by adding artificial Gaussian noise with standard deviation σ=λ/K to the image.
        
           This standard deviation is chosen in such a way that the expected squared L² norm of the noise is λ².
        
           One can specify other ranges for the noise level than the default range [0.02, 0.4] by using the paramters
           `noise_min` and `noise_max`.
        
        2. The α-shearlet transform of the distorted image is determined. 
        
        3. Hard thresholding is performed on the set of α-shearlet coefficients. The thresholding parameter (i.e., the cutoff value)
           c is chosen scale- and noise dependent via c = mσ, with m being a scale-dependent *multiplier*.
        
           Numerical experiments showed that good results are obtained by taking m=3 for all scales except the highest and m=4
           for the highest scale. This is the default choice made in `optimize_denoising`. If desired, this default choice can be
           modified using the parameter `thresh_multiplier`.
        
        4. The inverse α-shearlet transform of the thresholded coefficients is determined and the L²-error between this
           reconstruction and the original image is calculated. 
        
        5. The optimal value of α is the one for which the L²-error averaged over all images and all noise levels is the smallest. 
        
        The [plot](#tdp-plot) which is displayed by `optimize_denoising` depicts the mean error over all images as a function of
        the employed noise level λ for different values of α. In the present case, we are only considering one image (N=1).
        Clearly, shearlets (i.e., α=0.5) turn out to be optimal for our toy example.
        
        For an eye inspection, `optimize_denoising` also saves the noisy image and the reconstructions for the largest noise level
        (λ=0.4) and the different α values to the current working directory.
        A small zoomed part of these reconstructions - together with the same part of the original image and the noisy image - can
        be seen below: 
        
        <p align="center">
        <table>
            <tr>
                <td>
                    <p align="center">
                        <img src="./Review/OriginalImage.png" width="250" height="250" />
                    </p>
                    <p align="center">original image</p>
                </td>
                <td>
                    <p align="center">
                        <img src="./Review/NoisyImage.png" width="250" height="250" />
                    </p>
                    <p align="center">
                        noisy image
                    </p>
                </td>
            </tr>
        </table>
        </p>
        
        <p align="center">
        <table> 
        <tr>
            <td>
                <p align="center">
                    <img src="./Review/ReconstructionRidgelet.png" width="250" height="250" />
                </p>
                <p align="center">
                    α=0
                </p>
            </td>
            <td>
                <p align="center">
                    <img src="./Review/ReconstructionShearlet.png" width="250" height="250" />
                </p>
                <p align="center">
                    α=0.5
                </p>
            </td>
            <td>
                <p align="center">
                    <img src="./Review/ReconstructionWavelet.png" width="250" height="250" />
                </p>
                <p align="center">
                    α=1
                </p>
            </td>
        </tr> 
        </table>
        </p>
        
        
        ## <a name="sphere-valued"></a>The α-shearlet transform for functions defined on the sphere
        
        As described in the technical report, we use a *chart-based* approach for computing the
        α-shearlet transform of functions defined on the sphere. Precisely, we use the charts provided
        by the [HEALPix](http://healpix.jpl.nasa.gov/) [pixelization](https://en.wikipedia.org/wiki/HEALPix) of the sphere,
        which divides the sphere into 12 faces and provides cartesian charts for each of these faces. The crucial property
        of this pixelization is that each pixel has exactly the same spherical area.
        
        Now, given a function `f` defined on the sphere (in a discretized version, as a so-called **HEALPix map**), one can use
        the function `get_all_faces()` from `SphereTransform.py` to obtain a family of 12 cartesian images, each of which
        represents the function `f` restricted to one of the 12 faces of the sphere. As detailed in the technical report, analyzing
        each of these 12 cartesian images and concatenating the coefficients is equivalent to analyzing `f` using the sphere-based
        α-shearlets.
        
        Conversely, one needs a way to reconstruct (a slightly modified version of) `f`, given the (possibly thresholded or
        otherwise modified) α-shearlet coefficients. To this end, one first reconstructs each of the 12 cartesian images using
        the usual α-shearlet transform and then concatenates these to obtain a function `g` defined on the sphere, via the
        function `put_all_faces` defined in `SphereTransform.py`.
        
        Here, we just briefly indicate how `put_all_faces` can be used to obtain a plot of certain (randomly selected)
        α-shearlets on the sphere. Below, we use an alpha-shearlet system with 6 scales, but only use alpha-shearlets
        from the first three scales for plotting, since the alpha-shearlets on higher scales are very small/short
        and thus unsuitable for producing a nice plot.
        ```python
        >>> import healpy as hp
        >>> import numpy as np
        >>> import matplotlib.pyplot as plt
        >>> from SphereTransform import put_all_faces as Cartesian2Sphere
        >>> from AlphaTransform import AlphaShearletTransform as AST
        
        >>> width = height = 512
        >>> alpha = 0.5
        >>> num_scales = 6 # we use six scales, so that the shearlets on lower scales are big (good for plotting)
        
        >>> trafo = AST(width, height, [alpha] * num_scales, real=True)
        >>> all_shearlets = trafo.shearlets # get a list of all shearlets
        
        >>> cartesian_faces = np.empty((12, height, width))
        >>> # for each of the 12 faces, select a random shearlet from one of the first three scales
        >>> upper_bound = trafo.scale_slice(3).start
        >>> shearlet_indices = np.random.choice(np.arange(upper_bound), size=12, replace=False)
        >>> for i, index in enumerate(shearlet_indices):
                cartesian_faces[i] = all_shearlets[index]
                # normalize, so that the different shearlets are comparable in size
                max_val = np.max(np.abs(cartesian_faces[i]))
                cartesian_faces[i] /= max_val
        
        >>> # use HEALPix charts to push the 12 cartesian faces onto the sphere
        >>> sphere_shearlets = Cartesian2Sphere(cartesian_faces)
        >>> hp.mollview(sphere_shearlets, cbar=False, hold=True)
        >>> plt.title(r"Random $\alpha$-shearlets on the sphere", fontsize=20)
        >>> plt.show()
        ```
        
        The above code produces a plot similar to the following:
        <p align="center">
        <a name="sphere-plot"></a>
        <img src="./Review/sphere_alpha_shearlets.png" width="900"/>
        </p>
        
        
        
        ## Required packages
        
        In addition to Python 3, the software requires the following Python packages:
        
        * [numpy](http://www.numpy.org/)
        * [matplotlib](http://matplotlib.org/)
        * [numexpr](https://github.com/pydata/numexpr)
        * [pyfftw](https://github.com/pyFFTW/pyFFTW)
        * [tqdm](https://github.com/tqdm/tqdm)
        * [healpy](https://github.com/healpy/healpy) (only required for [SphereTransform.py](./SphereTransform.py))
        * [PIL](http://www.pythonware.com/products/pil/), the Python Imaging Library
        * [scipy.ndimage](https://www.scipy.org/)
        
Platform: UNKNOWN
Description-Content-Type: text/markdown
