diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 9e5e42b910de66622962fb56b71ebf7f45c3c758..b0c84773177ab8dda4180d53839e20174cce1d50 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -78,6 +78,7 @@ class Kernel(object):
# 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
@@ -88,7 +89,7 @@ class Kernel(object):
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) # TODO: v = np.swapaxes(u, 0,1)
+ self.v = coeff * get_elementary_phase(geometry, vv, uu, a)
# Calculate Fourier trafo of kernel components:
dim_combined = 3*np.array(dim_uv) - 1 # dim_uv + (2*dim_uv - 1) magnetisation + kernel
self.dim_fft = 2 ** np.ceil(np.log2(dim_combined)).astype(int) # next multiple of 2
@@ -98,8 +99,11 @@ class Kernel(object):
if numcore:
self.LOG.debug('numcore activated')
self._multiply_jacobi_method = self._multiply_jacobi_core
+ self._multiply_jacobi_T_method = self._multiply_jacobi_T_core
else:
+ self.LOG.debug('numcore deactivated')
self._multiply_jacobi_method = self._multiply_jacobi
+ self._multiply_jacobi_T_method = self._multiply_jacobi_T
self.LOG.debug('Created '+str(self))
def __repr__(self):
@@ -134,6 +138,8 @@ class Kernel(object):
'''
self.LOG.debug('Calling jac_dot')
+ assert len(vector) == 2 * self.size, \
+ 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.size)
return self._multiply_jacobi_method(vector)
def jac_T_dot(self, vector):
@@ -142,53 +148,53 @@ class Kernel(object):
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.
+ 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.
+ the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
'''
self.LOG.debug('Calling jac_dot_T')
- return self._multiply_jacobi_T(vector)
+ assert len(vector) == self.size, \
+ 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.size)
+ return self._multiply_jacobi_T_method(vector)
def _multiply_jacobi(self, vector):
self.LOG.debug('Calling _multiply_jacobi')
v_dim, u_dim = self.dim_uv
- size = np.prod(self.dim_uv)
- assert len(vector) == 2*size, 'vector size not compatible!'
- result = np.zeros(size)
- for s in range(size): # column-wise (two columns at a time, u- and v-component)
- i = s % u_dim
- j = int(s/u_dim)
- u_min = (u_dim-1) - i
+ assert len(vector) == 2 * self.size, \
+ 'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.size)
+ result = np.zeros(self.size)
+ # Iterate over all contributing pixels (numbered consecutively)
+ for s in range(self.size): # 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_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.u[v_min:v_max, u_min:u_max].reshape(-1) # u
- result -= vector[s+size] * self.v[v_min:v_max, u_min:u_max].reshape(-1) # v
+ result -= vector[s+self.size] * self.v[v_min:v_max, u_min:u_max].reshape(-1) # v
return result
def _multiply_jacobi_T(self, vector):
self.LOG.debug('Calling _multiply_jacobi_T')
v_dim, u_dim = self.dim_uv
- size = np.prod(self.dim_uv)
- assert len(vector) == size, \
- 'vector size not compatible! vector: {}, size: {}'.format(len(vector), size)
- result = np.zeros(2*size)
- for s in range(size): # row-wise (two rows at a time, u- and v-component)
- i = s % u_dim
- j = int(s/u_dim)
- u_min = (u_dim-1) - i
- u_max = (2*u_dim-1) - i
- v_min = (v_dim-1) - j
- v_max = (2*v_dim-1) - j
- result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1))
- result[s+size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))
+ result = np.zeros(2*self.size)
+ # Iterate over all contributing pixels (numbered consecutively):
+ for s in range(self.size): # 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)) # u
+ result[s+self.size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1)) # v
return result
def _multiply_jacobi_core(self, vector):
@@ -196,3 +202,29 @@ class Kernel(object):
result = np.zeros(np.prod(self.dim_uv))
core.multiply_jacobi_core(self.dim_uv[0], self.dim_uv[1], self.u, self.v, vector, result)
return result
+
+ def _multiply_jacobi_T_core(self, vector):
+ self.LOG.debug('Calling _multiply_jacobi_T_core')
+ result = np.zeros(2*np.prod(self.dim_uv))
+ core.multiply_jacobi_T_core(self.dim_uv[0], self.dim_uv[1], self.u, self.v, vector, result)
+ return result
+
+ def test(self, vector):
+ result = np.zeros(2*np.prod(self.dim_uv))
+ u_dim, v_dim = self.dim_uv
+ 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] * self.u[v, u]
+ result[r+self.size] -= vector[s] * self.v[v, u]
+ s += 1
+ r += 1
+ return result
+# TODO: delete
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index b7845ed76daea043829679da372773d4a8387f93..df365f9887bace0bb80dcaf9abcc01b7d54c0f59 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -2,13 +2,19 @@
"""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
@@ -506,3 +512,55 @@ class MagData(object):
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):
+ # TODO: Docstring!
+ 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))
+
+ 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')
+
+ 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))
+
+ tree.write(filename, pretty_print=True)
diff --git a/pyramid/numcore/kernel_core.pyx b/pyramid/numcore/kernel_core.pyx
index f0a7412f8ac73cdd78c0e0fb7820915dbf9c01eb..7eb5e84edfad64e1418b3beb206ff346d3ab4010 100644
--- a/pyramid/numcore/kernel_core.pyx
+++ b/pyramid/numcore/kernel_core.pyx
@@ -44,19 +44,69 @@ def multiply_jacobi_core(
iteration loop.
'''
- cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, ri, u, v, size
+ 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 pixel (numbered consecutively)
- # Go over all 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
- ri = 0 # Current result component (numbered consecutively)
+ 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[ri] += vector[s] * u_phi[v, u]
- result[ri] -= vector[s+size] * v_phi[v, u]
- ri += 1
+ 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 multiply_jacobi_T_core(
+ 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 in the :class:`~.Kernel` class.
+
+ 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?
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index 85b9f03c139cce4d52547b03645b6f4c8adcd3e2..a2f6f1a1300b71e161a3b787e47612bb0ff859cd 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -149,7 +149,7 @@ class PhaseMap(object):
def __repr__(self):
self.LOG.debug('Calling __repr__')
- return '%s(a=%r, phase=%r, unit=&r)' % \
+ return '%s(a=%r, phase=%r, unit=%r)' % \
(self.__class__, self.a, self.phase, self.unit)
def __str__(self):
@@ -180,7 +180,9 @@ class PhaseMap(object):
def __mul__(self, other): # self * other
self.LOG.debug('Calling __mul__')
- assert isinstance(other, Number), 'PhaseMap objects can only be multiplied by numbers!'
+ 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
@@ -397,7 +399,7 @@ class PhaseMap(object):
return axis
def display_holo(self, density=1, title='Holographic Contour Map',
- axis=None, grad_encode='dark', interpolation='none', show=True):
+ axis=None, grad_encode='bright', interpolation='none', show=True):
'''Display the color coded holography image.
Parameters
@@ -408,7 +410,7 @@ class PhaseMap(object):
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.
- grad_encode: {'dark', 'bright', 'color', 'none'}, optional
+ 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
@@ -458,8 +460,7 @@ class PhaseMap(object):
if axis is None:
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1, aspect='equal')
- # Plot the image on a black background and set axes:
- axis.patch.set_facecolor('black')
+ # 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)
@@ -478,7 +479,7 @@ class PhaseMap(object):
return axis
def display_combined(self, title='Combined Plot', cmap='RdBu', limit=None, norm=None,
- density=1, interpolation='none', grad_encode='dark'):
+ density=1, interpolation='none', grad_encode='bright'):
'''Display the phase map and the resulting color coded holography image in one plot.
Parameters
@@ -500,7 +501,7 @@ class PhaseMap(object):
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: {'dark', 'bright', 'color', 'none'}, optional
+ 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
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 20ec3d04c08efb9f73eb34ba4fac5d26457645a5..77ed1bed91e313ecac7705862c12d68b87094ae3 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -67,6 +67,9 @@ class PMAdapterFM(PhaseMapper):
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The 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
Elementary geometry which is used for the phase contribution of one pixel.
Default is 'disc'.
@@ -78,14 +81,15 @@ class PMAdapterFM(PhaseMapper):
LOG = logging.getLogger(__name__+'.PMAdapterFM')
- def __init__(self, a, projector, b_0=1, geometry='disc'):
+ def __init__(self, a, projector, b_0=1, numcore=True, geometry='disc'):
self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.b_0 = b_0
self.geometry = geometry
- self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, b_0, geometry))
+ self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, b_0,
+ numcore, geometry))
self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 6dca851d51f169cd9a0a57ad7e6f5f6f646caa76..758914840e59562fd22fd9909dca9e910554e805 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -291,7 +291,7 @@ class XTiltProjector(Projector):
Information about the projector as a string, e.g. for the use in plot titles.
'''
- return 'x-tilt: $(\phi = {:2.1f} \pi)$'.format(self.tilt)
+ return 'x-tilt: $(\phi = {:3.2f} \pi)$'.format(self.tilt/pi)
class YTiltProjector(Projector):
@@ -393,7 +393,7 @@ class YTiltProjector(Projector):
Information about the projector as a string, e.g. for the use in plot titles.
'''
- return 'y-tilt: $(\phi = {:2.1f} \pi)$'.format(self.tilt)
+ return 'y-tilt: $(\phi = {:3.2f} \pi)$'.format(self.tilt/pi)
class SimpleProjector(Projector):
@@ -459,4 +459,4 @@ class SimpleProjector(Projector):
Information about the projector as a string, e.g. for the use in plot titles.
'''
- return 'projected along {}-axis'.format(self.tilt)
+ return 'projected along {}-axis'.format(self.axis)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 04bdc208e316eeb06c0a7b0ceb89f0b8bf678fc0..c7a7e513b6d6b724fc8840572c509664a1608be9 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -38,10 +38,10 @@ class PrintIterator(object):
: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 : {2, 1, 0}, optional
- Parameter defining the verbosity 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.
+ 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
-----
@@ -52,7 +52,7 @@ class PrintIterator(object):
LOG = logging.getLogger(__name__+'.PrintIterator')
- def __init__(self, cost, verbosity=2):
+ def __init__(self, cost, verbosity):
self.LOG.debug('Calling __init__')
self.cost = cost
self.verbosity = verbosity
@@ -83,7 +83,7 @@ class PrintIterator(object):
return self.iteration
-def optimize_sparse_cg(data, verbosity=2):
+def optimize_sparse_cg(data, verbosity=0):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
@@ -93,10 +93,10 @@ def optimize_sparse_cg(data, verbosity=2):
: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.
+ 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
-------
@@ -113,7 +113,7 @@ def optimize_sparse_cg(data, verbosity=2):
# Optimize:
A = CFAdapterScipyCG(cost)
b = fwd_model.jac_T_dot(None, y)
- x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity))
+ x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity), maxiter=1000)
# Create and return fitting MagData object:
mag_opt = MagData(fwd_model.a, np.zeros((3,)+fwd_model.dim))
mag_opt.mag_vec = x_opt
@@ -150,7 +150,7 @@ def optimize_cg(data, first_guess):
cost = Costfunction(y, fwd_model)
# Optimize:
result = minimize(cost, x_0, method='Newton-CG', jac=cost.jac, hessp=cost.hess_dot,
- options={'maxiter': 200, 'disp': True})
+ 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))
diff --git a/pyramid/template.x3d b/pyramid/template.x3d
new file mode 100644
index 0000000000000000000000000000000000000000..1ab06c35d6d3aaf87748740071ce0c94a48552c6
--- /dev/null
+++ b/pyramid/template.x3d
@@ -0,0 +1,152 @@
+<!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>
diff --git a/scripts/collaborations/rueffner_file.py b/scripts/collaborations/rueffner_file.py
new file mode 100644
index 0000000000000000000000000000000000000000..5a5b501f988392c58d34e8f1c936097c8134f589
--- /dev/null
+++ b/scripts/collaborations/rueffner_file.py
@@ -0,0 +1,126 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jul 25 14:37:11 2014
+
+@author: Jan
+"""
+
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+from mayavi import mlab
+from pylab import griddata
+
+from tqdm import tqdm
+
+import tables
+
+import pyramid
+from pyramid.magdata import MagData
+from pyramid.projector import SimpleProjector
+from pyramid.phasemapper import PMConvolve
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+###################################################################################################
+PATH = '../../output/'
+
+#mag_file = netCDF4.Dataset(PATH+'dyn_0090mT_dyn.h5_dat.h5')
+#
+#print mag_file
+#
+#print mag_file.groups['data'].groups['fields']
+#
+#data = mag_file.groups['data']
+#
+#fields = data.groups['fields']
+
+
+#points = netCDF4.Dataset(PATH+'dyn_0090mT_dyn.h5_dat.h5').groups['mesh'].variables['points'][...]
+
+dim = (16, 182+40, 210+40)
+a = 1.
+b_0 = 1.1
+projector = SimpleProjector(dim)
+pm = PMConvolve(a, projector, b_0)
+
+h5file = tables.openFile(PATH+'dyn_0090mT_dyn.h5_dat.h5')
+points = h5file.root.mesh.points.read()
+
+axis = plt.figure().add_subplot(1, 1, 1, aspect='equal')
+axis.scatter(points[:, 0], points[:, 1])
+mlab.points3d(points[:, 0], points[:, 1], points[:, 2], mode='2dvertex')
+
+#data_old = h5file.root.data.fields.m.read(field='m_CoFeb')[0, ...]
+
+# Filling zeros:
+iz_x = np.concatenate([np.linspace(-74, -37, 20),
+ np.linspace(-37, 37, 20),
+ np.linspace(37, 74, 20),
+ np.linspace(74, 37, 20),
+ np.linspace(37, -37, 20),
+ np.linspace(-37, -74, 20)])
+iz_y = np.concatenate([np.linspace(0, 64, 20),
+ np.linspace(64, 64, 20),
+ np.linspace(64, 0, 20),
+ np.linspace(0, -64, 20),
+ np.linspace(-64, -64, 20),
+ np.linspace(-64, 0, 20)])
+iz_z = np.zeros(len(iz_x))
+
+# Set up grid:
+xs, ys, zs = np.arange(-105, 105), np.arange(-91, 91), np.arange(-8, 8)
+xx, yy = np.meshgrid(xs, ys)
+
+#test = []
+
+for t in np.arange(900, 1001, 5):
+ data = h5file.root.data.fields.m.read(field='m_CoFeb')[t, ...]
+# if data_old is not None:
+# np.testing.assert_equal(data, data_old)
+# data_old = np.copy(data)
+
+ zs_unique = np.unique(points[:, 2])
+
+ # Create empty magnitude:
+ magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
+
+ for i, z in tqdm(enumerate(zs), total=len(zs)):
+ # print z
+ # print a/2
+ # print np.abs(data[:, 2]-z)
+ z_slice = data[np.abs(points[:, 2]-z) <= a/2., :]
+ weights = 1 - np.abs(z_slice[:, 2]-z)*2/a
+ # print z_slice.shape, z
+ # print z_slice
+ # print weights
+ for j in range(3): # For all 3 components!
+ # if z <= 0:
+ grid_x = np.concatenate([z_slice[:, 0], iz_x])
+ grid_y = np.concatenate([z_slice[:, 1], iz_y])
+ grid_z = np.concatenate([weights*z_slice[:, j], iz_z])
+ # else:
+ # grid_x = z_slice[:, 0]
+ # grid_y = z_slice[:, 1]
+ # grid_z = weights*z_slice[:, 3+j]
+ grid = griddata(grid_x, grid_y, grid_z, xx, yy)
+ magnitude[j, i, :, :] = grid.filled(fill_value=0)
+
+ mag_data = MagData(a, magnitude)
+ mag_data.pad(20, 20, 0)
+ phase_map = pm(mag_data)
+ phase_map.display_combined(density=250, interpolation='bilinear', limit=0.25,
+ grad_encode='bright')[0]
+ plt.savefig(PATH+'rueffner/phase_map_t_'+str(t)+'.png')
+ plt.close('all')
+# mag_data.quiver_plot()
+# mag_data.save_to_x3d('rueffner.x3d')
+# mag_data.scale_down()
+# mag_data.quiver_plot3d()
diff --git a/scripts/collaborations/vtk_conversion.py b/scripts/collaborations/vtk_conversion.py
index ea5e99b492fdc21a8d90cc433f1af01626c0b4bb..d328fc41962ed30415773295010c936442e5c2e9 100644
--- a/scripts/collaborations/vtk_conversion.py
+++ b/scripts/collaborations/vtk_conversion.py
@@ -5,8 +5,10 @@
import os
import numpy as np
+from numpy import pi
import matplotlib.pyplot as plt
from pylab import griddata
+from mayavi import mlab
import pickle
import vtk
@@ -14,8 +16,9 @@ from tqdm import tqdm
import pyramid
from pyramid.magdata import MagData
-from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMAdapterFM
+from pyramid.projector import SimpleProjector, YTiltProjector, XTiltProjector
+from pyramid.phasemapper import PMAdapterFM, PMConvolve
+from pyramid.dataset import DataSet
import logging
import logging.config
@@ -26,9 +29,9 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
-PATH = '../../output/vtk data/tube_90x30x50_sat_edge_equil.gmr'
+PATH = '../../output/vtk data/longtube_withcap/CoFeB_tube_cap_4nm'
b_0 = 1.54
-gain = 12
+gain = 1
force_calculation = False
###################################################################################################
# Load vtk-data:
@@ -61,50 +64,110 @@ else:
with open(PATH+'.pickle') as pf:
data = pickle.load(pf)
# Scatter plot of all x-y-coordinates
-axis = plt.figure().add_subplot(1, 1, 1)
-axis.scatter(data[:, 0], data[:, 1])
+axis = plt.figure().add_subplot(1, 1, 1, aspect='equal')
+axis.scatter(data[data[:, 2] <= 0, 0], data[data[:, 2] <= 0, 1])
+axis
+mlab.points3d(data[:, 0], data[:, 1], data[:, 2], mode='2dvertex')
plt.show()
+
###################################################################################################
# Interpolate on regular grid:
if force_calculation or not os.path.exists(PATH+'.nc'):
- # Find unique z-slices:
- zs = np.unique(data[:, 2])
- # Determine the grid spacing:
- a = zs[1] - zs[0]
+
# Determine the size of object:
x_min, x_max = data[:, 0].min(), data[:, 0].max()
y_min, y_max = data[:, 1].min(), data[:, 1].max()
z_min, z_max = data[:, 2].min(), data[:, 2].max()
x_diff, y_diff, z_diff = np.ptp(data[:, 0]), np.ptp(data[:, 1]), np.ptp(data[:, 2])
x_cent, y_cent, z_cent = x_min+x_diff/2., y_min+y_diff/2., z_min+z_diff/2.
- # Create regular grid
+ # Filling zeros:
+ iz_x = np.concatenate([np.linspace(-2.95, -2.95, 20),
+ np.linspace(-2.95, 0, 20),
+ np.linspace(0, 2.95, 20),
+ np.linspace(2.95, 2.95, 20),
+ np.linspace(2.95, 0, 20),
+ np.linspace(0, -2.95, 20), ])
+ iz_y = np.concatenate([np.linspace(-1.70, 1.70, 20),
+ np.linspace(1.70, 3.45, 20),
+ np.linspace(3.45, 1.70, 20),
+ np.linspace(1.70, -1.70, 20),
+ np.linspace(-1.70, -3.45, 20),
+ np.linspace(-3.45, -1.70, 20), ])
+ # Find unique z-slices:
+ zs_unique = np.unique(data[:, 2])
+ # Determine the grid spacing (in 1/10 nm):
+ a = (zs_unique[1] - zs_unique[0])
+ # Set actual z-slices:
+ zs = np.arange(z_min, z_max, a)
+ # Create regular grid:
xs = np.arange(x_cent-x_diff, x_cent+x_diff, a)
ys = np.arange(y_cent-y_diff, y_cent+y_diff, a)
- o, p = np.meshgrid(xs, ys)
+ xx, yy = np.meshgrid(xs, ys)
# Create empty magnitude:
magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
- # WITH MASKING OF THE CENTER (SYMMETRIC):
- iz_x = np.concatenate([np.linspace(-4.95, -4.95, 50),
- np.linspace(-4.95, 0, 50),
- np.linspace(0, 4.95, 50),
- np.linspace(4.95, 4.95, 50),
- np.linspace(-4.95, 0, 50),
- np.linspace(0, 4.95, 50), ])
- iz_y = np.concatenate([np.linspace(-2.88, 2.88, 50),
- np.linspace(2.88, 5.7, 50),
- np.linspace(5.7, 2.88, 50),
- np.linspace(2.88, -2.88, 50),
- np.linspace(-2.88, -5.7, 50),
- np.linspace(-5.7, -2.88, 50), ])
for i, z in tqdm(enumerate(zs), total=len(zs)):
- subdata = data[data[:, 2] == z, :]
+# n = bisect.bisect(zs_unique, z)
+# z_lo, z_up = zs_unique[n], zs_unique[n+1]
+ z_slice = data[np.abs(data[:, 2]-z) <= a/2., :]
+ weights = 1 - np.abs(z_slice[:, 2]-z)*2/a
+# print n, z_slice.shape, z, z_lo, z_up
+# print z_slice
+# print weights
for j in range(3): # For all 3 components!
- gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
- np.concatenate([subdata[:, 1], iz_y]),
- np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
- o, p)
- magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
+ if z <= 0:
+ grid_x = np.concatenate([z_slice[:, 0], iz_x])
+ grid_y = np.concatenate([z_slice[:, 1], iz_y])
+ grid_z = np.concatenate([weights*z_slice[:, 3+j], np.zeros(len(iz_x))])
+ else:
+ grid_x = z_slice[:, 0]
+ grid_y = z_slice[:, 1]
+ grid_z = weights*z_slice[:, 3+j]
+ grid = griddata(grid_x, grid_y, grid_z, xx, yy)
+ magnitude[j, i, :, :] = grid.filled(fill_value=0)
+
+ a = a*10 # convert to nm
+# print a
+
+# for i, z in tqdm(enumerate(zs), total=len(zs)):
+# n = bisect.bisect(zs_unique, z)
+# z_lo, z_up = zs_unique[n], zs_unique[n+1]
+# slice_lo = data[data[:, 2] == z_lo, :]
+# slice_up = data[data[:, 2] == z_up, :]
+# print n, slice_up.shape, z, z_lo, z_up
+# print slice_up
+# if slice_up.shape[0] < 3:
+# continue
+# for j in range(3): # For all 3 components!
+# # Lower layer:
+# grid_lo = griddata(slice_lo[:, 0], slice_lo[:, 1], slice_lo[:, 3 + j], xx, yy)
+# # Upper layer:
+# grid_up = griddata(slice_up[:, 0], slice_up[:, 1], slice_up[:, 3 + j], xx, yy)
+# # Interpolated layer:
+# grid_interpol = (z-z_lo)/(z_up-z_lo)*grid_lo + (z_up-z)/(z_up-z_lo)*grid_up
+# magnitude[j, i
+
+# # WITH MASKING OF THE CENTER (SYMMETRIC):
+# iz_x = np.concatenate([np.linspace(-2.95, -2.95, 20),
+# np.linspace(-2.95, 0, 20),
+# np.linspace(0, 2.95, 20),
+# np.linspace(2.95, 2.95, 20),
+# np.linspace(2.95, 0, 20),
+# np.linspace(0, -2.95, 20), ])
+# iz_y = np.concatenate([np.linspace(-1.70, 1.70, 20),
+# np.linspace(1.70, 3.45, 20),
+# np.linspace(3.45, 1.70, 20),
+# np.linspace(1.70, -1.70, 20),
+# np.linspace(-1.70, -3.45, 20),
+# np.linspace(-3.45, -1.70, 20), ])
+# for i, z in tqdm(enumerate(zs), total=len(zs)):
+# subdata = data[data[:, 2] == z, :]
+# for j in range(3): # For all 3 components!
+# gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
+# np.concatenate([subdata[:, 1], iz_y]),
+# np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
+# o, p)
+# magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
# # WITH MASKING OF THE CENTER (ASYMMETRIC):
# iz_x = np.concatenate([np.linspace(-5.88, -5.88, 50),
@@ -131,16 +194,20 @@ if force_calculation or not os.path.exists(PATH+'.nc'):
# # WITHOUT MASKING OF THE CENTER:
# for i, z in tqdm(enumerate(zs), total=len(zs)):
# subdata = data[data[:, 2] == z, :]
+# print subdata.shape, z, zs_temp
+# if subdata.shape[0] < 3:
+# continue
# for j in range(3): # For all 3 components!
# gridded_subdata = griddata(subdata[:, 0], subdata[:, 1], subdata[:, 3 + j], o, p)
# magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
# Creating MagData object:
- mag_data = MagData(0.2*10, magnitude)
+ mag_data = MagData(a, np.pad(magnitude, ((0, 0), (0, 0), (0, 1), (0, 1)), mode='constant',
+ constant_values=0))
mag_data.save_to_netcdf4(PATH+'.nc')
else:
mag_data = MagData.load_from_netcdf4(PATH+'.nc')
-mag_data.quiver_plot()
+mag_data.quiver_plot(proj_axis='x')
###################################################################################################
# Phasemapping:
projector = SimpleProjector(mag_data.dim)
@@ -152,3 +219,28 @@ plt.savefig(PATH+'_{}T_cosx{}.png'.format(b_0, gain))
(-phase_map).display_combined(title=r'Combined Plot (B$_0$={} T, Cos x {})'.format(b_0, gain),
density=gain, interpolation='bilinear')
plt.savefig(PATH+'_{}T_cosx{}_smooth.png'.format(b_0, gain))
+mag_data.scale_down()
+mag_data.save_to_netcdf4(PATH+'_scaled.nc')
+mag_data_tip = MagData(mag_data.a, mag_data.magnitude[:, 350:, :, :])
+mag_data_tip.save_to_netcdf4(PATH+'_tip_scaled.nc')
+mag_data_tip.quiver_plot3d()
+mag_data_tip.pad(20, 20, 0)
+
+count = 16
+
+tilts_full = np.linspace(-pi/2, pi/2, num=count/2, endpoint=False)
+tilts_miss = np.linspace(-pi/3, pi/3, num=count/2, endpoint=False)
+
+projectors_y = [YTiltProjector(mag_data_tip.dim, tilt) for tilt in tilts_miss]
+projectors_x = [XTiltProjector(mag_data_tip.dim, tilt) for tilt in tilts_miss]
+projectors = np.concatenate((projectors_y, projectors_x))
+phasemappers = [PMConvolve(mag_data.a, proj, b_0) for proj in projectors]
+
+data_set = DataSet(mag_data_tip.a, mag_data_tip.dim[1:], b_0)
+
+gain = 10
+
+for i, pm in enumerate(phasemappers):
+ data_set.append((pm(mag_data_tip), projectors[i]))
+
+data_set.display()
diff --git a/scripts/collaborations/vtk_conversion_nanowire_segment.py b/scripts/collaborations/vtk_conversion_nanowire_segment.py
new file mode 100644
index 0000000000000000000000000000000000000000..ea5e99b492fdc21a8d90cc433f1af01626c0b4bb
--- /dev/null
+++ b/scripts/collaborations/vtk_conversion_nanowire_segment.py
@@ -0,0 +1,154 @@
+# -*- coding: utf-8 -*-
+"""Created on Fri Jan 24 11:17:11 2014 @author: Jan"""
+
+
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+from pylab import griddata
+
+import pickle
+import vtk
+from tqdm import tqdm
+
+import pyramid
+from pyramid.magdata import MagData
+from pyramid.projector import SimpleProjector
+from pyramid.phasemapper import PMAdapterFM
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+###################################################################################################
+PATH = '../../output/vtk data/tube_90x30x50_sat_edge_equil.gmr'
+b_0 = 1.54
+gain = 12
+force_calculation = False
+###################################################################################################
+# Load vtk-data:
+if force_calculation or not os.path.exists(PATH+'.pickle'):
+ # Setting up reader:
+ reader = vtk.vtkDataSetReader()
+ reader.SetFileName(PATH+'.vtk')
+ reader.ReadAllScalarsOn()
+ reader.ReadAllVectorsOn()
+ reader.Update()
+ # Getting output:
+ output = reader.GetOutput()
+ # Reading points and vectors:
+ size = output.GetNumberOfPoints()
+ vtk_points = output.GetPoints().GetData()
+ vtk_vectors = output.GetPointData().GetVectors()
+ # Converting points to numpy array:
+ point_array = np.zeros(vtk_points.GetSize())
+ vtk_points.ExportToVoidPointer(point_array)
+ point_array = np.reshape(point_array, (-1, 3))
+ # Converting vectors to numpy array:
+ vector_array = np.zeros(vtk_points.GetSize())
+ vtk_vectors.ExportToVoidPointer(vector_array)
+ vector_array = np.reshape(vector_array, (-1, 3))
+ # Combining data:
+ data = np.hstack((point_array, vector_array))
+ with open(PATH+'.pickle', 'w') as pf:
+ pickle.dump(data, pf)
+else:
+ with open(PATH+'.pickle') as pf:
+ data = pickle.load(pf)
+# Scatter plot of all x-y-coordinates
+axis = plt.figure().add_subplot(1, 1, 1)
+axis.scatter(data[:, 0], data[:, 1])
+plt.show()
+###################################################################################################
+# Interpolate on regular grid:
+if force_calculation or not os.path.exists(PATH+'.nc'):
+ # Find unique z-slices:
+ zs = np.unique(data[:, 2])
+ # Determine the grid spacing:
+ a = zs[1] - zs[0]
+ # Determine the size of object:
+ x_min, x_max = data[:, 0].min(), data[:, 0].max()
+ y_min, y_max = data[:, 1].min(), data[:, 1].max()
+ z_min, z_max = data[:, 2].min(), data[:, 2].max()
+ x_diff, y_diff, z_diff = np.ptp(data[:, 0]), np.ptp(data[:, 1]), np.ptp(data[:, 2])
+ x_cent, y_cent, z_cent = x_min+x_diff/2., y_min+y_diff/2., z_min+z_diff/2.
+ # Create regular grid
+ xs = np.arange(x_cent-x_diff, x_cent+x_diff, a)
+ ys = np.arange(y_cent-y_diff, y_cent+y_diff, a)
+ o, p = np.meshgrid(xs, ys)
+ # Create empty magnitude:
+ magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
+
+ # WITH MASKING OF THE CENTER (SYMMETRIC):
+ iz_x = np.concatenate([np.linspace(-4.95, -4.95, 50),
+ np.linspace(-4.95, 0, 50),
+ np.linspace(0, 4.95, 50),
+ np.linspace(4.95, 4.95, 50),
+ np.linspace(-4.95, 0, 50),
+ np.linspace(0, 4.95, 50), ])
+ iz_y = np.concatenate([np.linspace(-2.88, 2.88, 50),
+ np.linspace(2.88, 5.7, 50),
+ np.linspace(5.7, 2.88, 50),
+ np.linspace(2.88, -2.88, 50),
+ np.linspace(-2.88, -5.7, 50),
+ np.linspace(-5.7, -2.88, 50), ])
+ for i, z in tqdm(enumerate(zs), total=len(zs)):
+ subdata = data[data[:, 2] == z, :]
+ for j in range(3): # For all 3 components!
+ gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
+ np.concatenate([subdata[:, 1], iz_y]),
+ np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
+ o, p)
+ magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
+
+# # WITH MASKING OF THE CENTER (ASYMMETRIC):
+# iz_x = np.concatenate([np.linspace(-5.88, -5.88, 50),
+# np.linspace(-5.88, 0, 50),
+# np.linspace(0, 5.88, 50),
+# np.linspace(5.88, 5.88, 50),
+# np.linspace(5.88, 0, 50),
+# np.linspace(0, -5.88, 50),])
+# iz_y = np.concatenate([np.linspace(-2.10, 4.50, 50),
+# np.linspace(4.50, 7.90, 50),
+# np.linspace(7.90, 4.50, 50),
+# np.linspace(4.50, -2.10, 50),
+# np.linspace(-2.10, -5.50, 50),
+# np.linspace(-5.50, -2.10, 50), ])
+# for i, z in tqdm(enumerate(zs), total=len(zs)):
+# subdata = data[data[:, 2] == z, :]
+# for j in range(3): # For all 3 components!
+# gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
+# np.concatenate([subdata[:, 1], iz_y]),
+# np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
+# o, p)
+# magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
+
+# # WITHOUT MASKING OF THE CENTER:
+# for i, z in tqdm(enumerate(zs), total=len(zs)):
+# subdata = data[data[:, 2] == z, :]
+# for j in range(3): # For all 3 components!
+# gridded_subdata = griddata(subdata[:, 0], subdata[:, 1], subdata[:, 3 + j], o, p)
+# magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
+
+ # Creating MagData object:
+ mag_data = MagData(0.2*10, magnitude)
+ mag_data.save_to_netcdf4(PATH+'.nc')
+else:
+ mag_data = MagData.load_from_netcdf4(PATH+'.nc')
+mag_data.quiver_plot()
+###################################################################################################
+# Phasemapping:
+projector = SimpleProjector(mag_data.dim)
+phasemapper = PMAdapterFM(mag_data.a, projector)
+phase_map = phasemapper(mag_data)
+(-phase_map).display_combined(title=r'Combined Plot (B$_0$={} T, Cos x {})'.format(b_0, gain),
+ density=gain)
+plt.savefig(PATH+'_{}T_cosx{}.png'.format(b_0, gain))
+(-phase_map).display_combined(title=r'Combined Plot (B$_0$={} T, Cos x {})'.format(b_0, gain),
+ density=gain, interpolation='bilinear')
+plt.savefig(PATH+'_{}T_cosx{}_smooth.png'.format(b_0, gain))
diff --git a/scripts/create distributions/create_core_shell_sphere.py b/scripts/create distributions/create_core_shell_sphere.py
index 9c99de55d1ef62e219323f8c262c2d8ce09d1eb7..d2dbe8a8335c030dc9566c9666882b678aad44b3 100644
--- a/scripts/create distributions/create_core_shell_sphere.py
+++ b/scripts/create distributions/create_core_shell_sphere.py
@@ -28,7 +28,7 @@ if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_core_shell_sphere.txt'
-a = 50.0 # in nm
+a = 1.0 # in nm
density = 500
dim = (32, 32, 32)
center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_test.py b/scripts/reconstruction/reconstruction_sparse_cg_test.py
index 37045488c9ad3a0d03936d9025aed8a7c398fc10..c5a6fa1378030e5f877840238f66650f1b12e134 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_test.py
+++ b/scripts/reconstruction/reconstruction_sparse_cg_test.py
@@ -7,8 +7,11 @@ import os
import numpy as np
from numpy import pi
+from time import clock
+
import pyramid
from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
from pyramid.projector import YTiltProjector, XTiltProjector
from pyramid.phasemapper import PMConvolve
from pyramid.dataset import DataSet
@@ -33,14 +36,23 @@ print('--Generating input phase_maps')
a = 10.
b_0 = 1.
-dim = (16, 16, 16)
+dim = (8, 32, 32)
dim_uv = dim[1:3]
-count = 32
+count = 8
+
+#center = (int(dim[0]/2)-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+#mag_shape = mc.Shapes.disc(dim, center, dim[2]/4, dim[2]/2)
+#magnitude = mc.create_mag_dist_vortex(mag_shape)
+#mag_data = MagData(a, magnitude)
+#mag_data.quiver_plot3d()
-mag_shape = mc.Shapes.disc(dim, (7.5, 7.5, 7.5), dim[2]/4, dim[2]/4)
-magnitude = mc.create_mag_dist_vortex(mag_shape)
+vortex_shape = mc.Shapes.disc(dim, (3.5, 9.5, 9.5), 5, 4)
+sphere_shape = mc.Shapes.sphere(dim, (3.5, 22.5, 9.5), 3)
+slab_shape = mc.Shapes.slab(dim, (3.5, 15.5, 22.5), (5, 8, 8))
+mag_data = MagData(a, mc.create_mag_dist_vortex(vortex_shape, (3.5, 9.5, 9.5)))
+mag_data += MagData(a, mc.create_mag_dist_homog(sphere_shape, pi/4, pi/4))
+mag_data += MagData(a, mc.create_mag_dist_homog(slab_shape, -pi/6))
-mag_data = MagData(a, magnitude)
mag_data.quiver_plot3d()
tilts_full = np.linspace(-pi/2, pi/2, num=count/2, endpoint=False)
@@ -59,10 +71,15 @@ projectors_xy_miss.extend(projectors_y_miss)
###################################################################################################
projectors = projectors_xy_full
+noise = 0.0
###################################################################################################
phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
-phase_maps = [pm(mag_data) for pm in phasemappers]
+if noise != 0:
+ phase_maps = [pm(mag_data) + PhaseMap(a, np.random.normal(0, noise, dim_uv))
+ for pm in phasemappers]
+else:
+ phase_maps = [pm(mag_data) for pm in phasemappers]
###################################################################################################
print('--Setting up data collection')
@@ -78,10 +95,14 @@ kern = Kernel(data.a, data.dim_uv, data.b_0)
F = ForwardModel(data.projectors, kern)
C = Costfunction(y, F, lam)
-###################################################################################################
-print('--Test simple solver')
+data.display()
-mag_data_opt = rc.optimize_sparse_cg(data)
-mag_data_opt.quiver_plot3d()
-(mag_data_opt - mag_data).quiver_plot3d()
-#data.display(data.create_phase_maps(mag_data_opt))
+###################################################################################################
+#print('--Test simple solver')
+#
+#start = clock()
+#mag_data_opt = rc.optimize_sparse_cg(data, verbosity=1)
+#print 'Time:', str(clock()-start)
+#mag_data_opt.quiver_plot3d()
+#(mag_data_opt - mag_data).quiver_plot3d()
+##data.display(data.create_phase_maps(mag_data_opt))
diff --git a/scripts/test methods/compare_methods.py b/scripts/test methods/compare_methods.py
index f7e172f60784b23451cbb7097681ed34988f34af..c4a5fbb7912337acb124dc45bc8c1c7e6311ef99 100644
--- a/scripts/test methods/compare_methods.py
+++ b/scripts/test methods/compare_methods.py
@@ -30,7 +30,7 @@ b_0 = 1.0 # in T
a = 1.0 # in nm
phi = pi/4
numcore = False
-padding = 5
+padding = 10
density = 10
dim = (128, 128, 128) # in px (z, y, x)
diff --git a/tests/test_analytic.py b/tests/test_analytic.py
index fd2784c11416971dbbd011ce7a641433fbf86790..40b990b02db1625c0478ab5976efb695348f5055 100644
--- a/tests/test_analytic.py
+++ b/tests/test_analytic.py
@@ -11,6 +11,7 @@ import pyramid.analytic as an
class TestCaseAnalytic(unittest.TestCase):
+ """TestCase for the analytic module."""
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_analytic/')
@@ -29,12 +30,14 @@ class TestCaseAnalytic(unittest.TestCase):
self.radius = None
def test_phase_mag_slab(self):
+ '''Test of the phase_mag_slab method.'''
width = (self.dim[0]/2, self.dim[1]/2, self.dim[2]/2)
phase = an.phase_mag_slab(self.dim, self.res, self.phi, self.center, width)
reference = np.load(os.path.join(self.path, 'ref_phase_slab.npy'))
np.testing.assert_equal(phase, reference, 'Unexpected behavior in phase_mag_slab()')
def test_phase_mag_disc(self):
+ '''Test of the phase_mag_disc method.'''
radius = self.dim[2]/4
height = self.dim[2]/2
phase = an.phase_mag_disc(self.dim, self.res, self.phi, self.center, radius, height)
@@ -42,12 +45,14 @@ class TestCaseAnalytic(unittest.TestCase):
np.testing.assert_equal(phase, reference, 'Unexpected behavior in phase_mag_disc()')
def test_phase_mag_sphere(self):
+ '''Test of the phase_mag_sphere method.'''
radius = self.dim[2]/4
phase = an.phase_mag_sphere(self.dim, self.res, self.phi, self.center, radius)
reference = np.load(os.path.join(self.path, 'ref_phase_sphere.npy'))
np.testing.assert_equal(phase, reference, 'Unexpected behavior in phase_mag_sphere()')
def test_phase_mag_vortex(self):
+ '''Test of the phase_mag_vortex method.'''
radius = self.dim[2]/4
height = self.dim[2]/2
phase = an.phase_mag_vortex(self.dim, self.res, self.center, radius, height)
diff --git a/tests/test_compliance.py b/tests/test_compliance.py
index 99e6f060bbe605f9a3958fe54ea62ddb649f29d2..943c727765f858865959ea465f3c787f763610a1 100644
--- a/tests/test_compliance.py
+++ b/tests/test_compliance.py
@@ -13,7 +13,7 @@ import pep8
class TestCaseCompliance(unittest.TestCase):
- """Class for checking compliance of pyramid.""" # TODO: Docstring
+ """TestCase for checking the pep8 compliance of the pyramid package."""
def setUp(self):
self.path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] # Pyramid dir
@@ -31,7 +31,7 @@ class TestCaseCompliance(unittest.TestCase):
return filepaths
def test_pep8(self):
- # TODO: Docstring
+ '''Test for pep8 compliance.'''
files = self.get_files_to_check(os.path.join(self.path, 'pyramid')) \
+ self.get_files_to_check(os.path.join(self.path, 'scripts')) \
+ self.get_files_to_check(os.path.join(self.path, 'tests'))