Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • empyre/empyre
  • weber/empyre
  • wessels/empyre
  • bryan/empyre
4 results
Show changes
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 4386 deletions
pyramid/__init__.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Package for the creation and reconstruction of magnetic distributions and resulting phase maps.
Modules
-------
magcreator
Create simple magnetic distributions.
magdata
Class for the storage of magnetization data.
projector
Class for projecting given magnetization distributions.
kernel
Class for the kernel matrix representing one magnetized pixel.
phasemapper
Create magnetic and electric phase maps from magnetization data.
phasemap
Class for the storage of phase data.
analytic
Create phase maps for magnetic distributions with analytic solutions.
dataset
Class for collecting pairs of phase maps and corresponding projectors.
forwardmodel
Class which represents a phase mapping strategy.
costfunction
Class for the evaluation of the cost of a function.
reconstruction
Reconstruct magnetic distributions from given phasemaps.
Subpackages
-----------
numcore
Provides fast numerical functions for core routines.
"""
from _version import __version__
pyramid/_version.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Version string"""
__version__ = '0.1.0'
pyramid/analytic.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Create phase maps for magnetic distributions with analytic solutions.
This module provides methods for the calculation of the magnetic phase for simple geometries for
which the analytic solutions are known. These can be used for comparison with the phase
calculated by the functions from the :mod:`~pyramid.phasemapper` module.
"""
import numpy as np
from numpy import pi
from pyramid.phasemap import PhaseMap
import logging
LOG = logging.getLogger(__name__)
PHI_0 = -2067.83 # magnetic flux in T*nm²
def phase_mag_slab(dim, a, phi, center, width, b_0=1):
'''Calculate the analytic magnetic phase for a homogeneously magnetized slab.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the slab in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the slab in pixel coordinates `(z, y, x)`.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
LOG.debug('Calling phase_mag_slab')
# Function for the phase:
def phi_mag(x, y):
def F_0(x, y):
A = np.log(x**2 + y**2 + 1E-30)
B = np.arctan(x / (y+1E-30))
return x*A - 2*x + 2*y*B
return coeff * Lz * (- np.cos(phi) * (F_0(x-x0-Lx/2, y-y0-Ly/2)
- F_0(x-x0+Lx/2, y-y0-Ly/2)
- F_0(x-x0-Lx/2, y-y0+Ly/2)
+ F_0(x-x0+Lx/2, y-y0+Ly/2))
+ np.sin(phi) * (F_0(y-y0-Ly/2, x-x0-Lx/2)
- F_0(y-y0+Ly/2, x-x0-Lx/2)
- F_0(y-y0-Ly/2, x-x0+Lx/2)
+ F_0(y-y0+Ly/2, x-x0+Lx/2)))
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 define the center of a pixel,
x0 = a * (center[2] + 0.5) # hence: (cellindex + 0.5) * grid spacing
Lz, Ly, Lx = a * width[0], a * width[1], a * width[2]
coeff = b_0 / (4*PHI_0)
# Create grid:
x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, phi_mag(xx, yy))
def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
'''Calculate the analytic magnetic phase for a homogeneously magnetized disc.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
LOG.debug('Calling phase_mag_disc')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
r[center[1], center[2]] = 1E-30
result = coeff * Lz * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r <= R, 1, (R / r) ** 2)
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
Lz = a * height
R = a * radius
coeff = - pi * b_0 / (2*PHI_0)
# Create grid:
x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, phi_mag(xx, yy))
def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
'''Calculate the analytic magnetic phase for a homogeneously magnetized sphere.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the sphere in pixel coordinates `(z, y, x)`.
radius : float
The radius of the sphere in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
LOG.debug('Calling phase_mag_sphere')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
r[center[1], center[2]] = 1E-30
result = coeff * R ** 3 / r ** 2 * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r > R, 1, (1 - (1 - (r / R) ** 2) ** (3. / 2.)))
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
R = a * radius
coeff = - 2./3. * pi * b_0 / PHI_0
# Create grid:
x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, phi_mag(xx, yy))
def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
'''Calculate the analytic magnetic phase for a vortex state disc.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
a : float
The grid spacing in nm.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`, which is also the vortex center.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
LOG.debug('Calling phase_mag_vortex')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
result = coeff * np.where(r <= R, r - R, 0)
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
Lz = a * height
R = a * radius
coeff = pi * b_0 * Lz / PHI_0
# Create grid:
x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return PhaseMap(a, phi_mag(xx, yy))
pyramid/converter.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 19 08:48:45 2014
@author: Jan
""" # TODO: Docstring
# TODO: put into other class
# TODO: use 3 components (more complex)
# TODO: take masks into account
import numpy as np
class IndexConverter3components(object):
def __init__(self, dim):
self.dim = dim
self.size_3d = np.prod(dim)
self.size_2d = dim[2]*dim[1]
def ind_to_coord(self, ind):
m, remain = int(ind/self.size_3d), ind % self.size_3d
z, remain = int(remain/self.size_2d), remain%self.size_2d
y, x = int(remain/self.dim[2]), remain%self.dim[2]
coord = m, z, y, x
return coord
def coord_to_ind(self, coord):
z, y, x = coord
ind = [i*self.size_3d + z*self.size_2d + y*self.dim[2] + x for i in range(3)]
return ind
def find_neighbour_ind(self, coord):
z, y, x = coord
t, d = (z-1, y, x), (z+1, y, x) # t: top, d: down
f, b = (z, y-1, x), (z, y+1, x) # f: front, b: back
l, r = (z, y, x-1), (z, y, x+1) # l: left, r: right
neighbours = [self.coord_to_ind(i) for i in [t, d, f, b, l, r]]
return np.reshape(np.swapaxes(neighbours, 0, 1), (3, 3, 2))
class IndexConverter(object):
def __init__(self, dim):
self.dim = dim
self.size_2d = dim[2]*dim[1]
def ind_to_coord(self, ind):
z, remain = int(ind/self.size_2d), ind%self.size_2d
y, x = int(remain/self.dim[2]), remain%self.dim[2]
coord = z, y, x
return coord
def coord_to_ind(self, coord):
z, y, x = coord
ind = z*self.size_2d + y*self.dim[2] + x
return ind
def find_neighbour_ind(self, coord):
z, y, x = coord
t, d = (z-1, y, x), (z+1, y, x) # t: top, d: down
f, b = (z, y-1, x), (z, y+1, x) # f: front, b: back
l, r = (z, y, x-1), (z, y, x+1) # l: left, r: right
neighbours = [self.coord_to_ind(i) for i in [t, d, f, b, l, r]]
return neighbours
return np.reshape(neighbours, (3, 2))
pyramid/costfunction.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""This module provides the :class:`~.Costfunction` class which represents a strategy to calculate
the so called `cost` of a threedimensional magnetization distribution."""
import numpy as np
from scipy.sparse.linalg import LinearOperator
from scipy.sparse import eye
from pyramid.forwardmodel import ForwardModel
from pyramid.regularisator import ZeroOrderRegularisator
import logging
class Costfunction(object):
'''Class for calculating the cost of a 3D magnetic distributions in relation to 2D phase maps.
Represents a strategy for the calculation of the `cost` of a 3D magnetic distribution in
relation to two-dimensional phase maps. The `cost` is a measure for the difference of the
simulated phase maps from the magnetic distributions to the given set of phase maps and from
a priori knowledge represented by a :class:`~.Regularisator` object. Furthermore this class
provides convenient methods for the calculation of the derivative :func:`~.jac` or the product
with the Hessian matrix :func:`~.hess_dot` of the costfunction, which can be used by
optimizers. All required data should be given in a :class:`~DataSet` object.
Attributes
----------
data_set: :class:`~dataset.DataSet`
:class:`~dataset.DataSet` object, which stores all information for the cost calculation.
regularisator : :class:`~.Regularisator`
Regularisator class that's responsible for the regularisation term.
regularisator: :class:`~Regularisator`
y : :class:`~numpy.ndarray` (N=1)
Vector which lists all pixel values of all phase maps one after another.
fwd_model : :class:`~.ForwardModel`
The Forward model instance which should be used for the simulation of the phase maps which
will be compared to `y`.
Se_inv : :class:`~numpy.ndarray` (N=2), optional
Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
being the length of the targetvector y (vectorized phase map information).
'''
LOG = logging.getLogger(__name__+'.Costfunction')
def __init__(self, data_set, regularisator):
self.LOG.debug('Calling __init__')
self.data_set = data_set
self.fwd_model = ForwardModel(data_set)
self.regularisator = regularisator
# Extract important information:
self.y = data_set.phase_vec
self.Se_inv = data_set.Se_inv
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(data_set=%r, fwd_model=%r, regularisator=%r)' % \
(self.__class__, self.data_set, self.fwd_model, self.regularisator)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'Costfunction(data_set=%s, fwd_model=%s, regularisator=%s)' % \
(self.data_set, self.fwd_model, self.regularisator)
def __call__(self, x):
self.LOG.debug('Calling __call__')
delta_y = self.fwd_model(x)-self.y
return delta_y.dot(self.Se_inv.dot(delta_y)) + self.regularisator(x)
def jac(self, x):
'''Calculate the derivative of the costfunction for a given magnetization distribution.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Vectorized magnetization distribution, for which the Jacobi vector is calculated.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Jacobi vector which represents the cost derivative of all voxels of the magnetization.
'''
self.LOG.debug('Calling jac')
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x)-self.y))
+ self.regularisator.jac(x))
def hess_dot(self, x, vector):
'''Calculate the product of a `vector` with the Hessian matrix of the costfunction.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Vectorized magnetization distribution at which the Hessian is calculated. The Hessian
is constant in this case, thus `x` can be set to None (it is not used int the
computation). It is implemented for the case that in the future nonlinear problems
have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized magnetization distribution which is multiplied by the Hessian.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the input `vector` with the Hessian matrix of the costfunction.
'''
self.LOG.debug('Calling hess_dot')
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
+ self.regularisator.hess_dot(x, vector))
class CFAdapterScipyCG(LinearOperator):
'''Adapter class making the :class:`~.Costfunction` class accessible for scipy cg methods.
This class provides an adapter for the :class:`~.Costfunction` to be usable with the
:func:`~.scipy.sparse.linalg.cg` function. the :func:`~.matvec` function is overwritten to
implement a multiplication with the Hessian of the adapted costfunction. This is used in the
:func:`~pyramid.reconstruction.optimise_sparse_cg` function of the
:mod:`~pyramid.reconstruction` module.
Attributes
----------
cost : :class:`~.Costfunction`
Costfunction which should be made usable in the :func:`~.scipy.sparse.linalg.cg` function.
'''
# TODO: make obsolete!
LOG = logging.getLogger(__name__+'.CFAdapterScipyCG')
def __init__(self, cost):
self.LOG.debug('Calling __init__')
self.cost = cost
def matvec(self, vector):
'''Matrix-vector multiplication with the Hessian of the adapted costfunction.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vector which will be multiplied by the Hessian matrix provided by the adapted
costfunction.
'''
self.LOG.debug('Calling matvec')
return self.cost.hess_dot(None, vector)
@property
def shape(self):
return (self.cost.data_set.m, self.cost.data_set.m)
@property
def dtype(self):
return np.dtype("d")
pyramid/dataset.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""This module provides the :class:`~.DataSet` class for the collection of phase maps
and additional data like corresponding projectors."""
import numpy as np
from numbers import Number
import scipy.sparse as sp
import matplotlib.pyplot as plt
from pyramid.phasemap import PhaseMap
from pyramid.phasemapper import PhaseMapperRDFC
from pyramid.projector import Projector
from pyramid.kernel import Kernel
import logging
class DataSet(object):
'''Class for collecting phase maps and corresponding projectors.
Represents a collection of (e.g. experimentally derived) phase maps, stored as
:class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
objects. At creation, the grid spacing `a` and the dimension `dim` of the magnetization
distribution have to be given. Data can be added via the :func:`~.append` method, where
a :class:`~.PhaseMap`, a :class:`~.Projector` and additional info have to be given.
Attributes
----------
a: float
The grid spacing in nm.
dim: tuple (N=3)
Dimensions of the 3D magnetization distribution.
phase_maps:
A list of all stored :class:`~.PhaseMap` objects.
b_0: double
The saturation induction in `T`.
mask: :class:`~numpy.ndarray` (N=3), optional
A boolean mask which defines the magnetized volume in 3D.
projectors: list of :class:`~.Projector`
A list of all stored :class:`~.Projector` objects.
phase_maps: list of :class:`~.PhaseMap`
A list of all stored :class:`~.PhaseMap` objects.
phase_vec: :class:`~numpy.ndarray` (N=1)
The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
n: int
Size of the image space.
m: int
Size of the input space.
'''
LOG = logging.getLogger(__name__+'.DataSet')
@property
def n(self):
return np.sum([len(p.phase_vec) for p in self.phase_maps])
@property
def Se_inv(self):
# TODO: better implementation, maybe get-method? more flexible? input in append?
return sp.eye(self.n)
@property
def phase_vec(self):
return np.concatenate([p.phase_vec for p in self.phase_maps])
@property
def hook_points(self):
result = [0]
for i, phase_map in enumerate(self.phase_maps):
result.append(result[i]+np.prod(phase_map.dim_uv))
return result
@property
def phase_mappers(self):
dim_uv_list = np.unique([p.dim_uv for p in self.projectors])
kernel_list = [Kernel(self.a, tuple(dim_uv)) for dim_uv in dim_uv_list]
return {kernel.dim_uv: PhaseMapperRDFC(kernel) for kernel in kernel_list}
def __init__(self, a, dim, b_0=1, mask=None):
self.LOG.debug('Calling __init__')
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
assert isinstance(dim, tuple) and len(dim) == 3, \
'Dimension has to be a tuple of length 3!'
if mask is not None:
assert mask.shape == dim, 'Mask dimensions must match!'
self.m = 3 * np.sum(self.mask)
else:
self.m = 3 * np.prod(dim)
self.a = a
self.dim = dim
self.b_0 = b_0
self.mask = mask
self.phase_maps = []
self.projectors = []
self.LOG.debug('Created: '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(a=%r, dim=%r, b_0=%r)' % (self.__class__, self.a, self.dim, self.b_0)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'DataSet(a=%s, dim=%s, b_0=%s)' % (self.a, self.dim, self.b_0)
def append(self, phase_map, projector): # TODO: include Se_inv or 2D mask??
'''Appends a data pair of phase map and projection infos to the data collection.`
Parameters
----------
phase_map: :class:`~.PhaseMap`
A :class:`~.PhaseMap` object which should be added to the data collection.
projector: :class:`~.Projector`
A :class:`~.Projector` object which should be added to the data collection.
Returns
-------
None
'''
self.LOG.debug('Calling append')
assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
'Argument has to be a tuple of a PhaseMap and a Projector object!'
assert projector.dim == self.dim, '3D dimensions must match!'
assert phase_map.dim_uv == projector.dim_uv, 'Projection dimensions (dim_uv) must match!'
self.phase_maps.append(phase_map)
self.projectors.append(projector)
def create_phase_maps(self, mag_data):
'''Create a list of phasemaps with the projectors in the dataset for a given
:class:`~.MagData` object.
Parameters
----------
mag_data : :class:`~.MagData`
Magnetic distribution to which the projectors of the dataset should be applied.
Returns
-------
phase_maps : list of :class:`~.phasemap.PhaseMap`
A list of the phase maps resulting from the projections specified in the dataset.
'''
return [self.phase_mappers[projector.dim_uv](projector(mag_data))
for projector in self.projectors]
def display_phase(self, mag_data=None, title='Phase Map',
cmap='RdBu', limit=None, norm=None):
'''Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
Parameters
----------
mag_data : :class:`~.MagData`, optional
Magnetic distribution to which the projectors of the dataset should be applied. If not
given, the phase_maps in the dataset are used.
title : string, optional
The main part of the title of the plots. The default is 'Phase Map'. Additional
projector info is appended to this.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plots as a string.
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
Returns
-------
None
'''
self.LOG.debug('Calling display')
if mag_data is not None:
phase_maps = self.create_phase_maps(mag_data)
else:
phase_maps = self.phase_maps
[phase_map.display_phase('{} ({})'.format(title, self.projectors[i].get_info()),
cmap, limit, norm)
for (i, phase_map) in enumerate(phase_maps)]
plt.show()
def display_combined(self, mag_data=None, title='Combined Plot', cmap='RdBu', limit=None,
norm=None, gain=1, interpolation='none', grad_encode='bright'):
'''Display all phasemaps and the resulting color coded holography images.
Parameters
----------
mag_data : :class:`~.MagData`, optional
Magnetic distribution to which the projectors of the dataset should be applied. If not
given, the phase_maps in the dataset are used.
title : string, optional
The title of the plot. The default is 'Combined Plot'.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
gain : float, optional
The gain factor for determining the number of contour lines in the holographic
contour map. The default is 1.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
grad_encode: {'bright', 'dark', 'color', 'none'}, optional
Encoding mode of the phase gradient. 'none' produces a black-white image, 'color' just
encodes the direction (without gradient strength), 'dark' modulates the gradient
strength with a factor between 0 and 1 and 'bright' (which is the default) encodes
the gradient strength with color saturation.
Returns
-------
None
'''
self.LOG.debug('Calling display_combined')
if mag_data is not None:
phase_maps = self.create_phase_maps(mag_data)
else:
phase_maps = self.phase_maps
[phase_map.display_combined('{} ({})'.format(title, self.projectors[i].get_info()),
cmap, limit, norm, gain, interpolation, grad_encode)
for (i, phase_map) in enumerate(phase_maps)]
plt.show()
# TODO: method for constructing 3D mask from 2D masks?
pyramid/forwardmodel.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""This module provides the :class:`~.ForwardModel` class which represents a strategy to map a
threedimensional magnetization distribution onto a two-dimensional phase map."""
import numpy as np
from pyramid.magdata import MagData
import logging
class ForwardModel(object):
'''Class for mapping 3D magnetic distributions to 2D phase maps.
Represents a strategy for the mapping of a 3D magnetic distribution to two-dimensional
phase maps. Can handle a list of `projectors` of :class:`~.Projector` objects, which describe
different projection angles, so many phase_maps can be created from one magnetic distribution.
All required data should be given in a :class:`~DataSet` object.
Attributes
----------
data_set: :class:`~dataset.DataSet`
:class:`~dataset.DataSet` object, which stores all required information calculation.
projectors : list of :class:`~.Projector`
A list of all :class:`~.Projector` objects representing the projection directions.
kernel : :class:`~.Kernel`
A kernel which describes the phasemapping of the 2D projected magnetization distribution.
a : float
The grid spacing in nm.
dim : tuple (N=3)
Dimensions of the 3D magnetic distribution.
n: int
Size of the image space. Number of pixels of the 2-dimensional projected grid.
m: int
Size of the input space. Number of voxels of the 3-dimensional grid.
'''
LOG = logging.getLogger(__name__+'.ForwardModel')
def __init__(self, data_set):
self.LOG.debug('Calling __init__')
self.data_set = data_set
self.phase_mappers = data_set.phase_mappers
self.m = data_set.m
self.n = data_set.n
self.hook_points = data_set.hook_points
self.mag_data = MagData(data_set.a, np.zeros((3,)+data_set.dim))
self.LOG.debug('Creating '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(data_set=%r)' % (self.__class__, self.data_set)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'ForwardModel(data_set=%s)' % (self.data_set)
def __call__(self, x):
self.LOG.debug('Calling __call__')
self.mag_data.set_vector(x, self.data_set.mask)
# TODO: Multiprocessing
result = np.zeros(self.n)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
phase_map = self.phase_mappers[projector.dim_uv](projector(self.mag_data))
result[hp[i]:hp[i+1]] = phase_map.phase_vec
return np.reshape(result, -1)
def jac_dot(self, x, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Evaluation point of the jacobi-matrix. The Jacobi matrix is constant for a linear
problem, thus `x` can be set to None (it is not used int the computation). It is
implemented for the case that in the future nonlinear problems have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of the 3D magnetization distribution. First the `x`, then the `y` and
lastly the `z` components are listed.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the input
`vector`.
'''
self.LOG.debug('Calling jac_dot')
self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.n)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
mag_vec = self.mag_data.mag_vec
res = self.phase_mappers[projector.dim_uv].jac_dot(projector.jac_dot(mag_vec))
result[hp[i]:hp[i+1]] = res.flatten()
return result
def jac_T_dot(self, x, vector):
''''Calculate the product of the transposed Jacobi matrix with a given `vector`.
Parameters
----------
x : :class:`~numpy.ndarray` (N=1)
Evaluation point of the jacobi-matrix. The jacobi matrix is constant for a linear
problem, thus `x` can be set to None (it is not used int the computation). Is used
for the case that in the future nonlinear problems have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of all 2D phase maps one after another in one vector.
Returns
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the transposed Jacobi matrix (which is not explicitely calculated) with
the input `vector`.
'''
self.LOG.debug('Calling jac_T_dot')
result = np.zeros(self.m)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
vec = vector[hp[i]:hp[i+1]]
result += projector.jac_T_dot(self.phase_mappers[projector.dim_uv].jac_T_dot(vec))
self.mag_data.mag_vec = result
return self.mag_data.get_vector(self.data_set.mask)
pyramid/kernel.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""This module provides the :class:`~.Kernel` class, representing the phase contribution of one
single magnetized pixel."""
import numpy as np
import logging
PHI_0 = -2067.83 # magnetic flux in T*nm²
class Kernel(object):
'''Class for calculating kernel matrices for the phase calculation.
Represents the phase of a single magnetized pixel for two orthogonal directions (`u` and `v`),
which can be accessed via the corresponding attributes. The default elementary geometry is
`disc`, but can also be specified as the phase of a `slab` representation of a single
magnetized pixel. During the construction, a few attributes are calculated that are used in
the convolution during phase calculation in the different :class:`~Phasemapper` classes.
An instance of the :class:`~.Kernel` class can be called as a function with a `vector`,
which represents the projected magnetization onto a 2-dimensional grid.
Attributes
----------
a : float
The grid spacing in nm.
dim_uv : tuple of int (N=2), optional
Dimensions of the 2-dimensional projected magnetization grid from which the phase should
be calculated.
b_0 : float, optional
Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
numcore : boolean, optional
Boolean choosing if Cython enhanced routines from the :mod:`~.pyramid.numcore` module
should be used. Default is True.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
u : :class:`~numpy.ndarray` (N=3)
The phase contribution of one pixel magnetized in u-direction.
v : :class:`~numpy.ndarray` (N=3)
The phase contribution of one pixel magnetized in v-direction.
u_fft : :class:`~numpy.ndarray` (N=3)
The real FFT of the phase contribution of one pixel magnetized in u-direction.
v_fft : :class:`~numpy.ndarray` (N=3)
The real FFT of the phase contribution of one pixel magnetized in v-direction.
dim_fft : tuple of int (N=2)
Dimensions of the grid, which is used for the FFT. Calculated by adding the dimensions
`dim_uv` of the magnetization grid and the dimensions of the kernel (given by
``2*dim_uv-1``)
and increasing to the next multiple of 2 (for faster FFT).
slice_fft : tuple (N=2) of :class:`slice`
A tuple of :class:`slice` objects to extract the original field of view from the increased
size (`size_fft`) of the grid for the FFT-convolution.
'''
LOG = logging.getLogger(__name__+'.Kernel')
def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
self.LOG.debug('Calling __init__')
# Function for the phase of an elementary geometry:
def get_elementary_phase(geometry, n, m, a):
if geometry == 'disc':
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return m / (n**2 + m**2 + 1E-30) * in_or_out
elif geometry == 'slab':
def F_a(n, m):
A = np.log(a**2 * (n**2 + m**2))
B = np.arctan(n / m)
return n*A - 2*n + 2*m*B
return 0.5 * (F_a(n-0.5, m-0.5) - F_a(n+0.5, m-0.5)
- F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
# Set basic properties:
self.dim_uv = dim_uv # Size of the FOV, not the kernel (kernel is bigger)!
self.size = np.prod(dim_uv) # Pixel count
self.a = a
self.numcore = numcore
self.geometry = geometry
# Calculate kernel (single pixel phase):
coeff = -b_0 * a**2 / (2*PHI_0)
v_dim, u_dim = dim_uv
u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
uu, vv = np.meshgrid(u, v)
self.u = coeff * get_elementary_phase(geometry, uu, vv, a)
self.v = coeff * get_elementary_phase(geometry, vv, uu, a)
# Calculate Fourier trafo of kernel components:
dim_combined = 3*np.array(dim_uv) - 2 # (dim_uv-1) + dim_uv + (dim_uv-1) mag + kernel
self.dim_fft = 2 ** np.ceil(np.log2(dim_combined)).astype(int) # next multiple of 2
self.slice_fft = (slice(dim_uv[0]-1, 2*dim_uv[0]-1), slice(dim_uv[1]-1, 2*dim_uv[1]-1))
self.u_fft = np.fft.rfftn(self.u, self.dim_fft)
self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, numcore=%r, geometry=%r)' % \
(self.__class__, self.a, self.dim_uv, self.numcore, self.geometry)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'Kernel(a=%s, dim_uv=%s, numcore=%s, geometry=%s)' % \
(self.a, self.dim_uv, self.numcore, self.geometry)
pyramid/logging.ini deleted 100644 → 0
View file @ 04d9ed45
[formatters]
keys=console,file
[formatter_console]
format=%(asctime)s: %(levelname)s @ <%(name)s>:
--> %(message)s
datefmt=%H:%M:%S
class=logging.Formatter
[formatter_file]
format=%(asctime)s: %(levelname)-8s @ <%(name)s>: %(message)s
datefmt=%Y-%m-%d %H:%M:%S
class=logging.Formatter
[handlers]
keys=console,file
[handler_console]
class=logging.StreamHandler
level=WARNING
formatter=console
args=tuple()
[handler_file]
class=logging.FileHandler
level=DEBUG
formatter=file
args=('logfile.log', 'w')
[loggers]
keys=root
[logger_root]
level=DEBUG
handlers=console,file
\ No newline at end of file
pyramid/magcreator.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Create simple magnetic distributions.
The :mod:`~.magcreator` module is responsible for the creation of simple distributions of
magnetic moments. In the :class:`~.Shapes` class, you can find several general shapes for the
3-dimensional volume that should be magnetized (e.g. slabs, spheres, discs or single pixels).
These magnetic shapes are then used as input for the creating functions (or you could specify the
volume yourself as a 3-dimensional boolean matrix or a matrix with values in the range from 0 to 1,
which modifies the magnetization amplitude). The specified volume can either be magnetized
homogeneously with the :func:`~.create_mag_dist_homog` function by specifying the magnetization
direction, or in a vortex state with the :func:`~.create_mag_dist_vortex` by specifying the vortex
center.
"""
import numpy as np
from numpy import pi
import abc
import logging
LOG = logging.getLogger(__name__)
class Shapes(object):
'''Abstract class containing functions for generating magnetic shapes.
The :class:`~.Shapes` class is a collection of some methods that return a 3-dimensional
matrix that represents the magnetized volume and consists of values between 0 and 1.
This matrix is used in the functions of the :mod:`~.magcreator` module to create
:class:`~pyramid.magdata.MagData` objects which store the magnetic informations.
'''
__metaclass__ = abc.ABCMeta
LOG = logging.getLogger(__name__+'.Shapes')
@classmethod
def slab(cls, dim, center, width):
'''Create the shape of a slab.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the slab in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the slab in pixel coordinates `(z, y, x)`.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
cls.LOG.debug('Calling slab')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a tuple of length 3!'
mag_shape = np.array([[[abs(x - center[2]) <= width[2] / 2
and abs(y - center[1]) <= width[1] / 2
and abs(z - center[0]) <= width[0] / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def disc(cls, dim, center, radius, height, axis='z'):
'''Create the shape of a cylindrical disc in x-, y-, or z-direction.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
axis : {'z', 'y', 'x'}, optional
The orientation of the disc axis. The default is 'z'.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
cls.LOG.debug('Calling disc')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
if axis == 'z':
mag_shape = np.array([[[np.hypot(x - center[2], y - center[1]) <= radius
and abs(z - center[0]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
elif axis == 'y':
mag_shape = np.array([[[np.hypot(x - center[2], z - center[0]) <= radius
and abs(y - center[1]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
elif axis == 'x':
mag_shape = np.array([[[np.hypot(y - center[1], z - center[0]) <= radius
and abs(x - center[2]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def ellipse(cls, dim, center, width, height, axis='z'):
'''Create the shape of an elliptical cylinder in x-, y-, or z-direction.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the ellipse in pixel coordinates `(z, y, x)`.
width : tuple (N=2)
Length of the two axes of the ellipse.
height : float
The height of the ellipse in pixel coordinates.
axis : {'z', 'y', 'x'}, optional
The orientation of the ellipse axis. The default is 'z'.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
cls.LOG.debug('Calling ellipse')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (2,), 'Parameter width has to be a a tuple of length 2!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
if axis == 'z':
mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
(y-center[1])/(width[0]/2.)) <= 1
and abs(z - center[0]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
elif axis == 'y':
mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
(z-center[0])/(width[0]/2.)) <= 1
and abs(y-center[1]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
elif axis == 'x':
mag_shape = np.array([[[np.hypot((y-center[1])/(width[1]/2.),
(z-center[0])/(width[0]/2.)) <= 1
and abs(z - center[0]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def sphere(cls, dim, center, radius):
'''Create the shape of a sphere.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the sphere in pixel coordinates `(z, y, x)`.
radius : float
The radius of the sphere in pixel coordinates.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
cls.LOG.debug('Calling sphere')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
mag_shape = np.array([[[np.sqrt((x-center[2])**2
+ (y-center[1])**2
+ (z-center[0])**2) <= radius
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def ellipsoid(cls, dim, center, width):
'''Create the shape of an ellipsoid.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the ellipsoid in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the ellipsoid `(z, y, x)`.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
cls.LOG.debug('Calling ellipsoid')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a a tuple of length 3!'
mag_shape = np.array([[[np.sqrt((x-center[2])**2/(width[2]/2)**2
+ (y-center[1])**2/(width[1]/2)**2
+ (z-center[0])**2/(width[0]/2)**2) <= 1
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def filament(cls, dim, pos, axis='y'):
'''Create the shape of a filament.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
pos : tuple (N=2)
The position of the filament in pixel coordinates `(coord1, coord2)`.
`coord1` and `coord2` stand for the two axes, which are perpendicular to `axis`. For
the default case (`axis = y`), it is `(coord1, coord2) = (z, x)`.
axis : {'y', 'x', 'z'}, optional
The orientation of the filament axis. The default is 'y'.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
cls.LOG.debug('Calling filament')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pos) == (2,), 'Parameter pos has to be a tuple of length 2!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
mag_shape = np.zeros(dim)
if axis == 'z':
mag_shape[:, pos[0], pos[1]] = 1
elif axis == 'y':
mag_shape[pos[0], :, pos[1]] = 1
elif axis == 'x':
mag_shape[pos[0], pos[1], :] = 1
return mag_shape
@classmethod
def pixel(cls, dim, pixel):
'''Create the shape of a single pixel.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
pixel : tuple (N=3)
The coordinates of the pixel `(z, y, x)`.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
cls.LOG.debug('Calling pixel')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pixel) == (3,), 'Parameter pixel has to be a tuple of length 3!'
mag_shape = np.zeros(dim)
mag_shape[pixel] = 1
return mag_shape
def create_mag_dist_homog(mag_shape, phi, theta=pi/2, magnitude=1):
'''Create a 3-dimensional magnetic distribution of a homogeneously magnetized object.
Parameters
----------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see Shapes.* for examples).
phi : float
The azimuthal angle, describing the direction of the magnetized object.
theta : float, optional
The polar angle, describing the direction of the magnetized object.
The default is pi/2, which means, the z-component is zero.
magnitude : float, optional
The relative magnitude for the magnetic shape. The default is 1.
Returns
-------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The magnetic distribution as a tuple of the 3 components in
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
LOG.debug('Calling create_mag_dist_homog')
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
z_mag = np.ones(dim) * np.cos(theta) * mag_shape * magnitude
y_mag = np.ones(dim) * np.sin(theta) * np.sin(phi) * mag_shape * magnitude
x_mag = np.ones(dim) * np.sin(theta) * np.cos(phi) * mag_shape * magnitude
return np.array([x_mag, y_mag, z_mag])
def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
'''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
Parameters
----------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see :class:`.~Shapes` for examples).
center : tuple (N=2 or N=3), optional
The vortex center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis
is is discarded. Is set to the center of the field of view if not specified. The vortex
center has to be between two pixels.
magnitude : float, optional
The relative magnitude for the magnetic shape. The default is 1. Negative signs reverse
the vortex direction (left-handed instead of right-handed).
axis : {'z', 'y', 'x'}, optional
The orientation of the vortex axis. The default is 'z'.
Returns
-------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The magnetic distribution as a tuple of the 3 components in
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
LOG.debug('Calling create_mag_dist_vortex')
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
assert center is None or len(center) in {2, 3}, \
'Vortex center has to be defined in 3D or 2D or not at all!'
if center is None:
center = (int(dim[1]/2)-0.5, int(dim[2]/2)-0.5) # center has to be between (!) two pixels
if axis == 'z':
if len(center) == 3: # if a 3D-center is given, just take the x and y components
center = (center[1], center[2])
u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
v = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=0)
phi = np.tile(phi, (dim[0], 1, 1))
z_mag = np.zeros(dim)
y_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
x_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
elif axis == 'y':
if len(center) == 3: # if a 3D-center is given, just take the x and z components
center = (center[0], center[2])
u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=1)
phi = np.tile(phi, (1, dim[1], 1))
z_mag = np.ones(dim) * np.sin(phi) * mag_shape * magnitude
y_mag = np.zeros(dim)
x_mag = np.ones(dim) * np.cos(phi) * mag_shape * magnitude
elif axis == 'x':
if len(center) == 3: # if a 3D-center is given, just take the z and y components
center = (center[0], center[1])
u = np.linspace(-center[1], dim[1]-1-center[1], dim[1])
v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=2)
phi = np.tile(phi, (1, 1, dim[2]))
z_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
y_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
x_mag = np.zeros(dim)
return np.array([x_mag, y_mag, z_mag])
pyramid/magdata.py deleted 100644 → 0
View file @ 04d9ed45
This diff is collapsed.
pyramid/numcore/__init__.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Package for numerical core routines which enable fast computations.
Modules
-------
phasemapper_core
Provides numerical core routines for :class:`~.PhaseMapper` class.
kernel_core
Provides numerical core routines for the :class:`~.Kernel` class.
Notes
-----
Packages are written in `pyx`-format for the use with :mod:`~.cython`.
"""
pyramid/numcore/phasemapper_core.pyx deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Numerical core routines for the :mod:`~.pyramid.phasemapper` module.
Provides helper functions to speed up phase mapping calculations in the
:mod:`~pyramid.phasemapper` module, by using C-speed for the for-loops and by omitting boundary
and wraparound checks.
"""
import numpy as np
import math
cimport cython
cimport numpy as np
@cython.boundscheck(False)
@cython.wraparound(False)
def phasemapper_real_convolve(
unsigned int v_dim, unsigned int u_dim,
double[:, :] v_phi, double[:, :] u_phi,
double[:, :] v_mag, double[:, :] u_mag,
double[:, :] phase, float threshold):
'''Numerical core routine for the phase calculation using the real space approach.
Parameters
----------
v_dim, u_dim : int
Dimensions of the projection along the two major axes.
v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
v_mag, u_mag : :class:`~numpy.ndarray` (N=2)
Magnetization components in `u`- and `v`-direction.
phase : :class:`~numpy.ndarray` (N=2)
Matrix in which the resulting magnetic phase map should be stored.
threshold : float
The `threshold` determines which pixels contribute to the magnetic phase.
Returns
-------
None
'''
cdef unsigned int i, j, p, q, p_c, q_c
cdef double u_m, v_m
for j in range(v_dim):
for i in range(u_dim):
u_m = u_mag[j, i]
v_m = v_mag[j, i]
p_c = u_dim - 1 - i
q_c = v_dim - 1 - j
if abs(u_m) > threshold:
for q in range(v_dim):
for p in range(u_dim):
phase[q, p] += u_m * u_phi[q_c + q, p_c + p]
if abs(v_m) > threshold:
for q in range(v_dim):
for p in range(u_dim):
phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
@cython.boundscheck(False)
@cython.wraparound(False)
def jac_dot_real_convolve(
unsigned int v_dim, unsigned int u_dim,
double[:, :] u_phi, double[:, :] v_phi,
double[:] vector,
double[:] result):
'''Numerical core routine for the Jacobi matrix multiplication for the phase mapper.
Parameters
----------
v_dim, u_dim : int
Dimensions of the projection along the two major axes.
v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
vector : :class:`~numpy.ndarray` (N=1)
Input vector which should be multiplied by the Jacobi-matrix.
result : :class:`~numpy.ndarray` (N=1)
Vector in which the result of the multiplication should be stored.
Returns
-------
None
Notes
-----
The strategy involves iterating over the magnetization first and over the kernel in an inner
iteration loop.
'''
cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, r, u, v, size
size = u_dim * v_dim # Number of pixels
s = 0 # Current contributing pixel (numbered consecutively)
# Iterate over all contributingh pixels:
for j in range(v_dim):
for i in range(u_dim):
u_min = (u_dim-1) - i # u_max = u_min + u_dim
v_min = (v_dim-1) - j # v_max = v_min + v_dim
r = 0 # Current result component / affected pixel (numbered consecutively)
# Go over the current kernel cutout [v_min:v_max, u_min:u_max]:
for v in range(v_min, v_min + v_dim):
for u in range(u_min, u_min + u_dim):
result[r] += vector[s] * u_phi[v, u]
result[r] -= vector[s+size] * v_phi[v, u]
r += 1
s += 1
# TODO: linearize u and v, too?
@cython.boundscheck(False)
@cython.wraparound(False)
def jac_T_dot_real_convolve(
unsigned int v_dim, unsigned int u_dim,
double[:, :] u_phi, double[:, :] v_phi,
double[:] vector,
double[:] result):
'''Core routine for the transposed Jacobi multiplication for the phase mapper.
Parameters
----------
v_dim, u_dim : int
Dimensions of the projection along the two major axes.
v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
vector : :class:`~numpy.ndarray` (N=1)
Input vector which should be multiplied by the transposed Jacobi-matrix.
result : :class:`~numpy.ndarray` (N=1)
Vector in which the result of the multiplication should be stored.
Returns
-------
None
Notes
-----
The strategy involves iterating over the magnetization first and over the kernel in an inner
iteration loop.
'''
cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, r, u, v, size
size = u_dim * v_dim # Number of pixels
r = 0 # Current result component / contributing pixel (numbered consecutively)
# Iterate over all contributingh pixels:
for j in range(v_dim):
for i in range(u_dim):
u_min = (u_dim-1) - i # u_max = u_min + u_dim
v_min = (v_dim-1) - j # v_max = v_min + v_dim
s = 0 # Current affected pixel (numbered consecutively)
# Go over the current kernel cutout [v_min:v_max, u_min:u_max]:
for v in range(v_min, v_min + v_dim):
for u in range(u_min, u_min + u_dim):
result[r] += vector[s] * u_phi[v, u]
result[r+size] -= vector[s] * v_phi[v, u]
s += 1
r += 1
# TODO: linearize u and v, too?
pyramid/phasemap.py deleted 100644 → 0
View file @ 04d9ed45
This diff is collapsed.
pyramid/phasemapper.py deleted 100644 → 0
View file @ 04d9ed45
This diff is collapsed.
pyramid/projector.py deleted 100644 → 0
View file @ 04d9ed45
This diff is collapsed.
pyramid/reconstruction.py deleted 100644 → 0
View file @ 04d9ed45
This diff is collapsed.
pyramid/regularisator.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 18 09:24:58 2014
@author: Jan
""" # TODO: Docstring!
import abc
import numpy as np
from scipy.sparse import coo_matrix, csr_matrix
import jutil.norms as jnorm
from pyramid.converter import IndexConverter
import logging
# TODO: Fragen für Jörn: Macht es Sinn, x_a standardmäßig auf den Nullvektor zu setzen? Ansonsten
# besser im jeweiligen Konstruktor setzen, nicht im abstrakten!
# Wie kommt man genau an die Ableitungen (Normen sind nicht unproblematisch)?
# Wie implementiert man am besten verschiedene Normen?
class Regularisator(object):
# TODO: Docstring!
__metaclass__ = abc.ABCMeta
LOG = logging.getLogger(__name__+'.Regularisator')
@abc.abstractmethod
def __init__(self, norm, lam):
self.LOG.debug('Calling __init__')
self.norm = norm
self.lam = lam
self.LOG.debug('Created '+str(self))
def __call__(self, x):
self.LOG.debug('Calling __call__')
return self.lam * self.norm(x)
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(norm=%r, lam=%r)' % (self.__class__, self.norm, self.lam)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'Regularisator(norm=%s, lam=%s)' % (self.norm, self.lam)
def jac(self, x):
# TODO: Docstring!
self.LOG.debug('Calling jac')
return self.lam * self.norm.jac_dot(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
self.LOG.debug('Calling hess_dot')
return self.lam * self.norm.hess_dot(x, vector)
def hess_diag(self, x, vector):
# TODO: Docstring!
self.LOG.debug('Calling hess_diag')
return self.lam * self.norm.hess_diag(x, vector)
class NoneRegularisator(Regularisator):
# TODO: Docstring
# TODO: Necessary class? Use others with lam=0?
LOG = logging.getLogger(__name__+'.NoneRegularisator')
def __init__(self):
self.LOG.debug('Calling __init__')
self.norm = None
self.lam = 0
self.LOG.debug('Created '+str(self))
def __call__(self, x):
self.LOG.debug('Calling __call__')
return 0
def jac(self, x):
# TODO: Docstring!
self.LOG.debug('Calling jac')
return np.zeros_like(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
self.LOG.debug('Calling hess_dot')
return np.zeros_like(vector)
def hess_diag(self, x, vector):
# TODO: Docstring!
self.LOG.debug('Calling hess_diag')
return np.zeros_like(vector)
class ZeroOrderRegularisator(Regularisator):
# TODO: Docstring!
LOG = logging.getLogger(__name__+'.ZeroOrderRegularisator')
def __init__(self, lam):
self.LOG.debug('Calling __init__')
norm = jnorm.NormL2Square()
super(ZeroOrderRegularisator, self).__init__(norm, lam)
self.LOG.debug('Created '+str(self))
#class FirstOrderRegularisator(Regularisator):
# # TODO: Docstring!
#
# def __init__(self, fwd_model, lam, x_a=None):
# size_3d = fwd_model.size_3d
# dim = fwd_model.dim
# converter = IndexConverter(dim)
# row = []
# col = []
# data = []
#
# for i in range(size_3d):
# neighbours = converter.find_neighbour_ind(i)
#
# Sa_inv = csr_matrix(coo_matrix(data, (rows, columns)), shape=(3*size_3d, 3*size_3d))
#
# term2 = []
# for i in range(3):
# component = mag_data[i, ...]
# for j in range(3):
# if component.shape[j] > 1:
# term2.append(np.diff(component, axis=j).reshape(-1))
#
# super(FirstOrderRegularisator, self).__init__(Sa_inv, x_a)
# self.LOG.debug('Created '+str(self))
pyramid/template.x3d deleted 100644 → 0
View file @ 04d9ed45
This diff is collapsed.
scripts/collaborations/bielefeld.py deleted 100644 → 0
View file @ 04d9ed45
This diff is collapsed.