P_Sparsity.py 5.39 KB
Newer Older
1
2
3
from numpy import zeros_like, unravel_index, sum, max
import numpy as np
from proxtoolbox.proxoperators.proxoperator import ProxOperator
4
from proxtoolbox.utils.OrbitalTomog import tile_array, bin_array
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


__all__ = ['P_Sparsity', 'P_Sparsity_real', 'P_Sparsity_Symmetric', 'P_Sparsity_Symmetric_real',
           'P_Sparsity_Superpixel', 'P_Sparsity_Superpixel_real']


class P_Sparsity(ProxOperator):
    """
    Projection subroutine for projecting onto a sparsity constraint
    """

    def __init__(self, experiment):
        """
        Initialization

        Parameters
        ----------
        config : dict - Dictionary containing the problem configuration. It must contain the following mappings:
        'sparsity_parameter' : int
        """
        self.sparsity_parameter = experiment.sparsity_parameter

27
        if hasattr(experiment, 'use_sparsity_with_support') and experiment.use_sparsity_with_support:
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
            self.support = experiment.sparsity_support.real.astype(np.uint)
            assert (sum(self.support) > 0), 'Invalid empty sparsity_support given'
        else:
            self.support = 1

        if self.sparsity_parameter > 30 or len(experiment.u0.shape) != 2:
            def value_selection(original, indices, sparsity_parameter):
                idx_for_threshold = unravel_index(indices[-sparsity_parameter], original.shape)
                threshold_val = abs(original[idx_for_threshold].get())
                return (abs(original) >= threshold_val) * original
        else:
            def value_selection(original, indices, sparsity_parameter):
                out = zeros_like(original)
                hits = indices[-sparsity_parameter:].get()
                hit_idx = [unravel_index(hit, original.shape) for hit in hits]
                for _idx in hit_idx:
                    out[_idx[0], _idx[1]] = original[_idx[0], _idx[1]]
                return out

        self.value_selection = value_selection


50
    def eval(self, u, prox_idx=None):
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
        """
        Parameters
        ----------
        u : array_like - Input, array to be projected

        Returns
        -------
        p_Sparsity : array_like, the projection
        """
        u *= self.support  # apply support (simple 1 if no support)
        p_sparse = 0 * u
        sorting = np.argsort(abs(u), axis=None)  # gives indices of sorted array in ascending order
        indices = np.asarray([unravel_index(sorting[i], u.shape) for i in range(-1 * self.sparsity_parameter, 0)])
        p_sparse[indices[:, 0], indices[:, 1]] = u[indices[:, 0], indices[:, 1]]

        return p_sparse


class P_Sparsity_real(P_Sparsity):
    """
    Projection subroutine for projecting onto a combined sparsity and realness constraint

    order: realness, sparsity
    """

76
77
    def eval(self, u, prox_idx=None):
        return super(P_Sparsity_real, self).eval(u.real, prox_idx)
78
79
80
81
82
83
84
85
86
87
88
89
90
91


class P_Sparsity_Symmetric(P_Sparsity):
    """
    Projection subroutine for projecting onto a combined sparsity and symmetry constraint

    order: symmetry, sparsity
    """

    def __init__(self, experiment):
        super(P_Sparsity_Symmetric, self).__init__(experiment)
        self.symmetry = experiment.symmetry_type  # antisymmetric = -1, symmetric = 1
        self.symmetry_axis = experiment.symmetry_axis  # -1 for last, 0 for first, 'both', 'all' or None for full flip

92
    def eval(self, u, prox_idx=None):
93
94
95
96
97
        if self.symmetry_axis in ['both', 'all']:
            mirror = np.flip(u, axis=None)
        else:
            mirror = np.flip(u, axis=self.symmetry_axis)
        inp = (u + self.symmetry * mirror) / 2
98
        out = super(P_Sparsity_Symmetric, self).eval(inp, prox_idx)
99
100
101
102
103
104
105
106
107
108
        return out


class P_Sparsity_Symmetric_real(P_Sparsity_Symmetric):
    """
    Projection subroutine for projecting onto a combined symmetry, sparsity and realness constraint

    order: realness, symmetry, sparsity
    """

109
110
    def eval(self, u, prox_idx=None):
        out = super(P_Sparsity_Symmetric_real, self).eval(u.real, prox_idx)
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
        return out


class P_Sparsity_Superpixel(P_Sparsity):
    """
    Apply sparsity on superpixels, i.e. on the binned array
    """
    def __init__(self, experiment):
        super(P_Sparsity_Superpixel, self).__init__(experiment)
        # Set the superpixel size:
        self.superpixel_size = experiment.superpixel_size
        # Bin the support:
        if self.support is not 1:
            self.support = bin_array(self.support, self.superpixel_size)
            self.support = self.support / max(self.support)
        # Assert that the binning+upsampling conserves the array size
        test_shape = tile_array(bin_array(experiment.u0, self.superpixel_size, pad_zeros=False),
                                self.superpixel_size).shape
        assert test_shape == experiment.u0.shape, 'Given array size does not allow for binning'
        # TODO: allow for padding, then cut of the remainder after tile_array

132
    def eval(self, u, prox_idx=None):
133
        binned = bin_array(abs(u), self.superpixel_size)
134
        sparse_array = super(P_Sparsity_Superpixel, self).eval(binned, prox_idx)
135
136
137
138
139
140
141
142
        mask = tile_array(sparse_array, self.superpixel_size, normalize=True) > 0
        return np.where(mask, u, 0)


class P_Sparsity_Superpixel_real(P_Sparsity_Superpixel):
    """
    Apply real-valued sparsity on superpixels, i.e. on the binned array
    """
143
144
    def eval(self, u, prox_idx=None):
        return super(P_Sparsity_Superpixel_real, self).eval(u.real, prox_idx)
145