Newer
Older
# -*- 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
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.
dim : tuple (N=3)
Dimensions (z, y, x) of the magnetization distribution.
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.
__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.weight = weight
self.coeff = coeff
self.size_2d, self.size_3d = weight.shape
self.n = 3 * np.prod(dim)
self.m = 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)
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 '
@abc.abstractmethod
def get_info(self):
'''Get specific information about the projector as a string.
Parameters
----------
Returns
-------
info : string
Information about the projector as a string, e.g. for the use in plot titles.
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):
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
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)),
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__')
# 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)
# 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)),
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)