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
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
# -*- coding: utf-8 -*-
"""This module provides the :class:`~.MagData` class for storing of magnetization data."""
import os
import numpy as np
from numpy.linalg import norm
from scipy.ndimage.interpolation import zoom
import matplotlib.pyplot as plt
import matplotlib.cm as cmx
from matplotlib.ticker import MaxNLocator
from mayavi import mlab
from lxml import etree
from numbers import Number
import netCDF4
import logging
class MagData(object):
'''Class for storing magnetization data.
Represents 3-dimensional magnetic distributions with 3 components which are stored as a
2-dimensional numpy array in `magnitude`, but which can also be accessed as a vector via
`mag_vec`. :class:`~.MagData` objects support negation, arithmetic operators
(``+``, ``-``, ``*``) and their augmented counterparts (``+=``, ``-=``, ``*=``), with numbers
and other :class:`~.MagData` objects, if their dimensions and grid spacings match. It is
possible to load data from NetCDF4 or LLG (.txt) files or to save the data in these formats.
Plotting methods are also provided.
Attributes
----------
a: float
The grid spacing in nm.
dim: tuple (N=3)
Dimensions (z, y, x) of the grid.
magnitude: :class:`~numpy.ndarray` (N=4)
The `x`-, `y`- and `z`-component of the magnetization vector for every 3D-gridpoint
as a 4-dimensional numpy array (first dimension has to be 3, because of the 3 components).
mag_vec: :class:`~numpy.ndarray` (N=1)
Vector containing the magnetic distribution.
'''
LOG = logging.getLogger(__name__+'.MagData')
@property
def a(self):
return self._a
@a.setter
def a(self, a):
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
self._a = float(a)
@property
def dim(self):
return self._dim
@property
def magnitude(self):
return self._magnitude
@magnitude.setter
def magnitude(self, magnitude):
assert isinstance(magnitude, np.ndarray), 'Magnitude has to be a numpy array!'
assert len(magnitude.shape) == 4, 'Magnitude has to be 4-dimensional!'
assert magnitude.shape[0] == 3, 'First dimension of the magnitude has to be 3!'
self._magnitude = magnitude
self._dim = magnitude.shape[1:]
@property
def mag_vec(self):
return np.reshape(self.magnitude, -1)
@mag_vec.setter
def mag_vec(self, mag_vec):
assert isinstance(mag_vec, np.ndarray), 'Vector has to be a numpy array!'
assert np.size(mag_vec) == 3*np.prod(self.dim), 'Vector has to match magnitude dimensions!'
self.magnitude = mag_vec.reshape((3,)+self.dim)
def __init__(self, a, magnitude):
self.LOG.debug('Calling __init__')
self.a = a
self.magnitude = magnitude
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(a=%r, magnitude=%r)' % (self.__class__, self.a, self.magnitude)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'MagData(a=%s, dim=%s)' % (self.a, self.dim)
def __neg__(self): # -self
self.LOG.debug('Calling __neg__')
return MagData(self.a, -self.magnitude)
def __add__(self, other): # self + other
self.LOG.debug('Calling __add__')
assert isinstance(other, (MagData, Number)), \
'Only MagData objects and scalar numbers (as offsets) can be added/subtracted!'
if isinstance(other, MagData):
self.LOG.debug('Adding two MagData objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.magnitude.shape == (3,)+self.dim, \
'Added magnitude has to have the same dimensions!'
return MagData(self.a, self.magnitude+other.magnitude)
else: # other is a Number
self.LOG.debug('Adding an offset')
return MagData(self.a, self.magnitude+other)
def __sub__(self, other): # self - other
self.LOG.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
self.LOG.debug('Calling __mul__')
assert isinstance(other, Number), 'MagData objects can only be multiplied by numbers!'
return MagData(self.a, other*self.magnitude)
def __radd__(self, other): # other + self
self.LOG.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
self.LOG.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
self.LOG.debug('Calling __rmul__')
return self.__mul__(other)
def __iadd__(self, other): # self += other
self.LOG.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
self.LOG.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
self.LOG.debug('Calling __imul__')
return self.__mul__(other)
def copy(self):
'''Returns a copy of the :class:`~.MagData` object
Parameters
----------
None
Returns
-------
mag_data: :class:`~.MagData`
A copy of the :class:`~.MagData`.
'''
self.LOG.debug('Calling copy')
return MagData(self.a, self.magnitude.copy())
def scale_down(self, n=1):
'''Scale down the magnetic distribution by averaging over two pixels along each axis.
Parameters
----------
n : int, optional
Number of times the magnetic distribution is scaled down. The default is 1.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
Only possible, if each axis length is a power of 2!
'''
self.LOG.debug('Calling scale_down')
assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
assert np.all([d % 2**n == 0 for d in self.dim]), 'Dimensions must a multiples of 2!'
self.a = self.a * 2**n
for t in range(n):
# Create coarser grid for the magnetization:
self.magnitude = self.magnitude.reshape(
3, self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2).mean(axis=(6, 4, 2))
def scale_up(self, n=1, order=0):
'''Scale up the magnetic distribution using spline interpolation of the requested order.
Parameters
----------
n : int, optional
Power of 2 with which the grid is scaled. Default is 1, which means every axis is
increased by a factor of ``2**1 = 2``.
order : int, optional
The order of the spline interpolation, which has to be in the range between 0 and 5
and defaults to 0.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and grid spacing accordingly.
'''
self.LOG.debug('Calling scale_up')
assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
assert 5 > order >= 0 and isinstance(order, (int, long)), \
'order must be a positive integer between 0 and 5!'
self.a = self.a / 2**n
self.magnitude = np.array((zoom(self.magnitude[0], zoom=2**n, order=order),
zoom(self.magnitude[1], zoom=2**n, order=order),
zoom(self.magnitude[2], zoom=2**n, order=order)))
def pad(self, x_pad, y_pad, z_pad):
'''Pad the current magnetic distribution with zeros for each individual axis.
Parameters
----------
x_pad : int
Number of zeros which should be padded on both sides of the x-axis.
y_pad : int
Number of zeros which should be padded on both sides of the y-axis.
z_pad : int
Number of zeros which should be padded on both sides of the z-axis.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions accordingly.
'''
assert x_pad >= 0 and isinstance(x_pad, (int, long)), 'x_pad must be a positive integer!'
assert y_pad >= 0 and isinstance(y_pad, (int, long)), 'y_pad must be a positive integer!'
assert z_pad >= 0 and isinstance(z_pad, (int, long)), 'z_pad must be a positive integer!'
self.magnitude = np.pad(self.magnitude,
((0, 0), (z_pad, z_pad), (y_pad, y_pad), (x_pad, x_pad)),
mode='constant', constant_values=0)
def get_mask(self, threshold=0):
'''Mask all pixels where the amplitude of the magnetization lies above `threshold`.
Parameters
----------
threshold : float, optional
A pixel only gets masked, if it lies above this threshold . The default is 0.
Returns
-------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Mask of the pixels where the amplitude of the magnetization lies above `threshold`.
'''
return np.sqrt(np.sum(np.array(self.magnitude)**2, axis=0)) > threshold
def get_vector(self, mask):
'''Returns the magnetic components arranged in a vector, specified by a mask.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Masks the pixels from which the components should be taken.
Returns
-------
vector : :class:`~numpy.ndarray` (N=1)
The vector containing magnetization components of the specified pixels.
Order is: first all `x`-, then all `y`-, then all `z`-components.
'''
return np.reshape([self.magnitude[2][mask],
self.magnitude[1][mask],
self.magnitude[0][mask]], -1)
def set_vector(self, vector, mask=None):
'''Set the magnetic components of the masked pixels to the values specified by `vector`.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean), optional
Masks the pixels from which the components should be taken.
vector : :class:`~numpy.ndarray` (N=1)
The vector containing magnetization components of the specified pixels.
Order is: first all `x`-, then all `y-, then all `z`-components.
Returns
-------
None
'''
assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
count = np.size(vector)/3
if mask is not None:
self.magnitude[2][mask] = vector[:count] # x-component
self.magnitude[1][mask] = vector[count:2*count] # y-component
self.magnitude[0][mask] = vector[2*count:] # z-component
else:
self.mag_vec = vector
def save_to_llg(self, filename='..\output\magdata_output.txt'):
'''Save magnetization data in a file with LLG-format.
Parameters
----------
filename : string, optional
The name of the LLG-file in which to store the magnetization data.
The default is '..\output\magdata_output.txt'.
Returns
-------
None
'''
self.LOG.debug('Calling save_to_llg')
a = self.a * 1.0E-9 / 1.0E-2 # from nm to cm
# Create 3D meshgrid and reshape it and the magnetization into a list where x varies first:
zz, yy, xx = np.mgrid[a/2:(self.dim[0]*a-a/2):self.dim[0]*1j,
a/2:(self.dim[1]*a-a/2):self.dim[1]*1j,
a/2:(self.dim[2]*a-a/2):self.dim[2]*1j].reshape(3, -1)
x_vec, y_vec, z_vec = self.magnitude.reshape(3, -1)
# Save data to file:
data = np.array([xx, yy, zz, x_vec, y_vec, z_vec]).T
with open(filename, 'w') as mag_file:
mag_file.write('LLGFileCreator: %s\n' % filename.replace('.txt', ''))
mag_file.write(' %d %d %d\n' % (self.dim[2], self.dim[1], self.dim[0]))
mag_file.writelines('\n'.join(' '.join('{:7.6e}'.format(cell)
for cell in row) for row in data))
@classmethod
def load_from_llg(cls, filename):
'''Construct :class:`~.MagData` object from LLG-file.
Parameters
----------
filename : string
The name of the LLG-file from which to load the data.
Returns
-------
mag_data: :class:`~.MagData`
A :class:`~.MagData` object containing the loaded data.
'''
cls.LOG.debug('Calling load_from_llg')
SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
data = np.genfromtxt(filename, skip_header=2)
dim = tuple(np.genfromtxt(filename, dtype=int, skip_header=1, skip_footer=len(data[:, 0])))
a = (data[1, 0] - data[0, 0]) / SCALE
magnitude = data[:, 3:6].T.reshape((3,)+dim)
return MagData(a, magnitude)
def save_to_netcdf4(self, filename='..\output\magdata_output.nc'):
'''Save magnetization data in a file with NetCDF4-format.
Parameters
----------
filename : string, optional
The name of the NetCDF4-file in which to store the magnetization data.
The default is '..\output\magdata_output.nc'.
Returns
-------
None
'''
self.LOG.debug('Calling save_to_netcdf4')
mag_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
mag_file.a = self.a
mag_file.createDimension('comp', 3) # Number of components
mag_file.createDimension('z_dim', self.dim[0])
mag_file.createDimension('y_dim', self.dim[1])
mag_file.createDimension('x_dim', self.dim[2])
magnitude = mag_file.createVariable('magnitude', 'f', ('comp', 'z_dim', 'y_dim', 'x_dim'))
magnitude[...] = self.magnitude
mag_file.close()
@classmethod
def load_from_netcdf4(cls, filename):
'''Construct :class:`~.DataMag` object from NetCDF4-file.
Parameters
----------
filename : string
The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
Returns
-------
mag_data: :class:`~.MagData`
A :class:`~.MagData` object containing the loaded data.
'''
cls.LOG.debug('Calling copy')
mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = mag_file.a
magnitude = mag_file.variables['magnitude'][...]
mag_file.close()
return MagData(a, magnitude)
def quiver_plot(self, title='Magnetization Distribution', axis=None, proj_axis='z',
ax_slice=None, log=False, scaled=True, show=True):
'''Plot a slice of the magnetization as a quiver plot.
Parameters
----------
title : string, optional
The title for the plot.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
proj_axis : {'z', 'y', 'x'}, optional
The axis, from which a slice is plotted. The default is 'z'.
ax_slice : int, optional
The slice-index of the axis specified in `proj_axis`. Is set to the center of
`proj_axis` if not specified.
log : boolean, optional
Takes the Default is False.
scaled : boolean, optional
Normalizes the plotted arrows in respect to the highest one. Default is True.
show: bool, optional
A switch which determines if the plot is shown at the end of plotting. Default is True.
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
'''
self.LOG.debug('Calling quiver_plot')
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
self.LOG.debug('proj_axis == z')
if ax_slice is None:
self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[0]/2.)
plot_u = np.copy(self.magnitude[0][ax_slice, ...]) # x-component
plot_v = np.copy(self.magnitude[1][ax_slice, ...]) # y-component
u_label = 'x [px]'
v_label = 'y [px]'
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
self.LOG.debug('proj_axis == y')
if ax_slice is None:
self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[1]/2.)
plot_u = np.copy(self.magnitude[0][:, ax_slice, :]) # x-component
plot_v = np.copy(self.magnitude[2][:, ax_slice, :]) # z-component
u_label = 'x [px]'
v_label = 'z [px]'
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
self.LOG.debug('proj_axis == x')
if ax_slice is None:
self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[2]/2.)
plot_u = np.swapaxes(np.copy(self.magnitude[2][..., ax_slice]), 0, 1) # z-component
plot_v = np.swapaxes(np.copy(self.magnitude[1][..., ax_slice]), 0, 1) # y-component
u_label = 'z [px]'
v_label = 'y [px]'
# If no axis is specified, a new figure is created:
if axis is None:
self.LOG.debug('axis is None')
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
angles = np.angle(plot_u+1j*plot_v, deg=True)
# Take the logarithm of the arrows to clearly show directions (if specified):
if log:
cutoff = 10
amp = np.round(np.hypot(plot_u, plot_v), decimals=cutoff)
min_value = amp[np.nonzero(amp)].min()
plot_u = np.round(plot_u, decimals=cutoff) / min_value
plot_u = np.log10(np.abs(plot_u)+1) * np.sign(plot_u)
plot_v = np.round(plot_v, decimals=cutoff) / min_value
plot_v = np.log10(np.abs(plot_v)+1) * np.sign(plot_v)
# Scale the magnitude of the arrows to the highest one (if specified):
if scaled:
plot_u /= np.hypot(plot_u, plot_v).max()
plot_v /= np.hypot(plot_u, plot_v).max()
axis.quiver(plot_u, plot_v, pivot='middle', units='xy', angles=angles,
scale_units='xy', scale=1, headwidth=6, headlength=7)
axis.set_xlim(-1, np.shape(plot_u)[1])
axis.set_ylim(-1, np.shape(plot_u)[0])
axis.set_title(title, fontsize=18)
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
# Show Plot:
if show:
plt.show()
return axis
def quiver_plot3d(self, title='Magnetization Distribution'):
'''Plot the magnetization as 3D-vectors in a quiverplot.
Parameters
----------
None
Returns
-------
None
'''
self.LOG.debug('Calling quiver_plot3D')
a = self.a
dim = self.dim
# Create points and vector components as lists:
zz, yy, xx = np.mgrid[a/2:(dim[0]*a-a/2):dim[0]*1j,
a/2:(dim[1]*a-a/2):dim[1]*1j,
a/2:(dim[2]*a-a/2):dim[2]*1j]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
x_mag = np.reshape(self.magnitude[0], (-1))
y_mag = np.reshape(self.magnitude[1], (-1))
z_mag = np.reshape(self.magnitude[2], (-1))
# Plot them as vectors:
mlab.figure()
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow')
mlab.outline(plot)
mlab.axes(plot)
mlab.title(title, height=0.95, size=0.35)
mlab.colorbar(None, label_fmt='%.2f')
mlab.colorbar(None, orientation='vertical')
def save_to_x3d(self, filename='..\..\output\magdata_output.x3d', maximum=1):
'''Output the magnetization in the .x3d format for the Fraunhofer InstantReality Player.
Parameters
----------
None
Returns
-------
None
'''
self.LOG.debug('Calling save_to_x3d')
dim = self.dim
# Create points and vector components as lists:
zz, yy, xx = np.mgrid[0.5:(dim[0]-0.5):dim[0]*1j,
0.5:(dim[1]-0.5):dim[1]*1j,
0.5:(dim[2]-0.5):dim[2]*1j]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
x_mag = np.reshape(self.magnitude[0], (-1))
y_mag = np.reshape(self.magnitude[1], (-1))
z_mag = np.reshape(self.magnitude[2], (-1))
# Load template, load tree and write viewpoint information:
template = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'template.x3d')
parser = etree.XMLParser(remove_blank_text=True)
tree = etree.parse(template, parser)
scene = tree.find('Scene')
etree.SubElement(scene, 'Viewpoint', position='0 0 {}'.format(1.5*dim[0]),
fieldOfView='1')
# Write each "spin"-tag separately:
for i in range(np.prod(dim)):
mag = np.sqrt(x_mag[i]**2+y_mag[i]**2+z_mag[i]**2)
if mag != 0:
spin_position = (xx[i]-dim[2]/2., yy[i]-dim[1]/2., zz[i]-dim[0]/2.)
sx_ref = 0
sy_ref = 1
sz_ref = 0
rot_x = sy_ref*z_mag[i] - sz_ref*y_mag[i]
rot_y = sz_ref*x_mag[i] - sx_ref*z_mag[i]
rot_z = sx_ref*y_mag[i] - sy_ref*x_mag[i]
angle = np.arccos(y_mag[i]/mag)
if norm((rot_x, rot_y, rot_z)) < 1E-10:
rot_x, rot_y, rot_z = 1, 0, 0
spin_rotation = (rot_x, rot_y, rot_z, angle)
spin_color = cmx.RdYlGn(mag/maximum)[:3]
spin_scale = (1., 1., 1.)
spin = etree.SubElement(scene, 'ProtoInstance',
DEF='Spin {}'.format(i), name='Spin_Proto')
etree.SubElement(spin, 'fieldValue', name='spin_position',
value='{} {} {}'.format(*spin_position))
etree.SubElement(spin, 'fieldValue', name='spin_rotation',
value='{} {} {} {}'.format(*spin_rotation))
etree.SubElement(spin, 'fieldValue', name='spin_color',
value='{} {} {}'.format(*spin_color))
etree.SubElement(spin, 'fieldValue', name='spin_scale',
value='{} {} {}'.format(*spin_scale))
# Write the tree into the file in pretty print format:
tree.write(filename, pretty_print=True)
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
# -*- coding: utf-8 -*-
"""This module provides the :class:`~.PhaseMap` class for storing phase map data."""
import numpy as np
from numpy import pi
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import NullLocator, MaxNLocator, FuncFormatter
from PIL import Image
from numbers import Number
import netCDF4
import logging
class PhaseMap(object):
'''Class for storing phase map data.
Represents 2-dimensional phase maps. The phase information itself is stored as a 2-dimensional
matrix in `phase`, but can also be accessed as a vector via `phase_vec`. :class:`~.PhaseMap`
objects support negation, arithmetic operators (``+``, ``-``, ``*``) and their augmented
counterparts (``+=``, ``-=``, ``*=``), with numbers and other :class:`~.PhaseMap`
objects, if their dimensions and grid spacings match. It is possible to load data from NetCDF4
or textfiles or to save the data in these formats. Methods for plotting the phase or a
corresponding holographic contour map are provided. Holographic contour maps are created by
taking the cosine of the (optionally amplified) phase and encoding the direction of the
2-dimensional gradient via color. The directional encoding can be seen by using the
:func:`~.make_color_wheel` function. Use the :func:`~.display_combined` function to plot the
phase map and the holographic contour map next to each other.
Attributes
----------
a: float
The grid spacing in nm.
dim_uv: tuple (N=2)
Dimensions of the grid.
phase: :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
phase_vec: :class:`~numpy.ndarray` (N=2)
Vector containing the phase shift.
unit: {'rad', 'mrad'}, optional
Set the unit of the phase map. This is important for the :func:`display` function,
because the phase is scaled accordingly. Does not change the phase itself, which is
always in `rad`.
'''
LOG = logging.getLogger(__name__)
UNITDICT = {u'rad': 1E0,
u'mrad': 1E3,
u'µrad': 1E6}
CDICT = {'red': [(0.00, 1.0, 0.0),
(0.25, 1.0, 1.0),
(0.50, 1.0, 1.0),
(0.75, 0.0, 0.0),
(1.00, 0.0, 1.0)],
'green': [(0.00, 0.0, 0.0),
(0.25, 0.0, 0.0),
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 1.0)],
'blue': [(0.00, 1.0, 1.0),
(0.25, 0.0, 0.0),
(0.50, 0.0, 0.0),
(0.75, 0.0, 0.0),
(1.00, 1.0, 1.0)]}
CDICT_INV = {'red': [(0.00, 0.0, 1.0),
(0.25, 0.0, 0.0),
(0.50, 0.0, 0.0),
(0.75, 1.0, 1.0),
(1.00, 1.0, 0.0)],
'green': [(0.00, 1.0, 1.0),
(0.25, 1.0, 1.0),
(0.50, 0.0, 0.0),
(0.75, 0.0, 0.0),
(1.00, 1.0, 0.0)],
'blue': [(0.00, 0.0, 0.0),
(0.25, 1.0, 1.0),
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 0.0)]}
HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
HOLO_CMAP_INV = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT_INV, 256)
@property
def a(self):
return self._a
@a.setter
def a(self, a):
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
self._a = float(a)
@property
def dim_uv(self):
return self._dim_uv
@property
def phase(self):
return self._phase
@phase.setter
def phase(self, phase):
assert isinstance(phase, np.ndarray), 'Phase has to be a numpy array!'
assert len(phase.shape) == 2, 'Phase has to be 2-dimensional!'
self._phase = phase
self._dim_uv = phase.shape
@property
def phase_vec(self):
return np.reshape(self.phase, -1)
@phase_vec.setter
def phase_vec(self, phase_vec):
assert isinstance(phase_vec, np.ndarray), 'Vector has to be a numpy array!'
assert np.size(phase_vec) == np.prod(self.dim_uv), 'Vector size has to match phase!'
self.phase = phase_vec.reshape(self.dim_uv)
@property
def unit(self):
return self._unit
@unit.setter
def unit(self, unit):
assert unit in self.UNITDICT, 'Unit not supported!'
self._unit = unit
def __init__(self, a, phase, unit='rad'):
self.LOG.debug('Calling __init__')
self.a = a
self.phase = phase
self.unit = unit
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(a=%r, phase=%r, unit=%r)' % \
(self.__class__, self.a, self.phase, self.unit)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'PhaseMap(a=%s, dim_uv=%s)' % (self.a, self.dim_uv)
def __neg__(self): # -self
self.LOG.debug('Calling __neg__')
return PhaseMap(self.a, -self.phase, self.unit)
def __add__(self, other): # self + other
self.LOG.debug('Calling __add__')
assert isinstance(other, (PhaseMap, Number)), \
'Only PhaseMap objects and scalar numbers (as offsets) can be added/subtracted!'
if isinstance(other, PhaseMap):
self.LOG.debug('Adding two PhaseMap objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.phase.shape == self.dim_uv, \
'Added magnitude has to have the same dimensions!'
return PhaseMap(self.a, self.phase+other.phase, self.unit)
else: # other is a Number
self.LOG.debug('Adding an offset')
return PhaseMap(self.a, self.phase+other, self.unit)
def __sub__(self, other): # self - other
self.LOG.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
self.LOG.debug('Calling __mul__')
assert (isinstance(other, Number)
or (isinstance(other, np.ndarray) and other.shape == self.dim_uv)), \
'PhaseMap objects can only be multiplied by scalar numbers or fitting arrays!'
return PhaseMap(self.a, other*self.phase, self.unit)
def __radd__(self, other): # other + self
self.LOG.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
self.LOG.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
self.LOG.debug('Calling __rmul__')
return self.__mul__(other)
def __iadd__(self, other): # self += other
self.LOG.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
self.LOG.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
self.LOG.debug('Calling __imul__')
return self.__mul__(other)
def save_to_txt(self, filename='..\output\phasemap_output.txt'):
'''Save :class:`~.PhaseMap` data in a file with txt-format.
Parameters
----------
filename : string
The name of the file in which to store the phase map data.
The default is '..\output\phasemap_output.txt'.
Returns
-------
None
'''
self.LOG.debug('Calling save_to_txt')
with open(filename, 'w') as phase_file:
phase_file.write('{}\n'.format(filename.replace('.txt', '')))
phase_file.write('grid spacing = {} nm\n'.format(self.a))
np.savetxt(phase_file, self.phase, fmt='%7.6e', delimiter='\t')
@classmethod
def load_from_txt(cls, filename):
'''Construct :class:`~.PhaseMap` object from a human readable txt-file.
Parameters
----------
filename : string
The name of the file from which to load the data.
Returns
-------
phase_map : :class:`~.PhaseMap`
A :class:`~.PhaseMap` object containing the loaded data.
'''
cls.LOG.debug('Calling load_from_txt')
with open(filename, 'r') as phase_file:
phase_file.readline() # Headerline is not used
a = float(phase_file.readline()[15:-4])
phase = np.loadtxt(filename, delimiter='\t', skiprows=2)
return PhaseMap(a, phase)
def save_to_netcdf4(self, filename='..\output\phasemap_output.nc'):
'''Save :class:`~.PhaseMap` data in a file with NetCDF4-format.
Parameters
----------
filename : string, optional
The name of the NetCDF4-file in which to store the phase data.
The default is '..\output\phasemap_output.nc'.
Returns
-------
None
'''
self.LOG.debug('Calling save_to_netcdf4')
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
phase_file.a = self.a
phase_file.createDimension('v_dim', self.dim_uv[0])
phase_file.createDimension('u_dim', self.dim_uv[1])
phase = phase_file.createVariable('phase', 'f', ('v_dim', 'u_dim'))
phase[:] = self.phase
phase_file.close()
@classmethod
def load_from_netcdf4(cls, filename):
'''Construct :class:`~.PhaseMap` object from NetCDF4-file.
Parameters
----------
filename : string
The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
Returns
-------
phase_map: :class:`~.PhaseMap`
A :class:`~.PhaseMap` object containing the loaded data.
'''
cls.LOG.debug('Calling load_from_netcdf4')
phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = phase_file.a
phase = phase_file.variables['phase'][:]
phase_file.close()
return PhaseMap(a, phase)
def display_phase(self, title='Phase Map', cmap='RdBu', limit=None,
norm=None, axis=None, cbar=True, show=True):
'''Display the phasemap as a colormesh.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Phase Map'.
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.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
cbar : bool, optional
A switch determining if the colorbar should be plotted or not. Default is True.
show : bool, optional
A switch which determines if the plot is shown at the end of plotting.
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
'''
self.LOG.debug('Calling display_phase')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
if limit is None:
limit = np.max(np.abs(phase))
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Plot the phasemap:
im = axis.pcolormesh(phase, cmap=cmap, vmin=-limit, vmax=limit, norm=norm)
# Set the axes ticks and labels:
axis.set_xlim(0, self.dim_uv[1])
axis.set_ylim(0, self.dim_uv[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('u-axis [nm]', fontsize=15)
axis.set_ylabel('v-axis [nm]', fontsize=15)
# Add colorbar:
if cbar:
fig = plt.gcf()
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.82, 0.15, 0.02, 0.7])
cbar = fig.colorbar(im, cax=cbar_ax)
cbar.ax.tick_params(labelsize=14)
cbar.set_label(u'phase shift [{}]'.format(self.unit), fontsize=15)
# Show plot:
if show:
plt.show()
# Return plotting axis:
return axis
def display_phase3d(self, title='Phase Map', cmap='RdBu', show=True):
'''Display the phasemap as a 3-D surface with contourplots.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Phase Map'.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
show: bool, optional
A switch which determines if the plot is shown at the end of plotting. Default is True.
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
'''
self.LOG.debug('Calling display_phase3d')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
# Create figure and axis:
fig = plt.figure()
axis = Axes3D(fig)
# Plot surface and contours:
u = range(self.dim_uv[1])
v = range(self.dim_uv[0])
uu, vv = np.meshgrid(u, v)
axis.plot_surface(uu, vv, phase, rstride=4, cstride=4, alpha=0.7, cmap=cmap,
linewidth=0, antialiased=False)
axis.contourf(uu, vv, phase, 15, zdir='z', offset=np.min(phase), cmap=cmap)
axis.view_init(45, -135)
axis.set_xlabel('u-axis [px]')
axis.set_ylabel('v-axis [px]')
axis.set_zlabel('phase shift [{}]'.format(self.unit))
# Show Plot:
if show:
plt.show()
# Return plotting axis:
return axis
def display_holo(self, title=None, gain='auto', axis=None, grad_encode='bright',
interpolation='none', show=True):
'''Display the color coded holography image.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Contour Map (gain: %g)' % gain.
gain : float or 'auto', optional
The gain factor for determining the number of contour lines. The default is 'auto',
which means that the gain will be determined automatically to look pretty.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
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.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method. No interpolation is used in the default case.
show: bool, optional
A switch which determines if the plot is shown at the end of plotting. Default is True.
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
'''
self.LOG.debug('Calling display_holo')
# Calculate gain if 'auto' is selected:
if gain == 'auto':
gain = 4 * 2*pi/self.phase.max()
# Set title if not set:
if title is None:
title = 'Contour Map (gain: %.2g)' % gain
# Calculate the holography image intensity:
img_holo = (1 + np.cos(gain * self.phase)) / 2
# Calculate the phase gradients, expressed by magnitude and angle:
phase_grad_y, phase_grad_x = np.gradient(self.phase, self.a, self.a)
phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)
if phase_magnitude.max() != 0: # Take care of phase maps with only zeros
saturation = np.sin(phase_magnitude/phase_magnitude.max() * pi / 2)
phase_saturation = np.dstack((saturation,)*4)
# Color code the angle and create the holography image:
if grad_encode == 'dark':
self.LOG.debug('gradient encoding: dark')
rgba = self.HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * saturation.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'bright':
self.LOG.debug('gradient encoding: bright')
rgba = self.HOLO_CMAP(phase_angle)+(1-phase_saturation)*self.HOLO_CMAP_INV(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'color':
self.LOG.debug('gradient encoding: color')
rgba = self.HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'none':
self.LOG.debug('gradient encoding: none')
rgba = self.HOLO_CMAP(phase_angle)+self.HOLO_CMAP_INV(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
else:
raise AssertionError('Gradient encoding not recognized!')
holo_image = Image.fromarray(rgb)
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Plot the image and set axes:
axis.imshow(holo_image, origin='lower', interpolation=interpolation)
# Set the title and the axes labels:
axis.set_title(title)
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('u-axis [px]', fontsize=15)
axis.set_ylabel('v-axis [px]', fontsize=15)
axis.set_xlim(0, self.dim_uv[1])
axis.set_ylim(0, self.dim_uv[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
# Show Plot:
if show:
plt.show()
# Return plotting axis:
return axis
def display_combined(self, title='Combined Plot', cmap='RdBu', limit=None, norm=None,
gain='auto', interpolation='none', grad_encode='bright', show=True):
'''Display the phase map and the resulting color coded holography image in one plot.
Parameters
----------
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 or 'auto', optional
The gain factor for determining the number of contour lines. The default is 'auto',
which means that the gain will be determined automatically to look pretty.
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.
show: bool, optional
A switch which determines if the plot is shown at the end of plotting. Default is True.
Returns
-------
phase_axis, holo_axis: :class:`~matplotlib.axes.AxesSubplot`
The axes on which the graphs are plotted.
'''
self.LOG.debug('Calling display_combined')
# Create combined plot and set title:
fig = plt.figure(figsize=(16, 7))
fig.suptitle(title, fontsize=20)
# Plot holography image:
holo_axis = fig.add_subplot(1, 2, 1, aspect='equal')
self.display_holo(gain=gain, axis=holo_axis, interpolation=interpolation,
show=False, grad_encode=grad_encode)
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
self.display_phase(cmap='RdBu', limit=limit, norm=norm, axis=phase_axis, show=False)
# Show Plot:
if show:
plt.show()
# Return the plotting axes:
return phase_axis, holo_axis
@classmethod
def make_color_wheel(cls, show=True):
'''Display a color wheel to illustrate the color coding of the gradient direction.
Parameters
----------
show: bool, optional
A switch which determines if the plot is shown at the end of plotting. Default is True.
Returns
-------
None
'''
cls.LOG.debug('Calling make_color_wheel')
x = np.linspace(-256, 256, num=512)
y = np.linspace(-256, 256, num=512)
xx, yy = np.meshgrid(x, y)
r = np.sqrt(xx ** 2 + yy ** 2)
# Create the wheel:
color_wheel_magnitude = (1 - np.cos(r * pi/360)) / 2
color_wheel_magnitude *= 0 * (r > 256) + 1 * (r <= 256)
color_wheel_angle = (1 - np.arctan2(xx, -yy)/pi) / 2
# Color code the angle and create the holography image:
rgba = cls.HOLO_CMAP(color_wheel_angle)
rgb = (255.999 * color_wheel_magnitude.T * rgba[:, :, :3].T).T.astype(np.uint8)
color_wheel = Image.fromarray(rgb)
# Plot the color wheel:
fig = plt.figure(figsize=(4, 4))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.imshow(color_wheel, origin='lower')
axis.xaxis.set_major_locator(NullLocator())
axis.yaxis.set_major_locator(NullLocator())
# Show Plot:
if show:
plt.show()
pyramid/phasemapper.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""This module executes several forward models to calculate the magnetic or electric phase map from
a given projection of a 3-dimensional magnetic distribution (see :mod:`~pyramid.projector`).
For the magnetic phase map, an approach using real space and one using Fourier space is provided.
The electrostatic contribution is calculated by using the assumption of a mean inner potential.
"""
import numpy as np
from numpy import pi
import abc
import pyramid.numcore.phasemapper_core as nc
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.projector import SimpleProjector
from pyramid.kernel import Kernel
import logging
class PhaseMapper(object):
'''Abstract base class for the phase calculation from a 2-dimensional distribution.
The :class:`~.PhaseMapper-` class represents a strategy for the phasemapping of a
2-dimensional magnetic distribution with two components onto a scalar phase map.
:class:`~.Kernel` is an abstract base class and provides a unified interface which should be
subclassed with custom :func:`__init__` and :func:`__call__` functions. Concrete subclasses
can be called as a function and take a :class:`~.MagData` object as input and return a
:class:`~.PhaseMap` object.
'''
__metaclass__ = abc.ABCMeta
LOG = logging.getLogger(__name__+'.PhaseMapper')
@abc.abstractmethod
def __call__(self, mag_data):
self.LOG.debug('Calling __call__')
raise NotImplementedError()
@abc.abstractmethod
def jac_dot(self, vector):
self.LOG.debug('Calling jac_dot')
raise NotImplementedError()
@abc.abstractmethod
def jac_T_dot(self, vector):
self.LOG.debug('Calling jac_T_dot')
raise NotImplementedError()
class PhaseMapperRDFC(PhaseMapper):
'''Class representing a phase mapping strategy using real space discretization and Fourier
space convolution.
The :class:`~.PMConvolve` class represents a phase mapping strategy involving discretization in
real space. It utilizes the convolution in Fourier space, directly takes :class:`~.MagData`
objects and returns :class:`~.PhaseMap` objects.
Attributes
----------
a : float
The grid spacing in nm.
projector : :class:`~.Projector`
Projector which should be used for the projection of the 3-dimensional magnetization
distribution.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
geometry : {'disc', 'slab'}, optional
Elementary geometry which is used for the phase contribution of one pixel.
Default is 'disc'.
n: int
Size of the image space.
m: int
Size of the input space.
'''
LOG = logging.getLogger(__name__+'.PMConvolve')
def __init__(self, kernel):
self.LOG.debug('Calling __init__')
self.kernel = kernel
self.n = np.prod(kernel.dim_uv)
self.m = 2 * self.n
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(kernel=%r)' % (self.__class__, self.kernel)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'PhaseMapperRDFC(kernel=%s)' % (self.kernel)
def __call__(self, mag_data):
self.LOG.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.kernel.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert mag_data.dim[1:3] == self.kernel.dim_uv, 'Dimensions do not match!'
# Process input parameters:
u_mag, v_mag = mag_data.magnitude[0:2, 0, ...]
return PhaseMap(mag_data.a, self._convolve(u_mag, v_mag))
def _convolve(self, u_mag, v_mag):
# Fourier transform the projected magnetisation:
u_mag_fft = np.fft.rfftn(u_mag, self.kernel.dim_fft)
v_mag_fft = np.fft.rfftn(v_mag, self.kernel.dim_fft)
# Convolve the magnetization with the kernel in Fourier space:
phase_fft = u_mag_fft*self.kernel.u_fft - v_mag_fft*self.kernel.v_fft
# Return the result:
return np.fft.irfftn(phase_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
(row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
``N**2`` elements to the `v`-component of the magnetization.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
self.LOG.debug('Calling jac_dot')
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
u_mag, v_mag = np.reshape(vector, (2,)+self.kernel.dim_uv)
return self._convolve(u_mag, v_mag)
def jac_T_dot(self, vector):
'''Calculate the product of the transposed Jacobi matrix with a given `vector`.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vector with ``N**2`` entries which represents a matrix with dimensions like a scalar
phasemap.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the transposed Jacobi matrix (which is not explicitely calculated) with
the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
'''
self.LOG.debug('Calling jac_T_dot')
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
# TODO: directly? ask Jörn again!
result = np.zeros(self.m)
nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1],
self.kernel.u, self.kernel.v, vector, result)
return result
class PhaseMapperRDRC(PhaseMapper):
'''Class representing a phase mapping strategy using real space discretization.
The :class:`~.PMReal` class represents a phase mapping strategy involving discretization in
real space. It directly takes :class:`~.MagData` objects and returns :class:`~.PhaseMap`
objects.
Attributes
----------
a : float
The grid spacing in nm.
projector : :class:`~.Projector`
Projector which should be used for the projection of the 3-dimensional magnetization
distribution.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
threshold : float, optional
Threshold determining for which voxels the phase contribution will be calculated. The
default is 0, which means that all voxels with non-zero magnetization will contribute.
Should be above noise level.
geometry : {'disc', 'slab'}, optional
Elementary geometry which is used for the phase contribution of one pixel.
Default is 'disc'.
numcore : boolean, optional
Boolean choosing if Cython enhanced routines from the :mod:`~.pyramid.numcore` module
should be used. Default is True.
'''
LOG = logging.getLogger(__name__+'.PMReal')
def __init__(self, kernel, threshold=0, numcore=True):
self.LOG.debug('Calling __init__')
self.kernel = kernel
self.threshold = threshold
self.numcore = numcore
self.n = np.prod(kernel.dim_uv)
self.m = 2 * self.n
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(kernel=%r, threshold=%r, numcore=%r)' % \
(self.__class__, self.kernel, self.threshold, self.numcore)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'PhaseMapperRDRC(kernel=%s, threshold=%s, numcore=%s)' % \
(self.kernel, self.threshold, self.numcore)
def __call__(self, mag_data):
self.LOG.debug('Calling __call__')
dim_uv = self.kernel.dim_uv
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.kernel.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert mag_data.dim[1:3] == dim_uv, 'Dimensions do not match!'
# Process input parameters:
u_mag, v_mag = mag_data.magnitude[0:2, 0, ...]
# Get kernel (lookup-tables for the phase of one pixel):
u_phi = self.kernel.u
v_phi = self.kernel.v
# Calculation of the phase:
phase = np.zeros(dim_uv)
if self.numcore:
nc.phasemapper_real_convolve(dim_uv[0], dim_uv[1], v_phi, u_phi,
v_mag, u_mag, phase, self.threshold)
else:
for j in range(dim_uv[0]):
for i in range(dim_uv[1]):
u_phase = u_phi[dim_uv[0]-1-j:(2*dim_uv[0]-1)-j,
dim_uv[1]-1-i:(2*dim_uv[1]-1)-i]
if abs(u_mag[j, i]) > self.threshold:
phase += u_mag[j, i] * u_phase
v_phase = v_phi[dim_uv[0]-1-j:(2*dim_uv[0]-1)-j,
dim_uv[1]-1-i:(2*dim_uv[1]-1)-i]
if abs(v_mag[j, i]) > self.threshold:
phase -= v_mag[j, i] * v_phase
# Return the phase:
return PhaseMap(mag_data.a, phase)
def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
(row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
``N**2`` elements to the `v`-component of the magnetization.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
self.LOG.debug('Calling jac_dot')
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
v_dim, u_dim = self.kernel.dim_uv
result = np.zeros(self.n)
if self.numcore:
nc.jac_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result)
else:
# Iterate over all contributing pixels (numbered consecutively)
for s in range(self.n): # column-wise (two columns at a time, u- and v-component)
i = s % u_dim # u-coordinate of current contributing pixel
j = int(s/u_dim) # v-coordinate of current ccontributing pixel
u_min = (u_dim-1) - i # u_dim-1: center of the kernel
u_max = (2*u_dim-1) - i # = u_min + u_dim
v_min = (v_dim-1) - j # v_dim-1: center of the kernel
v_max = (2*v_dim-1) - j # = v_min + v_dim
result += vector[s] * self.kernel.u[v_min:v_max, u_min:u_max].reshape(-1)
result -= vector[s+self.n] * self.kernel.v[v_min:v_max, u_min:u_max].reshape(-1)
return result
def jac_T_dot(self, vector):
'''Calculate the product of the transposed Jacobi matrix with a given `vector`.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vector with ``N**2`` entries which represents a matrix with dimensions like a scalar
phasemap.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the transposed Jacobi matrix (which is not explicitely calculated) with
the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
'''
self.LOG.debug('Calling jac_T_dot')
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
v_dim, u_dim = self.dim_uv
result = np.zeros(self.m)
if self.numcore:
nc.jac_T_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result)
else:
# Iterate over all contributing pixels (numbered consecutively):
for s in range(self.n): # row-wise (two rows at a time, u- and v-component)
i = s % u_dim # u-coordinate of current contributing pixel
j = int(s/u_dim) # v-coordinate of current contributing pixel
u_min = (u_dim-1) - i # u_dim-1: center of the kernel
u_max = (2*u_dim-1) - i # = u_min + u_dim
v_min = (v_dim-1) - j # v_dim-1: center of the kernel
v_max = (2*v_dim-1) - j # = v_min + v_dim
result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1))
result[s+self.n] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))
return result
class PhaseMapperFDFC(PhaseMapper):
'''Class representing a phase mapping strategy using a discretization in Fourier space.
The :class:`~.PMFourier` class represents a phase mapping strategy involving discretization in
Fourier space. It utilizes the Fourier transforms, which are inherently calculated in the
:class:`~.Kernel` class and directly takes :class:`~.MagData` objects and returns
:class:`~.PhaseMap` objects.
Attributes
----------
a : float
The grid spacing in nm.
projector : :class:`~.Projector`
Projector which should be used for the projection of the 3-dimensional magnetization
distribution.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
padding : int, optional
Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
'''
LOG = logging.getLogger(__name__+'.PMFourier')
PHI_0 = -2067.83 # magnetic flux in T*nm²
def __init__(self, a, dim_uv, b_0=1, padding=0):
self.LOG.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.b_0 = b_0
self.padding = padding
self.n = np.prod(dim_uv)
self.m = 2 * self.n
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, b_0=%r, padding=%r)' % \
(self.__class__, self.a, self.dim_uv, self.b_0, self.padding)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'PhaseMapperFDFC(a=%s, dim_uv=%s, b_0=%s, padding=%s)' % \
(self.a, self.dim_uv, self.b_0, self.padding)
def __call__(self, mag_data):
self.LOG.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert mag_data.dim[1:3] == self.dim_uv, 'Dimensions do not match!'
v_dim, u_dim = self.dim_uv
u_mag, v_mag = mag_data.magnitude[0:2, 0, ...]
# Create zero padded matrices:
u_pad = int(u_dim/2 * self.padding)
v_pad = int(v_dim/2 * self.padding)
u_mag_pad = np.pad(u_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant', constant_values=0)
v_mag_pad = np.pad(v_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant', constant_values=0)
# Fourier transform of the two components:
u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_pad), axes=0)
v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_pad), axes=0)
# Calculate the Fourier transform of the phase:
f_nyq = 0.5 / self.a # nyquist frequency
f_u = np.linspace(0, f_nyq, u_mag_fft.shape[1])
f_v = np.linspace(-f_nyq, f_nyq, u_mag_fft.shape[0], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
coeff = (1j*self.b_0*self.a) / (2*self.PHI_0)
phase_fft = coeff * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
phase_pad = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
phase = phase_pad[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
return PhaseMap(mag_data.a, phase)
def jac_dot(self, vector):
self.LOG.debug('Calling jac_dot')
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
(row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
``N**2`` elements to the `v`-component of the magnetization.
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
self.LOG.debug('Calling jac_dot')
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv))
magnitude_proj = self.jac_dot(vector).reshape((2, )+self.dim_uv)
mag_proj.magnitude[1:3, 0, ...] = magnitude_proj
# TODO: instead call common subroutine operating on u_mag, v_mag with __call__?
return self(mag_proj).phase_vec
def jac_T_dot(self, vector):
self.LOG.debug('Calling jac_T_dot')
raise NotImplementedError()
# TODO: Implement!
class PhaseMapperElectric(PhaseMapper):
'''Class representing a phase mapping strategy for the electrostatic contribution.
The :class:`~.PMElectric` class represents a phase mapping strategy for the electrostatic
contribution to the electron phase shift which results e.g. from the mean inner potential in
certain samples and which is sensitive to properties of the electron microscope. It directly
takes :class:`~.MagData` objects and returns :class:`~.PhaseMap` objects.
Attributes
----------
a : float
The grid spacing in nm.
projector : :class:`~.Projector`
Projector which should be used for the projection of the 3-dimensional magnetization
distribution.
v_0 : float, optional
The mean inner potential of the specimen in V. The default is 1.
v_acc : float, optional
The acceleration voltage of the electron microscope in V. The default is 30000.
threshold : float, optional
Threshold for the recognition of the specimen location. The default is 0, which means that
all voxels with non-zero magnetization will contribute. Should be above noise level.
'''
LOG = logging.getLogger(__name__+'.PMElectric')
H_BAR = 6.626E-34 # Planck constant in J*s
M_E = 9.109E-31 # electron mass in kg
Q_E = 1.602E-19 # electron charge in C
C = 2.998E8 # speed of light in m/s
def __init__(self, a, dim_uv, v_0=1, v_acc=30000, threshold=0):
self.LOG.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.v_0 = v_0
self.v_acc = v_acc
self.threshold = threshold
self.m = np.prod(self.dim_uv)
self.n = np.prod(self.dim_uv)
# Coefficient calculation:
H_BAR, M_E, Q_E, C = self.H_BAR, self.M_E, self.Q_E, self.C
lam = H_BAR / np.sqrt(2 * M_E * Q_E * v_acc * (1 + Q_E*v_acc / (2*M_E*C**2)))
C_e = 2*pi*Q_E/lam * (Q_E*v_acc + M_E*C**2) / (Q_E*v_acc * (Q_E*v_acc + 2*M_E*C**2))
self.coeff = v_0 * C_e * a * 1E-9
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, v_0=%r, v_acc=%r, threshold=%r)' % \
(self.__class__, self.a, self.dim_uv, self.v_0, self.v_acc, self.threshold)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'PhaseMapperElectric(a=%s, dim_uv=%s, v_0=%s, v_acc=%s, threshold=%s)' % \
(self.a, self.dim_uv, self.v_0, self.v_acc, self.threshold)
def __call__(self, mag_data):
self.LOG.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert mag_data.dim[1:3] == self.dim_uv, 'Dimensions do not match!'
return self.coeff * mag_data.get_mask(self.threshold)[0, ...].reshape(self.dim_uv)
# Calculate mask:
mask = mag_data.get_mask(self.threshold)
# Project and calculate phase:
# TODO: check if projector manages scalar multiplication (problem with MagData) and fix!
projection = self.projector(mask.reshape(-1)).reshape(self.projector.dim_uv)
phase = self.coeff * projection
return PhaseMap(mag_data.a, phase)
def jac_dot(self, vector):
self.LOG.debug('Calling jac_dot')
raise NotImplementedError()
# TODO: Implement?
def jac_T_dot(self, vector):
self.LOG.debug('Calling jac_T_dot')
raise NotImplementedError()
# TODO: Implement?
def pm(mag_data, axis='z', dim_uv=None, b_0=1):
'''Convenience function for fast phase mapping.
Parameters
----------
mag_data : :class:`~.MagData`
A :class:`~.MagData` object, from which the projected phase map should be calculated.
axis: {'z', 'y', 'x'}, optional
Axis along which the :class:`.~SimpleProjector` projects the magnetic distribution.
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.
Returns
-------
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
'''
projector = SimpleProjector(mag_data.dim, axis=axis, dim_uv=dim_uv)
phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
return phasemapper(projector(mag_data))
pyramid/projector.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""This module provides the abstract base class :class:`~.Projector` and concrete subclasses for
projections of vector and scalar fields."""
import numpy as np
from numpy import pi
import abc
import itertools
from scipy.sparse import coo_matrix, csr_matrix
from pyramid.magdata import MagData
import logging
class Projector(object):
'''Abstract base class representing a projection function.
The :class:`~.Projector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid. :class:`~.Projector` is an abstract base
class and provides a unified interface which should be subclassed with a custom
:func:`__init__` function, which should call the parent :func:`__init__` method. Concrete
subclasses can be called as a function and take a `vector` as argument which contains the
3-dimensional field. The output is the projected field, given as a `vector`. Depending on the
length of the input and the given dimensions `dim` at construction time, vector or scalar
projection is choosen intelligently.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
dim_uv : tuple (N=2)
Dimensions (v, u) of the projected grid.
size_3d : int
Number of voxels of the 3-dimensional grid.
size_2d : int
Number of pixels of the 2-dimensional projected grid.
weight : :class:`~scipy.sparse.csr_matrix` (N=2)
The weight matrix containing the weighting coefficients for the 3D to 2D mapping.
coeff : list (N=2)
List containing the six weighting coefficients describing the influence of the 3 components
of a 3-dimensional vector field on the 2 projected components.
n: int
Size of the image space.
m: int
Size of the input space.
'''
__metaclass__ = abc.ABCMeta
LOG = logging.getLogger(__name__+'.Projector')
@abc.abstractmethod
def __init__(self, dim, dim_uv, weight, coeff):
self.LOG.debug('Calling __init__')
self.dim = dim
self.dim_uv = dim_uv
self.weight = weight
self.coeff = coeff
self.size_2d, self.size_3d = weight.shape
self.m = 3 * np.prod(dim)
self.n = 2 * np.prod(dim_uv)
self.LOG.debug('Created '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(dim=%r, dim_uv=%r, weight=%r, coeff=%r)' % \
(self.__class__, self.dim, self.dim_uv, self.weight, self.coeff)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'Projector(dim=%s, dim_uv=%s, coeff=%s)' % (self.dim, self.dim_uv, self.coeff)
def __call__(self, mag_data):
self.LOG.debug('Calling as function')
mag_proj = MagData(mag_data.a, np.zeros((3, 1)+self.dim_uv))
magnitude_proj = self.jac_dot(mag_data.mag_vec).reshape((2, )+self.dim_uv)
mag_proj.magnitude[0:2, 0, ...] = magnitude_proj
return mag_proj
def _vector_field_projection(self, vector):
self.LOG.debug('Calling _vector_field_projection')
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(2*size_2d)
# Go over all possible component projections (z, y, x) to (u, v):
if self.coeff[0][0] != 0: # x to u
result[:size_2d] += self.coeff[0][0] * self.weight.dot(vector[:size_3d])
if self.coeff[0][1] != 0: # y to u
result[:size_2d] += self.coeff[0][1] * self.weight.dot(vector[size_3d:2*size_3d])
if self.coeff[0][2] != 0: # z to u
result[:size_2d] += self.coeff[0][2] * self.weight.dot(vector[2*size_3d:])
if self.coeff[1][0] != 0: # x to v
result[size_2d:] += self.coeff[1][0] * self.weight.dot(vector[:size_3d])
if self.coeff[1][1] != 0: # y to v
result[size_2d:] += self.coeff[1][1] * self.weight.dot(vector[size_3d:2*size_3d])
if self.coeff[1][2] != 0: # z to v
result[size_2d:] += self.coeff[1][2] * self.weight.dot(vector[2*size_3d:])
return result
def _vector_field_projection_T(self, vector):
self.LOG.debug('Calling _vector_field_projection_T')
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(3*size_3d)
# Go over all possible component projections (u, v) to (z, y, x):
if self.coeff[0][0] != 0: # u to x
result[:size_3d] += self.coeff[0][0] * self.weight.T.dot(vector[:size_2d])
if self.coeff[0][1] != 0: # u to y
result[size_3d:2*size_3d] += self.coeff[0][1] * self.weight.T.dot(vector[:size_2d])
if self.coeff[0][2] != 0: # u to z
result[2*size_3d:] += self.coeff[0][2] * self.weight.T.dot(vector[:size_2d])
if self.coeff[1][0] != 0: # v to x
result[:size_3d] += self.coeff[1][0] * self.weight.T.dot(vector[size_2d:])
if self.coeff[1][1] != 0: # v to y
result[size_3d:2*size_3d] += self.coeff[1][1] * self.weight.T.dot(vector[size_2d:])
if self.coeff[1][2] != 0: # v to z
result[2*size_3d:] += self.coeff[1][2] * self.weight.T.dot(vector[size_2d:])
return result
def _scalar_field_projection(self, vector):
self.LOG.debug('Calling _scalar_field_projection')
return np.array(self.weight.dot(vector))
def _scalar_field_projection_T(self, vector):
self.LOG.debug('Calling _scalar_field_projection_T')
return np.array(self.weight.T.dot(vector))
def jac_dot(self, vector):
'''Multiply a `vector` with the jacobi matrix of this :class:`~.Projector` object.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vector containing the field which should be projected. Must have the same or 3 times
the size of `size_3d` of the projector for scalar and vector projection, respectively.
Returns
-------
proj_vector : :class:`~numpy.ndarray` (N=1)
Vector containing the projected field of the 2-dimensional grid. The length is
always`size_2d`.
'''
self.LOG.debug('Calling jac_dot')
if len(vector) == 3*self.size_3d: # mode == 'vector'
self.LOG.debug('mode == vector')
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
self.LOG.debug('mode == scalar')
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for '
'vector- or scalar-field-projection!')
def jac_T_dot(self, vector):
'''Multiply a `vector` with the transp. jacobi matrix of this :class:`~.Projector` object.
Parameters
----------
vector : :class:`~numpy.ndarray` (N=1)
Vector containing the field which should be projected. Must have the same or 2 times
the size of `size_2d` of the projector for scalar and vector projection, respectively.
Returns
-------
proj_vector : :class:`~numpy.ndarray` (N=1)
Vector containing the multiplication of the input with the transposed jacobi matrix
of the :class:`~.Projector` object.
'''
self.LOG.debug('Calling jac_T_dot')
if len(vector) == 2*self.size_2d: # mode == 'vector'
self.LOG.debug('mode == vector')
return self._vector_field_projection_T(vector)
elif len(vector) == self.size_2d: # mode == 'scalar'
self.LOG.debug('mode == scalar')
return self._scalar_field_projection_T(vector)
else:
raise AssertionError('Vector size has to be suited either for '
'vector- or scalar-field-projection!')
@abc.abstractmethod
def get_info(self):
'''Get specific information about the projector as a string.
Parameters
----------
None
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
'''
raise NotImplementedError()
class XTiltProjector(Projector):
'''Class representing a projection function with a tilt around the x-axis.
The :class:`~.XTiltProjector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
:class:`~.Projector`.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
tilt : float
Angle in `rad` describing the tilt of the beam direction relative to the x-axis.
dim_uv : tuple (N=2), optional
Dimensions (v, u) of the projection. If not set defaults to the (y, x)-dimensions.
'''
LOG = logging.getLogger(__name__+'.XTiltProjector')
def __init__(self, dim, tilt, dim_uv=None):
def get_position(p, m, b, size):
self.LOG.debug('Calling get_position')
y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
def get_impact(pos, r, size):
self.LOG.debug('Calling get_impact')
return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
if 0 <= x < size]
def get_weight(delta, rho): # use circles to represent the voxels
self.LOG.debug('Calling get_weight')
lo, up = delta-rho, delta+rho
# Upper boundary:
if up >= 1:
w_up = 0.5
else:
w_up = (up*np.sqrt(1-up**2) + np.arctan(up/np.sqrt(1-up**2))) / pi
# Lower boundary:
if lo <= -1:
w_lo = -0.5
else:
w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
return w_up - w_lo
self.LOG.debug('Calling __init__')
self.tilt = tilt
# Set starting variables:
# length along projection (proj, z), perpendicular (perp, y) and rotation (rot, x) axis:
dim_proj, dim_perp, dim_rot = dim
if dim_uv is None:
dim_uv = (max(dim_perp, dim_proj), dim_rot) # x-y-plane
dim_v, dim_u = dim_uv # y, x
assert dim_v >= dim_perp and dim_u >= dim_rot, 'Projected dimensions are too small!'
size_2d = np.prod(dim_uv)
size_3d = np.prod(dim)
# Creating coordinate list of all voxels:
voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-y-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj/2., dim_perp/2.)
m = np.where(tilt <= pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
b = center[0] - m * center[1]
positions = get_position(voxels, m, b, dim_v)
# Calculate weight-matrix:
r = 1/np.sqrt(np.pi) # radius of the voxel circle
rho = 0.5 / r
row = []
col = []
data = []
# one slice:
for i, voxel in enumerate(voxels):
impacts = get_impact(positions[i], r, dim_v) # impact along projected y-axis
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
col.append(voxel[0]*dim_rot*dim_perp + voxel[1]*dim_rot) # 0: z, 1: y
row.append(impact*dim_u+int((dim_u-dim_rot)/2))
data.append(get_weight(delta, rho))
# All other slices (along x):
columns = col
rows = row
for s in np.arange(1, dim_rot):
columns = np.hstack((np.array(columns), np.array(col)+s))
rows = np.hstack((np.array(rows), np.array(row)+s))
# Calculate weight matrix and coefficients for jacobi matrix:
weight = csr_matrix(coo_matrix((np.tile(data, dim_rot), (rows, columns)),
shape=(size_2d, size_3d)))
coeff = [[1, 0, 0], [0, np.cos(tilt), np.sin(tilt)]]
super(XTiltProjector, self).__init__(dim, dim_uv, weight, coeff)
self.LOG.debug('Created '+str(self))
def get_info(self):
'''Get specific information about the projector as a string.
Parameters
----------
None
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
'''
return 'x-tilt: $\phi = {:3.2f} \pi$'.format(self.tilt/pi)
class YTiltProjector(Projector):
'''Class representing a projection function with a tilt around the y-axis.
The :class:`~.YTiltProjector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
:class:`~.Projector`.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
tilt : float
Angle in `rad` describing the tilt of the beam direction relative to the y-axis.
dim_uv : tuple (N=2), optional
Dimensions (v, u) of the projection. If not set defaults to the (y, x)-dimensions.
'''
LOG = logging.getLogger(__name__+'.YTiltProjector')
def __init__(self, dim, tilt, dim_uv=None):
def get_position(p, m, b, size):
y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
def get_impact(pos, r, size):
return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
if 0 <= x < size]
def get_weight(delta, rho): # use circles to represent the voxels
lo, up = delta-rho, delta+rho
# Upper boundary:
if up >= 1:
w_up = 0.5
else:
w_up = (up*np.sqrt(1-up**2) + np.arctan(up/np.sqrt(1-up**2))) / pi
# Lower boundary:
if lo <= -1:
w_lo = -0.5
else:
w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
return w_up - w_lo
self.LOG.debug('Calling __init__')
self.tilt = tilt
# Set starting variables:
# length along projection (proj, z), rotation (rot, y) and perpendicular (perp, x) axis:
dim_proj, dim_rot, dim_perp = dim
if dim_uv is None:
dim_uv = (dim_rot, max(dim_perp, dim_proj)) # x-y-plane
dim_v, dim_u = dim_uv # y, x
assert dim_v >= dim_rot and dim_u >= dim_perp, 'Projected dimensions are too small!'
size_2d = np.prod(dim_uv)
size_3d = np.prod(dim)
# Creating coordinate list of all voxels:
voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-x-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj/2., dim_perp/2.)
m = np.where(tilt <= pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
b = center[0] - m * center[1]
positions = get_position(voxels, m, b, dim_u)
# Calculate weight-matrix:
r = 1/np.sqrt(np.pi) # radius of the voxel circle
rho = 0.5 / r
row = []
col = []
data = []
# one slice:
for i, voxel in enumerate(voxels):
impacts = get_impact(positions[i], r, dim_u) # impact along projected x-axis
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
col.append(voxel[0]*dim_perp*dim_rot + voxel[1]) # 0: z, 1: x
row.append(impact+int((dim_v-dim_rot)/2)*dim_u)
data.append(get_weight(delta, rho))
# All other slices (along y):
columns = col
rows = row
for s in np.arange(1, dim_rot):
columns = np.hstack((np.array(columns), np.array(col)+s*dim_perp))
rows = np.hstack((np.array(rows), np.array(row)+s*dim_u))
# Calculate weight matrix and coefficients for jacobi matrix:
weight = csr_matrix(coo_matrix((np.tile(data, dim_rot), (rows, columns)),
shape=(size_2d, size_3d)))
coeff = [[np.cos(tilt), 0, np.sin(tilt)], [0, 1, 0]]
super(YTiltProjector, self).__init__(dim, dim_uv, weight, coeff)
self.LOG.debug('Created '+str(self))
def get_info(self):
'''Get specific information about the projector as a string.
Parameters
----------
None
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
'''
return 'y-tilt: $\phi = {:3.2f} \pi$'.format(self.tilt/pi)
class SimpleProjector(Projector):
'''Class representing a projection function along one of the major axes.
The :class:`~.SimpleProjector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
:class:`~.Projector`.
Attributes
----------
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
axis : {'z', 'y', 'x'}, optional
Main axis along which the magnetic distribution is projected (given as a string). Defaults
to the z-axis.
dim_uv : tuple (N=2), optional
Dimensions (v, u) of the projection. If not set it uses the 3D default dimensions.
'''
LOG = logging.getLogger(__name__+'.SimpleProjector')
AXIS_DICT = {'z': (0, 1, 2), 'y': (1, 0, 2), 'x': (2, 1, 0)} # (0:z, 1:y, 2:x) -> (proj, v, u)
def __init__(self, dim, axis='z', dim_uv=None):
self.LOG.debug('Calling __init__')
assert axis in {'z', 'y', 'x'}, 'Projection axis has to be x, y or z (given as a string)!'
proj, v, u = self.AXIS_DICT[axis]
dim_proj, dim_v, dim_u = dim[proj], dim[v], dim[u]
dim_z, dim_y, dim_x = dim
size_2d = dim_u * dim_v
size_3d = np.prod(dim)
data = np.repeat(1, size_3d) # size_3d ones in the matrix (each voxel is projected)
indptr = np.arange(0, size_3d+1, dim_proj) # each row has dim_proj ones
if axis == 'z':
self.LOG.debug('Projecting along the z-axis')
coeff = [[1, 0, 0], [0, 1, 0]]
indices = np.array([np.arange(row, size_3d, size_2d)
for row in range(size_2d)]).reshape(-1)
elif axis == 'y':
self.LOG.debug('Projection along the y-axis')
coeff = [[1, 0, 0], [0, 0, 1]]
indices = np.array([np.arange(row%dim_x, dim_x*dim_y, dim_x)+int(row/dim_x)*dim_x*dim_y
for row in range(size_2d)]).reshape(-1)
elif axis == 'x':
self.LOG.debug('Projection along the x-axis')
coeff = [[0, 0, 1], [0, 1, 0]] # Caution, coordinate switch: u, v --> z, y (not y, z!)
# indices = np.array([np.arange(dim_proj) + row*dim_proj
# for row in range(size_2d)]).reshape(-1) # this is u, v --> y, z
indices = np.array([np.arange(dim_x) + (row%dim_z)*dim_x*dim_y + int(row/dim_z)*dim_x
for row in range(size_2d)]).reshape(-1)
if dim_uv is not None:
indptr = indptr.tolist() # convert to use insert() and append()
d_v, d_u = int((dim_uv[0]-dim_v)/2), int((dim_uv[1]-dim_u)/2) # padding in u and v
indptr[-1:-1] = [indptr[-1]] * d_v*dim_uv[1] # append empty rows at the end
for i in np.arange(dim_v, 0, -1): # all slices in between
u, l = i*dim_u, (i-1)*dim_u+1 # upper / lower slice end
indptr[u:u] = [indptr[u]] * d_u # end of the slice
indptr[l:l] = [indptr[l]] * d_u # start of the slice
indptr[0:0] = [0] * d_v*dim_uv[1] # insert empty rows at the beginning
size_2d = np.prod(dim_uv) # increase size_2d
# Make sure dim_uv is defined (used for the assertion)
if dim_uv is None:
dim_uv = dim_v, dim_u
assert dim_uv[0] >= dim_v and dim_uv[1] >= dim_u, 'Projected dimensions are too small!'
# Create weight-matrix:
weight = csr_matrix((data, indices, indptr), shape=(size_2d, size_3d))
super(SimpleProjector, self).__init__(dim, dim_uv, weight, coeff)
self.LOG.debug('Created '+str(self))
def get_info(self):
'''Get specific information about the projector as a string.
Parameters
----------
None
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
'''
return 'projected along {}-axis'.format(self.axis)
pyramid/reconstruction.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Reconstruct magnetic distributions from given phasemaps.
This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
objects) by using several model based reconstruction algorithms which use the forward model
provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
the distribution.
"""
import numpy as np
from scipy.sparse.linalg import cg
from scipy.optimize import minimize, leastsq
from pyramid.kernel import Kernel
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PhaseMapperRDFC
from pyramid.forwardmodel import ForwardModel
from pyramid.costfunction import Costfunction, CFAdapterScipyCG
from pyramid.magdata import MagData
import logging
LOG = logging.getLogger(__name__)
class PrintIterator(object):
'''Iterator class which is responsible to give feedback during reconstruction iterations.
Parameters
----------
cost : :class:`~.Costfunction`
:class:`~.Costfunction` class for outputting the `cost` of the current magnetization
distribution. This should decrease per iteration if the algorithm converges and is only
printed for a `verbosity` of 2.
verbosity : {0, 1, 2}, optional
Parameter defining the verbosity of the output. `2` will show the current number of the
iteration and the cost of the current distribution. `1` will just show the iteration
number and `0` will prevent output all together.
Notes
-----
Normally this class should not be used by the user and is instantiated whithin the
:mod:`~.reconstruction` module itself.
'''
LOG = logging.getLogger(__name__+'.PrintIterator')
def __init__(self, cost, verbosity):
self.LOG.debug('Calling __init__')
self.cost = cost
self.verbosity = verbosity
assert verbosity in {0, 1, 2}, 'verbosity has to be set to 0, 1 or 2!'
self.iteration = 0
self.LOG.debug('Created '+str(self))
def __call__(self, xk):
self.LOG.debug('Calling __call__')
if self.verbosity == 0:
return
print 'iteration #', self.next(),
if self.verbosity > 1:
print 'cost =', self.cost(xk)
else:
print ''
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(cost=%r, verbosity=%r)' % (self.__class__, self.cost, self.verbosity)
def __str__(self):
self.LOG.debug('Calling __str__')
return 'PrintIterator(cost=%s, verbosity=%s)' % (self.cost, self.verbosity)
def next(self):
self.iteration += 1
return self.iteration
def optimize_sparse_cg(data, regularisator=None, maxiter=1000, verbosity=0):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
Parameters
----------
data : :class:`~.DataSet`
:class:`~.DataSet` object containing all phase maps in :class:`~.PhaseMap` objects and all
projection directions in :class:`~.Projector` objects. These provide the essential
information for the reconstruction.
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). Defaults to
an appropriate unity matrix if none is provided.
regularisator : :class:`~.Regularisator`, optional
Regularisator class that's responsible for the regularisation term. Defaults to zero
order Tikhonov if none is provided.
maxiter : int
Maximum number of iterations.
verbosity : {0, 1, 2}, optional
Parameter defining the verposity of the output. `2` will show the current number of the
iteration and the cost of the current distribution. `1` will just show the iteration
number and `0` is the default and will prevent output all together.
Returns
-------
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
'''
LOG.debug('Calling optimize_sparse_cg')
# Set up necessary objects:
y = data.phase_vec
fwd_model = ForwardModel(data)
cost = Costfunction(data, regularisator)
# Optimize:
A = CFAdapterScipyCG(cost)
b = fwd_model.jac_T_dot(None, y)
x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity), maxiter=maxiter)
# Create and return fitting MagData object:
mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
mag_opt.mag_vec = x_opt
return mag_opt
def optimize_cg(data, first_guess):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
Parameters
----------
data : :class:`~.DataSet`
:class:`~.DataSet` object containing all phase maps in :class:`~.PhaseMap` objects and all
projection directions in :class:`~.Projector` objects. These provide the essential
information for the reconstruction.
verbosity : {2, 1, 0}, optional
Parameter defining the verposity of the output. `2` is the default and will show the
current number of the iteration and the cost of the current distribution. `2` will just
show the iteration number and `0` will prevent output all together.
Returns
-------
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
'''
LOG.debug('Calling optimize_cg')
mag_0 = first_guess
x_0 = first_guess.mag_vec
y = data.phase_vec
kernel = Kernel(data.a, data.dim_uv)
fwd_model = ForwardModel(data.projectors, kernel, data.b_0)
cost = Costfunction(y, fwd_model)
# Optimize:
result = minimize(cost, x_0, method='Newton-CG', jac=cost.jac, hessp=cost.hess_dot,
options={'maxiter': 1000, 'disp': True})
# Create optimized MagData object:
x_opt = result.x
mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
mag_opt.mag_vec = x_opt
return mag_opt
def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
'''Reconstruct a magnetic distribution for a 2-D problem with known pixel locations.
Parameters
----------
phase_map : :class:`~pyramid.phasemap.PhaseMap`
A :class:`~pyramid.phasemap.PhaseMap` object, representing the phase from which to
reconstruct the magnetic distribution.
mask : :class:`~numpy.ndarray` (N=3)
A boolean matrix (or a matrix consisting of ones and zeros), representing the
positions of the magnetized voxels in 3 dimensions.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
lam : float, optional
The regularisation parameter. Defaults to 1E-4.
order : int {0, 1}, optional
order of the regularisation function. Default is 0 for a Tikhonov regularisation of order
zero. A first order regularisation, which uses the derivative is available with value 1.
Returns
-------
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
Notes
-----
Only works for a single phase_map, if the positions of the magnetized voxels are known and
for slice thickness of 1 (constraint for the `z`-dimension).
'''
# Read in parameters:
y_m = phase_map.phase_vec # Measured phase map as a vector
a = phase_map.a # Grid spacing
dim = mask.shape # Dimensions of the mag. distr.
count = mask.sum() # Number of pixels with magnetization
# Create empty MagData object for the reconstruction:
mag_data_rec = MagData(a, np.zeros((3,)+dim))
# Function that returns the phase map for a magnetic configuration x:
def F(x):
mag_data_rec.set_vector(x, mask)
proj = SimpleProjector(dim)
phase_map = PhaseMapperRDFC(Kernel(a, proj.dim_uv, b_0))(proj(mag_data_rec))
return phase_map.phase_vec
# Cost function of order zero which should be minimized:
def J_0(x_i):
y_i = F(x_i)
term1 = (y_i - y_m)
term2 = lam * x_i
return np.concatenate([term1, term2])
# First order cost function which should be minimized:
def J_1(x_i):
y_i = F(x_i)
term1 = (y_i - y_m)
mag_data = mag_data_rec.magnitude
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))
term2 = lam * np.concatenate(term2)
return np.concatenate([term1, term2])
J_DICT = [J_0, J_1] # list of cost-functions with different regularisations
# Reconstruct the magnetization components:
x_rec, _ = leastsq(J_DICT[order], np.zeros(3*count))
mag_data_rec.set_vector(x_rec, mask)
return mag_data_rec
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
<!DOCTYPE X3D PUBLIC 'ISO//Web3D//DTD X3D 3.2//EN' 'http://www.web3d.org/specifications/x3d-3.2.dtd'>
<X3D xmlns:cfn='_xml/_xsd/CommonFunctions' xmlns:cvs='CVSInfo' xmlns:fn='http://www.w3.org/2005/xpath-functions' xmlns:lfn='local_functions' xmlns:ss='urn:schemas-microsoft-com:office:spreadsheet' xmlns:v2kfn='_xml/_xsd/V2Kfunctions' xmlns:x3d='http://www.web3d.org/specifications/x3d-3.2.xsd' xmlns:xlink='http://www.w3.org/1999/xlink' xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:xsd='http://www.w3.org/2001/XMLSchema-instance' profile='Immersive' version='3.2' xsd:noNamespaceSchemaLocation='http://www.web3d.org/specifications/x3d-3.2.xsd'>
<Engine DEF='engine'>
<RenderJob DEF='render'>
<WindowGroup>
<Window position='10 50' size='933,700' fullScreen='false'/>
</WindowGroup>
</RenderJob>
</Engine>
<Scene>
<Background skyColor='0 0 0'/>
<ProximitySensor DEF='HereIAm' enabled='true' size='10000000 10000000 10000000'/>
<ProtoDeclare name='Spin_Proto'>
<ProtoInterface>
<field name='spin_position' accessType='inputOnly' type='SFVec3f'/>
<field name='spin_scale' accessType='inputOnly' type='SFVec3f'/>
<field name='spin_rotation' accessType='inputOnly' type='SFRotation'/>
<field name='spin_color' accessType='inputOnly' type='SFColor'/>
</ProtoInterface>
<ProtoBody>
<Group>
<Transform DEF='Spin'>
<IS>
<connect nodeField='scale' protoField='spin_scale'/>
<connect nodeField='translation' protoField='spin_position'/>
<connect nodeField='rotation' protoField='spin_rotation'/>
</IS>
<Transform translation='0 0.375 0'>
<Shape>
<Cone bottomRadius='0.25' bottom='true' solid='true' height='0.25' side='true'/>
<Appearance>
<Material DEF='ConeColor'>
<IS>
<connect nodeField='diffuseColor' protoField='spin_color'/>
</IS>
</Material>
</Appearance>
</Shape>
</Transform>
<Transform translation='0 -0.125 0'>
<Shape>
<Cylinder radius='0.1' height='0.75'/>
<Appearance>
<Material DEF='CylinderColor'>
<IS>
<connect nodeField='diffuseColor' protoField='spin_color'/>
</IS>
</Material>
</Appearance>
</Shape>
</Transform>
</Transform>
</Group>
</ProtoBody>
</ProtoDeclare>
<Transform DEF='ArrowsHUD'>
<Transform DEF='Arrows' translation='-4.5 -3 -10' scale='0.75 0.75 0.75'>
<Group DEF='ArrowGreen'>
<Shape>
<Cylinder DEF='ArrowCylinder' radius='0.025' top='false'/>
<Appearance DEF='Green'>
<Material diffuseColor='0 1 0'/>
</Appearance>
</Shape>
<Transform translation='0 1 0'>
<Shape>
<Cone DEF='ArrowCone' bottomRadius='0.05' height='0.1'/>
<Appearance USE='Green'/>
</Shape>
</Transform>
<Transform translation='0 1.1 0'>
<Billboard>
<Shape>
<Appearance DEF='LabelAppearance'>
<Material diffuseColor='1 1 0'/>
</Appearance>
<Text string='Y'>
<FontStyle DEF='LabelFont' size='0.25'/>
</Text>
</Shape>
</Billboard>
</Transform>
</Group>
<Group DEF='ArrowBlue'>
<Transform rotation='0 0 1 -1.57079'>
<Shape>
<Cylinder USE='ArrowCylinder'/>
<Appearance DEF='Blue'>
<Material diffuseColor='0 0 1'/>
</Appearance>
</Shape>
<Transform translation='0 1 0'>
<Shape>
<Cone USE='ArrowCone'/>
<Appearance USE='Blue'/>
</Shape>
</Transform>
<Transform rotation='0 0 1 1.57079' translation='0.072 1.1 0'>
<Billboard>
<Shape>
<Appearance USE='LabelAppearance'/>
<Text string='X'>
<FontStyle USE='LabelFont'/>
</Text>
</Shape>
</Billboard>
</Transform>
</Transform>
</Group>
<Group DEF='ArrowRed'>
<Transform rotation='1 0 0 1.57079'>
<Shape>
<Cylinder USE='ArrowCylinder'/>
<Appearance DEF='Red'>
<Material diffuseColor='1 0 0'/>
</Appearance>
</Shape>
<Transform translation='0 1 0'>
<Shape>
<Cone USE='ArrowCone'/>
<Appearance USE='Red'/>
</Shape>
</Transform>
<Transform rotation='1 0 0 -1.57079' translation='0 1.1 0.072'>
<Billboard>
<Shape>
<Appearance USE='LabelAppearance'/>
<Text string='Z'>
<FontStyle USE='LabelFont'/>
</Text>
</Shape>
</Billboard>
</Transform>
</Transform>
</Group>
</Transform>
</Transform>
<Script DEF='ChangeRot' >
<field accessType='inputOnly' type='SFRotation' name='orientation'/>
<field accessType='outputOnly' name='CoordRot' type='SFRotation'/>
<![CDATA[javascript:
function orientation(orientation){
CoordRot= new SFRotation(orientation[0], orientation[1], orientation[2], -orientation[3]);
}
]]>
<ROUTE fromNode='HereIAm' fromField='orientation_changed' toNode='ArrowsHUD' toField='rotation'/>
<ROUTE fromNode='HereIAm' fromField='position_changed' toNode='ArrowsHUD' toField='translation'/>
<ROUTE fromNode='HereIAm' fromField='orientation_changed' toNode='ChangeRot' toField='orientation'/>
<ROUTE fromNode='ChangeRot' fromField='CoordRot' toNode='Arrows' toField='set_rotation'/>
</Script>
</Scene>
</X3D>
scripts/collaborations/bielefeld.py deleted 100644 → 0
View file @ 04d9ed45
# -*- coding: utf-8 -*-
"""Created on Thu Oct 02 11:53:25 2014 @author: Jan"""
import os
import numpy as np
import matplotlib.pyplot as plt
import pyramid
from pyramid.magdata import MagData
from pyramid.phasemapper import pm
import logging
import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
PATH = '../../output/bielefeld/'
def load_from_llg(filename):
'''Custom for Bielefeld files'''
SCALE = 1 # .0E-9 / 1.0E-2 # From cm to nm
data = np.genfromtxt(filename, skip_header=2)
dim = tuple(np.genfromtxt(filename, dtype=int, skip_header=1, skip_footer=len(data[:, 0])))
a = (data[1, 0] - data[0, 0]) / SCALE
magnitude = data[:, 3:6].T.reshape((3,)+dim)
# import pdb; pdb.set_trace()
return MagData(a, magnitude)
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
## 2DCoFe:
mag_data_2DCoFe = MagData.load_from_llg(PATH+'magnetic_distribution_2DCoFe.txt')
mag_data_2DCoFe.quiver_plot(proj_axis='x')
plt.savefig(PATH+'magnetic_distribution_2DCoFe.png')
mag_data_2DCoFe.quiver_plot3d()
phase_map_2DCoFe = pm(mag_data_2DCoFe, axis='x')
phase_map_2DCoFe.display_combined(gain='auto')
# AH21FIN:
mag_data_AH21FIN = MagData.load_from_llg(PATH+'magnetic_distribution_AH21FIN.txt')
mag_data_AH21FIN.quiver_plot(proj_axis='x')
plt.savefig(PATH+'magnetic_distribution_AH21FIN.png')
mag_data_AH21FIN.quiver_plot3d()
phase_map_AH21FIN = pm(mag_data_AH21FIN, axis='x')
phase_map_AH21FIN.display_combined(gain='auto')
# AH41FIN:
mag_data_AH41FIN = MagData.load_from_llg(PATH+'magnetic_distribution_AH41FIN.txt')
mag_data_AH41FIN.quiver_plot(proj_axis='x')
plt.savefig(PATH+'magnetic_distribution_AH41FIN.png')
mag_data_AH41FIN.quiver_plot3d()
phase_map_AH41FIN = pm(mag_data_AH41FIN, axis='x')
phase_map_AH41FIN.display_combined(gain='auto')
# Chain:
mag_data_chain = load_from_llg(PATH+'magnetic_distribution_Chain.txt')
mag_data_chain.quiver_plot(proj_axis='x')
plt.savefig(PATH+'magnetic_distribution_Chain.png')
mag_data_chain.quiver_plot3d()
phase_map_chain = pm(mag_data_chain, axis='x')
phase_map_chain.display_combined(gain='auto')
# Cylinder:
mag_data_cyl = load_from_llg(PATH+'magnetic_distribution_Cylinder.txt')
mag_data_cyl.quiver_plot(proj_axis='z')
plt.savefig(PATH+'magnetic_distribution_Cylinder.png')
mag_data_cyl.quiver_plot3d()
phase_map_cyl = pm(mag_data_cyl, axis='z')
phase_map_cyl.display_combined(gain='auto')
# Ring:
mag_data_ring = load_from_llg(PATH+'magnetic_distribution_ring.txt')
mag_data_ring.quiver_plot(proj_axis='x')
plt.savefig(PATH+'magnetic_distribution_ring.png')
mag_data_ring.quiver_plot3d()
phase_map_ring = pm(mag_data_ring, axis='x')
phase_map_ring.display_combined(gain='auto')