Commit 92aa7ea5 authored by Philipp Arras's avatar Philipp Arras


parent 2ca43cc8
......@@ -73,7 +73,7 @@ from .minimization.kl_energy import KL_Energy
from .sugar import *
from .plot import Plot
from .library.smooth_linear_amplitude import SLAmplitude
from .library.smooth_linear_amplitude import (SLAmplitude, CepstrumOperator)
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator,
......@@ -32,9 +32,9 @@ class LogRGSpace(StructuredDomain):
Number of grid points or numbers of gridpoints along each axis.
bindistances : float or tuple of float
Distance between two grid points along each axis. These are
measured on logarithmic scale and are constant therfore.
measured on logarithmic scale and are constant therefore.
t_0 : float or tuple of float
Coordinate of pixel ndim*(1,).
harmonic : bool, optional
Whether the space represents a grid in position or harmonic space.
Default: False.
......@@ -17,8 +17,8 @@
import numpy as np
from ..domain_tuple import DomainTuple
from import PowerSpace
from import UnstructuredDomain
from ..field import Field
from ..operators.exp_transform import ExpTransform
from ..operators.offset_operator import OffsetOperator
......@@ -32,39 +32,70 @@ def _ceps_kernel(dof_space, k, a, k0):
return a**2/(1 + (k/k0)**2)**2
def _create_cepstrum_amplitude_field(domain, cepstrum):
dim = len(domain.shape)
shape = domain.shape
q_array = domain.get_k_array()
def CepstrumOperator(target, a, k0):
'''Turns a white Gaussian random field into a smooth field on a LogRGSpace.
Composed out of three operators:
sym @ qht @ diag(sqrt_ceps),
where sym is a :class:`SymmetrizingOperator`, qht is a :class:`QHTOperator`
and ceps is the so-called cepstrum:
.. math::
\\mathrm{sqrt\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
These operators are combined in this fashion in order to generate:
- A field which is smooth, i.e. second derivatives are punished (note
that the sqrt-cepstrum is essentially proportional to 1/k**2).
- A field which is symmetric around the pixel in the middle of the space.
This is result of the :class:`SymmetrizingOperator` and needed in order to
decouple the degrees of freedom at the beginning and the end of the
amplitude whenever :class:`CepstrumOperator` is used as in
FIXME The prior on the zero mode is ...
target : LogRGSpace
Target domain of the operator, needs to be non-harmonic and
a : float
Strength of smoothness prior (positive only). FIXME
k0 : float
Cutoff of smothness prior in quefrency space (positive only). FIXME
a, k0 = float(a), float(k0)
target = DomainTuple.make(target)
if a <= 0 or k0 <= 0:
raise ValueError
if len(target) > 1 or target[0].harmonic or len(target[0].shape) > 1:
raise TypeError
qht = QHTOperator(target)
dom = qht.domain[0]
sym = SymmetrizingOperator(target)
# Compute cepstrum field
dim = len(dom.shape)
shape = dom.shape
q_array = dom.get_k_array()
# Fill all non-zero modes
no_zero_modes = (slice(1, None),)*dim
ks = q_array[(slice(None),) + no_zero_modes]
cepstrum_field = np.zeros(shape)
cepstrum_field[no_zero_modes] = cepstrum(ks)
cepstrum_field[no_zero_modes] = _ceps_kernel(dom, ks, a, k0)
# Fill zero-mode subspaces
for i in range(dim):
fst_dims = (slice(None),)*i
sl = fst_dims + (slice(1, None),)
sl2 = fst_dims + (0,)
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
return Field.from_global_data(domain, cepstrum_field)
def CepstrumOperator(domain, a, k0):
.. math::
C(k) = \\left(\\frac{a}{1+(k/k0)^2}\\right)^2
if a <= 0 or k0 <= 0:
raise ValueError
cepstrum = Field.from_global_data(dom, cepstrum_field)
qht = QHTOperator(target=domain)
dof_space = qht.domain[0]
sym = SymmetrizingOperator(domain)
kern = lambda k: _ceps_kernel(dof_space, k, a, k0)
cepstrum = _create_cepstrum_amplitude_field(dof_space, kern)
return sym @ qht @ makeOp(cepstrum.sqrt())
......@@ -26,9 +26,9 @@ from .linear_operator import LinearOperator
class QHTOperator(LinearOperator):
"""Does a Hartley transform on LogRGSpace
This operator takes a field on a LogRGSpace and transforms it
according to the Hartley transform. The zero modes are not transformed
because they are infinitely far away.
This operator takes a field on a LogRGSpace and performs an Hartley
transform. The zero modes are not transformed because they are infinitely
far away in LogRGSpaces.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment