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
Showing
with 0 additions and 2710 deletions
# -*- 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 magnetizatin data.
projector
Create projections of a given magnetization distribution.
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.
holoimage
Create holographic contour maps from a given phase map.
reconstructor
Reconstruct magnetic distributions from given phasemaps.
Subpackages
-----------
numcore
Provides fast numerical functions for core routines.
"""
# -*- 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
PHI_0 = -2067.83 # magnetic flux in T*nm²
def phase_mag_slab(dim, res, 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)`.
res : float
The resolution of the grid (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.
'''
# 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 = res * (center[1] + 0.5) # y0, x0 define the center of a pixel,
x0 = res * (center[2] + 0.5) # hence: (cellindex + 0.5) * resolution
Lz, Ly, Lx = res * width[0], res * width[1], res * width[2]
coeff = b_0 / (4*PHI_0)
# Create grid:
x = np.linspace(res/2, x_dim*res-res/2, num=x_dim)
y = np.linspace(res/2, y_dim*res-res/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
def phase_mag_disc(dim, res, 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)`.
res : float
The resolution of the grid (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.
'''
# 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 = res * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
x0 = res * (center[2] + 0.5) # hence: cellindex + 0.5
Lz = res * height
R = res * radius
coeff = - pi * b_0 / (2*PHI_0)
# Create grid:
x = np.linspace(res/2, x_dim*res-res/2, num=x_dim)
y = np.linspace(res/2, y_dim*res-res/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
def phase_mag_sphere(dim, res, 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)`.
res : float
The resolution of the grid (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.
'''
# 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 = res * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
x0 = res * (center[2] + 0.5) # hence: cellindex + 0.5
R = res * radius
coeff = - 2./3. * pi * b_0 / PHI_0
# Create grid:
x = np.linspace(res / 2, x_dim * res - res / 2, num=x_dim)
y = np.linspace(res / 2, y_dim * res - res / 2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
def phase_mag_vortex(dim, res, 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)`.
res : float
The resolution of the grid (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.
'''
# 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 = res * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
x0 = res * (center[2] + 0.5) # hence: cellindex + 0.5
Lz = res * height
R = res * radius
coeff = pi * b_0 * Lz / PHI_0
# Create grid:
x = np.linspace(res/2, x_dim*res-res/2, num=x_dim)
y = np.linspace(res/2, y_dim*res-res/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
# -*- coding: utf-8 -*-
"""Create holographic contour maps from a given phase map.
This module converts phase maps into holographic contour map. This basically means 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:`~.holoimage` function to create these holographic contour maps. It is
possible to use these as input for the :func:`~.display` to plot them, or just pass the
:class:`~pyramid.phasemap.PhaseMap` object to the :func:`~.display_combined` function to plot
the phase map and the holographic contour map next to each other.
"""
import numpy as np
from numpy import pi
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import NullLocator
from PIL import Image
from pyramid.phasemap import PhaseMap
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)]}
HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
def holo_image(phase_map, density=1):
'''Create a holographic contour map from a :class:`~pyramid.phasemap.PhaseMap` object.
Parameters
----------
phase_map : :class:`~pyramid.phasemap.PhaseMap`
A :class:`~pyramid.phasemap.PhaseMap` object storing the phase information.
density : float, optional
The gain factor for determining the number of contour lines. The default is 1.
Returns
-------
holo_image : :class:`~PIL.image`
The resulting holographic contour map with color encoded gradient.
'''
assert isinstance(phase_map, PhaseMap), 'phase_map has to be a PhaseMap object!'
# Extract the phase (considering the unit):
phase = phase_map.phase
# Calculate the holography image intensity:
img_holo = (1 + np.cos(density * phase)) / 2
# Calculate the phase gradients, expressed by magnitude and angle:
phase_grad_y, phase_grad_x = np.gradient(phase, phase_map.res, phase_map.res)
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:
phase_magnitude = np.sin(phase_magnitude/phase_magnitude.max() * pi / 2)
# Color code the angle and create the holography image:
rgba = HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * phase_magnitude.T * rgba[:, :, :3].T).T.astype(np.uint8)
holo_image = Image.fromarray(rgb)
return holo_image
def make_color_wheel():
'''Display a color wheel to illustrate the color coding of the gradient direction.
Parameters
----------
None
Returns
-------
None
'''
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 = 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())
def display(holo_image, title='Holographic Contour Map', axis=None, interpolation='none'):
'''Display the color coded holography image.
Parameters
----------
holo_image : :class:`~PIL.image`
The resulting holographic contour map with color encoded gradient.
title : string, optional
The title of the plot. The default is 'Holographic Contour Map'.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method. No interpolation is used in the default case.
Returns
-------
None
'''
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1, 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)
plt.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('x-axis [px]', fontsize=15)
axis.set_ylabel('y-axis [px]', fontsize=15)
def display_combined(phase_map, density=1, title='Combined Plot', interpolation='none'):
'''Display a given phase map and the resulting color coded holography image in one plot.
Parameters
----------
phase_map : :class:`~pyramid.phasemap.PhaseMap`
A :class:`~pyramid.phasemap.PhaseMap` object storing the phase information.
density : float, optional
The gain factor for determining the number of contour lines. The default is 1.
title : string, optional
The title of the plot. The default is 'Combined Plot'.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
Returns
-------
None
'''
# 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')
display(holo_image(phase_map, density), axis=holo_axis, interpolation=interpolation)
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
phase_map.display(axis=phase_axis)
# -*- 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
class Shapes:
'''Class containing functions for generating magnetic shapes.
The :class:`~.Shapes` class is a collection of some mehtods 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.
'''
@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.
'''
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 zylindrical 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.
'''
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 == 'z' or axis == 'y' or axis == '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 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.
'''
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 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.
'''
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 == 'z' or axis == 'y' or axis == '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.
'''
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 homogeneous 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.
'''
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 z_mag, y_mag, x_mag
def create_mag_dist_vortex(mag_shape, center=None, 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 `(x, y)` or 3D `(x, y, z)`, where `z` is discarded.
Is set to the center of the field of view if not specified.
Vortex center has to be between two pixels.
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.
'''
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
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 len(center) == 3: # if a 3D-center is given, just take the x and y components
center = (center[1], center[2])
assert len(center) == 2, 'Vortex center has to be defined in 3D or 2D or not at all!'
x = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
y = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
xx, yy = np.meshgrid(x, y)
phi = np.arctan2(yy, xx) - pi/2
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
return z_mag, y_mag, x_mag
# -*- coding: utf-8 -*-
"""Class for the storage of magnetizatin data.
This module provides the :class:`~.MagData` class whose instances can be used to store
magnetization data for 3 components for a 3-dimensional grid. It is possible to load data from
NetCDF4 or LLG (.txt) files or to save the data in these formats. Also plotting methods are
provided. See :class:`~.MagData` for further information.
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from mayavi import mlab
import netCDF4
class MagData:
'''Class for storing magnetization data.
Represents 3-dimensional magnetic distributions with 3 components which are stored as a
tuple of size 3 in `magnitude`. The object can be created empty by omitting the `magnitude`
paramter. The `magnitude` can be added later by using the :func:`~.add_magnitude` function.
This is useful, if the `magnitude` is more complex and more than one magnetized object should
be represented by the :class:`~.MagData` object, which can be added one after another by the
:func:`~.add_magnitude` function. The dimensions `dim` of the grid will be set as soon as the
magnitude is specified. However, `res` has to be always specified at construction time.
Attributes
----------
res : float
The resolution of the grid (grid spacing in nm)
dim : tuple (N=3)
Dimensions of the grid.
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The `z`-, `y`- and `x`-component of the magnetization vector for every 3D-gridpoint
as a tuple.
'''
def __init__(self, res, magnitude=None):
'''Constructor for a :class:`~.MagData` object for storing magnetization data.
Parameters
----------
res : float
The resolution of the grid (grid spacing) in nm.
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3), optional
The `z`-, `y`- and `x`-component of the magnetization vector for every 3D-gridpoint
as a tuple. Is zero everywhere if not specified.
Returns
-------
mag_data : :class:`~.MagData`
The 3D magnetic distribution as a :class:`~.MagData` object.
'''
if magnitude is not None:
dim = np.shape(magnitude[0])
assert len(dim) == 3, 'Magnitude has to be defined for a 3-dimensional grid!'
assert np.shape(magnitude[1]) == np.shape(magnitude[2]) == dim, \
'Dimensions of the magnitude components do not match!'
self.magnitude = magnitude
self.dim = dim
else:
self.magnitude = None
self.dim = None
self.res = res
def add_magnitude(self, magnitude):
'''Add a given magnitude to the magnitude of the :class:`~.MagData`.
Parameters
----------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The `z`-, `y`- and `x`-component of the magnetization vector for every 3D-gridpoint
as a tuple. If the :class:`~.MagData` object already has a `magnitude`, the added
one has to match the dimensions `dim`, otherwise the added magnitude sets `dim`.
Returns
-------
None
'''
if self.magnitude is not None: # Add magnitude to existing one
assert np.shape(magnitude) == (3,) + self.dim, \
'Added magnitude has to have the same dimensions!'
z_mag, y_mag, x_mag = self.magnitude
z_new, y_new, x_new = magnitude
z_mag += z_new
y_mag += y_new
x_mag += x_new
self.magnitude = (z_mag, y_mag, x_mag)
else: # Magnitude is set for the first time and the dimensions are calculated
dim = np.shape(magnitude[0])
assert len(dim) == 3, 'Magnitude has to be defined for a 3-dimensional grid!'
assert np.shape(magnitude[1]) == np.shape(magnitude[2]) == dim, \
'Dimensions of the magnitude components do not match!'
self.magnitude = magnitude
self.dim = dim
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.concatenate([self.magnitude[2][mask],
self.magnitude[1][mask],
self.magnitude[0][mask]])
def set_vector(self, mask, vector):
'''Set the magnetic components of the masked pixels to the values specified by `vector`.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean)
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
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
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 resolution accordingly.
Only possible, if each axis length is a power of 2!
'''
assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
assert self.dim[0] % 2**n == 0 and self.dim[1] % 2**n == 0 and self.dim[2] % 2**n == 0, \
'For downscaling, every dimension must be a multiple of 2!'
for t in range(n):
# Create coarser grid for the magnetization:
z_mag = self.magnitude[0].reshape(self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2)
y_mag = self.magnitude[1].reshape(self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2)
x_mag = self.magnitude[2].reshape(self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2)
self.dim = (self.dim[0]/2, self.dim[1]/2, self.dim[2]/2)
self.res = self.res * 2
self.magnitude = (z_mag.mean(axis=5).mean(axis=3).mean(axis=1),
y_mag.mean(axis=5).mean(axis=3).mean(axis=1),
x_mag.mean(axis=5).mean(axis=3).mean(axis=1))
@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.
'''
SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
data = np.genfromtxt(filename, skip_header=2)
x_dim, y_dim, z_dim = np.genfromtxt(filename, dtype=int, skip_header=1,
skip_footer=len(data[:, 0]))
res = (data[1, 0] - data[0, 0]) / SCALE
# Reshape in Python and Igor is different, Python fills rows first, Igor columns:
x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim) for i in range(3, 6)]
return MagData(res, (z_mag, y_mag, x_mag))
def save_to_llg(self, filename='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 'magdata_output.txt'.
Returns
-------
None
'''
dim = self.dim
res = self.res * 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[res/2:(dim[0]*res-res/2):dim[0]*1j,
res/2:(dim[1]*res-res/2):dim[1]*1j,
res/2:(dim[2]*res-res/2):dim[2]*1j]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
x_mag = np.reshape(self.magnitude[2], (-1))
y_mag = np.reshape(self.magnitude[1], (-1))
z_mag = np.reshape(self.magnitude[0], (-1))
# Save data to file:
data = np.array([xx, yy, zz, x_mag, y_mag, z_mag]).T
with open(filename, 'w') as mag_file:
mag_file.write('LLGFileCreator: %s\n' % filename.replace('.txt', ''))
mag_file.write(' %d %d %d\n' % (dim[2], dim[1], dim[0]))
mag_file.writelines('\n'.join(' '.join('{:7.6e}'.format(cell)
for cell in row) for row in data))
@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.
'''
mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
res = mag_file.res
z_mag = mag_file.variables['z_mag'][:]
y_mag = mag_file.variables['y_mag'][:]
x_mag = mag_file.variables['x_mag'][:]
mag_file.close()
return MagData(res, (z_mag, y_mag, x_mag))
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 'magdata_output.nc'.
Returns
-------
None
'''
mag_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
mag_file.res = self.res
mag_file.createDimension('z_dim', self.dim[0])
mag_file.createDimension('y_dim', self.dim[1])
mag_file.createDimension('x_dim', self.dim[2])
z_mag = mag_file.createVariable('z_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
y_mag = mag_file.createVariable('y_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
x_mag = mag_file.createVariable('x_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
z_mag[:] = self.magnitude[0]
y_mag[:] = self.magnitude[1]
x_mag[:] = self.magnitude[2]
print mag_file
mag_file.close()
def quiver_plot(self, title='Magnetic Distribution', filename=None, axis=None,
proj_axis='z', ax_slice=None):
'''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.
filename : string, optional
The filename, specifying the location where the image is saved. If not specified,
the image is shown instead.
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.
Returns
-------
None
'''
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
if ax_slice is None:
ax_slice = int(self.dim[0]/2)
mag_slice_u = self.magnitude[2][ax_slice, ...]
mag_slice_v = self.magnitude[1][ax_slice, ...]
u_label = 'x [px]'
v_label = 'y [px]'
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
if ax_slice is None:
ax_slice = int(self.dim[1]/2)
mag_slice_u = self.magnitude[2][:, ax_slice, :]
mag_slice_v = self.magnitude[0][:, ax_slice, :]
u_label = 'x [px]'
v_label = 'z [px]'
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
if ax_slice is None:
ax_slice = int(self.dim[2]/2)
mag_slice_u = self.magnitude[1][..., ax_slice]
mag_slice_v = self.magnitude[0][..., ax_slice]
u_label = 'y [px]'
v_label = 'z [px]'
# 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, aspect='equal')
axis.quiver(mag_slice_u, mag_slice_v, pivot='middle', angles='xy', scale_units='xy',
scale=1, headwidth=6, headlength=7)
axis.set_xlim(-1, np.shape(mag_slice_u)[1])
axis.set_ylim(-1, np.shape(mag_slice_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=8, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=8, integer=True))
plt.show()
def quiver_plot3d(self):
'''Plot the magnetization as 3D-vectors in a quiverplot.
Parameters
----------
None
Returns
-------
None
'''
res = self.res
dim = self.dim
# Create points and vector components as lists:
zz, yy, xx = np.mgrid[res/2:(dim[0]*res-res/2):dim[0]*1j,
res/2:(dim[1]*res-res/2):dim[1]*1j,
res/2:(dim[2]*res-res/2):dim[2]*1j]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
x_mag = np.reshape(self.magnitude[2], (-1))
y_mag = np.reshape(self.magnitude[1], (-1))
z_mag = np.reshape(self.magnitude[0], (-1))
# Plot them as vectors:
mlab.figure()
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow', scale_factor=10.0)
mlab.outline(plot)
mlab.axes(plot)
mlab.colorbar()
"""Package for numerical core routines which enable fast computations.
Modules
-------
phase_mag_real
Provides numerical core routines for the real space approach in
:func:`~pyramid.phasemapper.phase_mag_real` of module :mod:`~pyramid.phasemapper`.
Notes
-----
Packages are written in `pyx`-format for the use with :mod:`~cython`.
"""
from phase_mag_real import *
# -*- coding: utf-8 -*-
"""Numerical core routines for the phase calculation using the real space approach.
Provides a helper function to speed up :func:`~pyramid.phasemapper.phase_mag_real` of module
:mod:`~pyramid.phasemapper`, 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 phase_mag_real_core(
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]
# -*- coding: utf-8 -*-
"""Class for the storage of phase data.
This module provides the :class:`~.PhaseMap` class whose instances can be used to store
phase data for a 2-dimensional grid. It is possible to load data from NetCDF4 or LLG (.txt) files
or to save the data in these formats. Also plotting methods are provided. See :class:`~.PhaseMap`
for further information.
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import netCDF4
class PhaseMap:
'''Class for storing phase map data.
Represents 2-dimensional phase maps. the phase information itself is stored in `phase`.
The dimensions `dim` of the grid with resolution `res` will be calculated at construction
time, but `res` has to be specified.
Attributes
----------
res : float
The resolution of the grid (grid spacing in nm)
dim : tuple (N=2)
Dimensions of the grid.
phase : :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
unit : {'rad', 'mrad', 'µrad'}, 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`.
'''
UNITDICT = {'rad': 1E0,
'mrad': 1E3,
'µrad': 1E6}
def __init__(self, res, phase, unit='rad'):
'''Constructor for a :class:`~.PhaseMap` object for storing phase data.
Parameters
----------
res : float
The resolution of the grid (grid spacing) in nm.
phase : :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
unit : {'rad', 'mrad', 'µrad'}, 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`.
Returns
-------
phase_map : :class:`~.PhaseMap`
The 2D phase shift as a :class:`~.PhaseMap` object.
'''
dim = np.shape(phase)
assert len(dim) == 2, 'Phasemap has to be 2-dimensional!'
self.res = res
self.dim = dim
self.unit = unit
self.phase = phase
def set_unit(self, unit):
'''Set the unit for the phase map.
Parameters
----------
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`.
Returns
-------
None
'''
assert unit in ['rad', 'mrad']
self.unit = unit
@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.
'''
with open(filename, 'r') as phase_file:
phase_file.readline() # Headerline is not used
res = float(phase_file.readline()[13:-4])
phase = np.loadtxt(filename, delimiter='\t', skiprows=2)
return PhaseMap(res, phase)
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 'phasemap_output.txt'.
Returns
-------
None
'''
with open(filename, 'w') as phase_file:
phase_file.write('{}\n'.format(filename.replace('.txt', '')))
phase_file.write('resolution = {} nm\n'.format(self.res))
np.savetxt(phase_file, self.phase, fmt='%7.6e', delimiter='\t')
@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.
'''
phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
res = phase_file.res
phase = phase_file.variables['phase'][:]
phase_file.close()
return PhaseMap(res, 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 'phasemap_output.nc'.
Returns
-------
None
'''
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
phase_file.res = self.res
phase_file.createDimension('v_dim', self.dim[0])
phase_file.createDimension('u_dim', self.dim[1])
phase = phase_file.createVariable('phase', 'f', ('v_dim', 'u_dim'))
phase[:] = self.phase
print phase_file
phase_file.close()
def display(self, title='Phase Map', cmap='RdBu', limit=None, norm=None, axis=None):
'''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 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.
Returns
-------
None
'''
# 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, 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, np.shape(phase)[1])
axis.set_ylim(0, np.shape(phase)[0])
ticks = (axis.get_xticks()*self.res).astype(int)
axis.set_xticklabels(ticks)
ticks = (axis.get_yticks()*self.res).astype(int)
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_yticklabels(ticks)
axis.set_title(title, fontsize=18)
axis.set_xlabel('x [nm]', fontsize=15)
axis.set_ylabel('y [nm]', fontsize=15)
# Add colorbar:
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('phase shift [{}]'.format(self.unit), fontsize=15)
# Show plot:
plt.show()
def display3d(self, title='Phase Map', cmap='RdBu'):
'''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'.
Returns
-------
None
'''
# 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[1])
v = range(self.dim[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('x-axis [px]')
axis.set_ylabel('y-axis [px]')
axis.set_zlabel('phase shift [{}]'.format(self.unit))
# Show Plot:
plt.show()
This diff is collapsed.
# -*- coding: utf-8 -*-
"""Create projections of a given magnetization distribution.
This module creates 2-dimensional projections from 3-dimensional magnetic distributions, which
are stored in :class:`~pyramid.magdata.MagData` objects. Either simple projections along the
major axes are possible (:func:`~.simple_axis_projection`), or projections along arbitrary
directions with the use of transfer functions (work in progress).
"""
import time
import numpy as np
from numpy import pi
import itertools
from pyramid.magdata import MagData
def simple_axis_projection(mag_data, axis='z', threshold=0):
'''
Project a magnetization distribution along one of the main axes of the 3D-grid.
Parameters
----------
mag_data : :class:`~pyramid.magdata.MagData`
A :class:`~pyramid.magdata.MagData` object storing the magnetization distribution,
which should be projected.
axis : {'z', 'y', 'x'}, optional
The projection direction as a string.
threshold : float, optional
A pixel only gets masked, if it lies above this threshold. The default is 0.
Returns
-------
projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
of the magnetization and the thickness projection for the resulting 2D-grid. The latter
has to be multiplied with the resolution for a value in nm.
'''
assert isinstance(mag_data, MagData), 'Parameter mag_data has to be a MagData object!'
assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as string)!'
if axis == 'z':
projection = (mag_data.magnitude[1].sum(0), # y_mag -> v_mag
mag_data.magnitude[2].sum(0), # x_mag -> u_mag
mag_data.get_mask(threshold).sum(0)) # thickness profile
elif axis == 'y':
projection = (mag_data.magnitude[0].sum(1), # z_mag -> v_mag
mag_data.magnitude[2].sum(1), # x_mag -> u_mag
mag_data.get_mask(threshold).sum(1)) # thickness profile
elif axis == 'x':
projection = (mag_data.magnitude[0].sum(2), # z_mag -> v_mag
mag_data.magnitude[1].sum(2), # y_mag -> u_mag
mag_data.get_mask(threshold).sum(2)) # thickness profile
return projection
def single_tilt_projection(mag_data, tilt=0, threshold=0):
# TODO: Docstring!!!
assert isinstance(mag_data, MagData), 'Parameter mag_data has to be a MagData object!'
# assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as string)!'
# Set starting variables:
dim = (mag_data.dim[0], mag_data.dim[2])
z_mag, y_mag, x_mag = mag_data.magnitude
mask = mag_data.get_mask()
projection = (np.zeros((mag_data.dim[1], mag_data.dim[2])),
np.zeros((mag_data.dim[1], mag_data.dim[2])),
np.zeros((mag_data.dim[1], mag_data.dim[2])))
def get_position(p, m, b, size):
x, y = 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):
a, b = delta-rho, delta+rho
if a >= 1 or b <= -1: # TODO: Should not be necessary!
print 'Outside of bounds:', delta
return 0
# Upper boundary:
if b >= 1:
w_b = 0.5
else:
w_b = (b*np.sqrt(1-b**2) + np.arctan(b/np.sqrt(1-b**2))) / pi
# Lower boundary:
if a <= -1:
w_a = -0.5
else:
w_a = (a*np.sqrt(1-a**2) + np.arctan(a/np.sqrt(1-a**2))) / pi
return w_b - w_a
# Creating coordinate list of all voxels:
xi = range(dim[0])
yj = range(dim[1])
ii, jj = np.meshgrid(xi, yj)
voxels = list(itertools.product(yj, xi))
# Calculate positions along the projected pixel coordinate system:
direction = (-np.cos(tilt), np.sin(tilt))
center = (dim[0]/2., dim[1]/2.)
m = direction[0]/(direction[1]+1E-30)
b = center[0] - m * center[1]
positions = get_position(voxels, m, b, dim[0])
# Calculate weights:
r = 1/np.sqrt(np.pi) # radius of the voxel circle
rho = 0.5 / r # TODO: ratio of radii
weights = []
for i, voxel in enumerate(voxels):
voxel_weights = []
impacts = get_impact(positions[i], r, dim[0])
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
voxel_weights.append((impact, get_weight(delta, rho)))
weights.append((voxel, voxel_weights))
# Calculate projection with the calculated weights for the voxels:
for i, weight in enumerate(weights):
voxel = weights[i][0]
voxel_weights = weights[i][1]
for voxel_weight in voxel_weights:
pixel, weight = voxel_weight
# Component parallel to tilt axis (':' goes over all slices):
projection[0][:, pixel] += weight * y_mag[voxel[0], :, voxel[1]]
# Component perpendicular to tilt axis:
projection[1][:, pixel] += weight * (x_mag[voxel[0], :, voxel[1]]*np.cos(tilt)
+ z_mag[voxel[0], :, voxel[1]]*np.sin(tilt))
# Thickness profile:
projection[2][:, pixel] += weight * mask[voxel[0], :, voxel[1]]
return projection
# -*- 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.
So far, only a simple least square algorithm for known pixel locations for 2-dimensional problems
is implemented (:func:`~.reconstruct_simple_leastsq`), but more complex solutions are planned.
"""
import numpy as np
from scipy.optimize import leastsq
import pyramid.projector as pj
import pyramid.phasemapper as pm
from pyramid.magdata import MagData
def reconstruct_simple_leastsq(phase_map, mask, b_0=1):
'''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.
Returns
-------
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a
:class:`~pyramid.magdata.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.reshape(-1) # Measured phase map as a vector
res = phase_map.res # Resolution
dim = mask.shape # Dimensions of the mag. distr.
count = mask.sum() # Number of pixels with magnetization
lam = 1e-6 # Regularisation parameter
# Create empty MagData object for the reconstruction:
mag_data_rec = MagData(res, (np.zeros(dim), np.zeros(dim), np.zeros(dim)))
# Function that returns the phase map for a magnetic configuration x:
def F(x):
mag_data_rec.set_vector(mask, x)
phase = pm.phase_mag_real(res, pj.simple_axis_projection(mag_data_rec), 'slab', b_0)
return phase.reshape(-1)
# Cost function which should be minimized:
def J(x_i):
y_i = F(x_i)
term1 = (y_i - y_m)
term2 = lam * x_i
return np.concatenate([term1, term2])
# Reconstruct the magnetization components:
x_rec, _ = leastsq(J, np.zeros(3*count))
mag_data_rec.set_vector(mask, x_rec)
return mag_data_rec
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.