Skip to content
Snippets Groups Projects
Commit ec080940 authored by Jan Caron's avatar Jan Caron
Browse files

First simple reconstruction approach 

documentation: minor changes and updates of docstrings
parent 355a7e41
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 588 additions and 251 deletions
pyramid/__init__.py
+ 1
− 1
View file @ ec080940
......@@ -8,7 +8,7 @@ magcreator
magdata
Class for the storage of magnetization data.
projector
Class for projections of given magnetization distribution.
Class for projecting given magnetization distributions.
kernel
Class for the kernel matrix representing one magnetized pixel.
phasemapper
......
pyramid/costfunction.py
+ 12
− 3
View file @ ec080940
......@@ -26,20 +26,29 @@ class Costfunction:
y = self.y
F = self.F
Se_inv = self.Se_inv
return (F(x)-y).dot(Se_inv.dot(F(x)-y))
# print 'Costfunction - __call__ - input: ', len(x)
result = (F(x)-y).dot(Se_inv.dot(F(x)-y))
# print 'Costfunction - __call__ - output:', result
return result
def jac(self, x):
# TODO: Docstring!
y = self.y
F = self.F
Se_inv = self.Se_inv
return F.jac_T_dot(x, (Se_inv.dot(F(x)-y)))
# print 'Costfunction - jac - input: ', len(x)
result = F.jac_T_dot(x, Se_inv.dot(F(x)-y))
# print 'Costfunction - jac - output:', len(result)
return result
def hess_dot(self, x, vector):
# TODO: Docstring!
F = self.F
Se_inv = self.Se_inv
return F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector)))
# print 'Costfunction - hess_dot - input: ', len(vector)
result = F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector)))
# print 'Costfunction - hess_dot - output:', len(result)
return result
# TODO: Murks!
......
pyramid/datacollection.py
+ 17
− 14
View file @ ec080940
......@@ -17,8 +17,8 @@ class DataCollection(object):
'''Class for collecting phase maps and corresponding projectors.
Represents a collection of (e.g. experimentally derived) phase maps, stored as
:class:´~.PhaseMap´ objects and corresponding projectors stored as :class:`~.Projector`
objects. At creation, the grid spacing `a` and the dimension `dim` of the projected grid.
:class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
objects. At creation, the grid spacing `a` and the dimension `dim_uv` of the projected grid.
Data can be added via the :func:`~.append` method, where a :class:`~.PhaseMap` and a
:class:`~.Projector` have to be given as tuple argument.
......@@ -26,7 +26,7 @@ class DataCollection(object):
----------
a: float
The grid spacing in nm.
dim: tuple (N=2)
dim_uv: tuple (N=2)
Dimensions of the projected grid.
phase_maps:
A list of all stored :class:`~.PhaseMap` objects.
......@@ -42,8 +42,8 @@ class DataCollection(object):
return self._a
@property
def dim(self):
return self._dim
def dim_uv(self):
return self._dim_uv
@property
def phase_vec(self):
......@@ -57,25 +57,28 @@ class DataCollection(object):
def projectors(self):
return [d[1] for d in self.data]
def __init__(self, a, dim):
def __init__(self, a, dim_uv, b_0):
self.log = logging.getLogger(__name__)
self.log.info('Creating '+str(self))
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
self._a = a
assert isinstance(dim, tuple) and len(dim)==2, 'Dimension has to be a tuple of length 2!'
assert isinstance(dim_uv, tuple) and len(dim_uv)==2, \
'Dimension has to be a tuple of length 2!'
assert np.all([])
self.dim = dim
self._dim_uv = dim_uv
self.b_0 = b_0
self.data = []
# TODO: make it work
# self.log.info('Created:', str(self))
def __repr__(self):
self.log.info('Calling __repr__')
return '%s(a=%r, dim=%r)' % (self.__class__, self.a, self.dim)
return '%s(a=%r, dim_uv=%r)' % (self.__class__, self.a, self.dim_uv)
def __str__(self):
self.log.info('Calling __str__')
return 'DataCollection(a=%s, dim=%s, data_count=%s)' % \
(self.__class__, self.a, self.dim, len(self.phase_maps))
return 'DataCollection(%s, a=%s, dim_uv=%s, data_count=%s)' % \
(self.__class__, self.a, self.dim_uv, len(self.data))
def append(self, (phase_map, projector)):
'''Appends a data pair of phase map and projection infos to the data collection.`
......@@ -94,6 +97,6 @@ class DataCollection(object):
self.log.info('Calling append')
assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
'Argument has to be a tuple of a PhaseMap and a Projector object!'
assert phase_map.dim == self.dim, 'Added phasemap must have the same dimension!'
assert (projector.dim_v, projector.dim_u) == self.dim, 'Projector dimensions must match!'
assert phase_map.dim_uv == self.dim_uv, 'Added phasemap must have the same dimension!'
assert projector.dim_uv == self.dim_uv, 'Projector dimensions must match!'
self.data.append((phase_map, projector))
pyramid/forwardmodel.py
+ 23
− 8
View file @ ec080940
......@@ -17,7 +17,7 @@ class ForwardModel:
@property
def projectors(self):
return self._projectors
@projectors.setter
def projectors(self, projectors):
assert np.all([isinstance(projector, Projector) for projector in projectors]), \
......@@ -33,27 +33,42 @@ class ForwardModel:
assert isinstance(kernel, Kernel), 'A Kernel object has to be provided!'
self._kernel = kernel
def __init__(self,projectors, kernel):
def __init__(self, projectors, kernel, b_0):
# TODO: Docstring!
self.kernel = kernel
self.b_0 = kernel.b_0
self.b_0 = b_0
self.a = kernel.a
self.dim = kernel.dim
self.dim_uv = kernel.dim_uv
self.projectors = projectors
def __call__(self, x):
# TODO: Docstring!
# print 'FWD Model - __call__ - input: ', len(x)
result = [self.kernel.jac_dot(projector.jac_dot(x)) for projector in self.projectors]
return np.reshape(result, -1)
result = self.b_0 * np.reshape(result, -1)
# print 'FWD Model - __call__ - output:', len(result)
return result
# TODO: jac_dot ausschreiben!!!
def jac_dot(self, x, vector):
# TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
return self(vector) # The jacobi-matrix does not depend on x in a linear problem
# print 'FWD Model - jac_dot - input: ', len(vector)
result = self(vector)
# print 'FWD Model - jac_dot - output:', len(result)
return result # The jacobi-matrix does not depend on x in a linear problem
def jac_T_dot(self, x, vector):
# TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
# The jacobi-matrix does not depend on x in a linear problem
result = [projector.jac_T_dot(self.kernel.jac_T_dot(x)) for projector in self.projectors]
return np.reshape(result, -1)
# print 'FWD Model - jac_T_dot - input: ', len(vector)
size_3d = self.projectors[0].size_3d
size_2d = np.prod(self.dim_uv)
result = np.zeros(3*size_3d)
for (i, projector) in enumerate(self.projectors):
result += projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
# result = [projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
# for (i, projector) in enumerate(self.projectors)]
result = np.reshape(result, -1)
# print 'FWD Model - jac_T_dot - output:', len(result)
return result
pyramid/kernel.py
+ 22
− 24
View file @ ec080940
......@@ -51,14 +51,12 @@ class Kernel(object):
Attributes
----------
dim : tuple (N=2)
dim_uv : tuple (N=2)
Dimensions of the projected magnetization grid.
a : float
The grid spacing in nm.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
b_0 : float, optional
The saturation magnetic induction. Default is 1.
u : :class:`~numpy.ndarray` (N=3)
The phase contribution of one pixel magnetized in u-direction.
v : :class:`~numpy.ndarray` (N=3)
......@@ -69,7 +67,8 @@ class Kernel(object):
The real FFT of the phase contribution of one pixel magnetized in v-direction.
dim_fft : tuple (N=2)
Dimensions of the grid, which is used for the FFT. Calculated by adding the dimensions
`dim` of the magnetization grid and the dimensions of the kernel (given by ``2*dim-1``)
`dim_uv` of the magnetization grid and the dimensions of the kernel (given by
``2*dim_uv-1``)
and increasing to the next multiple of 2 (for faster FFT).
slice_fft : tuple (N=2) of :class:`slice`
A tuple of :class:`slice` objects to extract the original field of view from the increased
......@@ -77,24 +76,21 @@ class Kernel(object):
'''# TODO: Can be used for several PhaseMappers via the fft arguments or via calling!
def __init__(self, a, dim, b_0=1, numcore=True, geometry='disc'):
def __init__(self, a, dim_uv, numcore=True, geometry='disc'):
'''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
Parameters
----------
a : float
The grid spacing in nm.
dim : tuple (N=2)
dim_uv : tuple (N=2)
Dimensions of the projected magnetization grid.
b_0 : float, optional
The saturation magnetic induction. Default is 1.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
''' # TODO: Docstring
self.log = logging.getLogger(__name__)
self.log.info('Calling __init__')
# TODO: Check if b_0 has an influence or is forgotten
# Function for the phase of an elementary geometry:
def get_elementary_phase(geometry, n, m, a):
if geometry == 'disc':
......@@ -108,23 +104,22 @@ class Kernel(object):
return 0.5 * (F_a(n-0.5, m-0.5) - F_a(n+0.5, m-0.5)
- F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
# Set basic properties:
self.dim = dim # !!! size of the FOV, not the kernel (kernel is bigger)!
self.dim_uv = dim_uv # !!! size of the FOV, not the kernel (kernel is bigger)!
self.a = a
self.numcore = numcore
self.geometry = geometry
self.b_0 = b_0
# Calculate kernel (single pixel phase):
coeff = -a**2 * b_0 / (2*PHI_0)
v_dim, u_dim = dim
coeff = -a**2 / (2*PHI_0)
v_dim, u_dim = dim_uv
u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
uu, vv = np.meshgrid(u, v)
self.u = coeff * get_elementary_phase(geometry, uu, vv, a)
self.v = coeff * get_elementary_phase(geometry, vv, uu, a)
# Calculate Fourier trafo of kernel components:
dim_combined = 3*np.array(dim) - 1 # dim + (2*dim - 1) magnetisation + kernel
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
self.slice_fft = (slice(dim[0]-1, 2*dim[0]-1), slice(dim[1]-1, 2*dim[1]-1))
self.slice_fft = (slice(dim_uv[0]-1, 2*dim_uv[0]-1), slice(dim_uv[1]-1, 2*dim_uv[1]-1))
self.u_fft = np.fft.rfftn(self.u, self.dim_fft)
self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
self.log.info('Created '+str(self))
......@@ -132,6 +127,7 @@ class Kernel(object):
def __call__(self, x):
'''Test'''
self.log.info('Calling __call__')
# print 'Kernel - __call__:', len(x)
if self.numcore:
return self._multiply_jacobi_core(x)
else:
......@@ -140,12 +136,13 @@ class Kernel(object):
def __repr__(self):
self.log.info('Calling __repr__')
return '%s(a=%r, dim=%r, b_0=%r, numcore=%r, geometry=%r)' % \
(self.__class__, self.a, self.dim, self.b_0, self.numcore, self.geometry)
return '%s(a=%r, dim_uv=%r, numcore=%r, geometry=%r)' % \
(self.__class__, self.a, self.dim_uv, self.numcore, self.geometry)
def jac_dot(self, vector):
'''TEST'''# TODO: Docstring
self.log.info('Calling jac_dot')
# print 'Kernel - jac_dot:', len(vector)
if self.numcore:
return self._multiply_jacobi_core(vector)
else:
......@@ -154,6 +151,7 @@ class Kernel(object):
def jac_T_dot(self, vector):
# TODO: Docstring
self.log.info('Calling jac_dot_T')
# print 'Kernel - jac_T_dot:', len(vector)
return self._multiply_jacobi_T(vector)
def _multiply_jacobi(self, vector):
......@@ -173,8 +171,8 @@ class Kernel(object):
'''# TODO: move!
self.log.info('Calling _multiply_jacobi')
v_dim, u_dim = self.dim
size = v_dim * u_dim
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)
......@@ -206,9 +204,9 @@ class Kernel(object):
'''# TODO: move!
self.log.info('Calling _multiply_jacobi_T')
v_dim, u_dim = self.dim
size = v_dim * u_dim
assert len(vector) == size, 'vector size not compatible!'
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
......@@ -223,6 +221,6 @@ class Kernel(object):
def _multiply_jacobi_core(self, vector):
self.log.info('Calling _multiply_jacobi_core')
result = np.zeros(np.prod(self.dim))
core.multiply_jacobi_core(self.dim[0], self.dim[1], self.u, self.v, vector, result)
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
pyramid/logfile.log 0 → 100644
+ 0
− 0
View file @ ec080940
pyramid/optimizer.py
+ 9
− 6
View file @ ec080940
......@@ -34,21 +34,24 @@ from pyramid.magdata import MagData
def optimize_cg(self, data_collection, first_guess):
def optimize_cg(data_collection, first_guess):
# TODO: where should data_collection and first_guess go? here or __init__?
data = data_collection
data = data_collection # TODO: name the parameter 'data' directly...
mag_0 = first_guess
x_0 = first_guess.mag_vec
y = data.phase_vec
kern = Kernel(data.dim, data.a, data.b_0)
F = ForwardModel(data, kern)
kern = Kernel(data.a, data.dim_uv)
F = ForwardModel(data.projectors, kern, data.b_0)
C = Costfunction(y, F)
result = minimize(C, x_0, method='CG', jac=C.jac, hessp=C.hess_dot)
result = minimize(C, x_0, method='Newton-CG', jac=C.jac, hessp=C.hess_dot,
options={'maxiter':200, 'disp':True})
x_opt = result.x
mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim)).mag_vec = x_opt
mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
mag_opt.mag_vec = x_opt
return mag_opt
......
pyramid/phasemap.py
+ 20
− 23
View file @ ec080940
......@@ -38,7 +38,7 @@ class PhaseMap(object):
----------
a: float
The grid spacing in nm.
dim: tuple (N=2)
dim_uv: tuple (N=2)
Dimensions of the grid.
phase: :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
......@@ -107,8 +107,8 @@ class PhaseMap(object):
self._a = a
@property
def dim(self):
return self._dim
def dim_uv(self):
return self._dim_uv
@property
def phase(self):
......@@ -119,7 +119,7 @@ class PhaseMap(object):
assert isinstance(phase, np.ndarray), 'Phase has to be a numpy array!'
assert len(phase.shape) == 2, 'Phase has to be 2-dimensional!'
self._phase = phase
self._dim = phase.shape
self._dim_uv = phase.shape
@property
def phase_vec(self):
......@@ -128,8 +128,8 @@ class PhaseMap(object):
@phase_vec.setter
def phase_vec(self, phase_vec):
assert isinstance(phase_vec, np.ndarray), 'Vector has to be a numpy array!'
assert np.size(phase_vec) == np.prod(self.dim), 'Vector size has to match phase!'
self.phase = phase_vec.reshape(self.dim)
assert np.size(phase_vec) == np.prod(self.dim_uv), 'Vector size has to match phase!'
self.phase = phase_vec.reshape(self.dim_uv)
@property
def unit(self):
......@@ -154,7 +154,7 @@ class PhaseMap(object):
def __str__(self):
self.log.info('Calling __str__')
return 'PhaseMap(a=%s, dim=%s)' % (self.a, self.dim)
return 'PhaseMap(a=%s, dim_uv=%s)' % (self.a, self.dim_uv)
def __neg__(self): # -self
self.log.info('Calling __neg__')
......@@ -167,7 +167,7 @@ class PhaseMap(object):
if isinstance(other, PhaseMap):
self.log.info('Adding two PhaseMap objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.phase.shape == self.dim, \
assert other.phase.shape == self.dim_uv, \
'Added magnitude has to have the same dimensions!'
return PhaseMap(self.a, self.phase+other.phase, self.unit)
else: # other is a Number
......@@ -266,8 +266,8 @@ class PhaseMap(object):
self.log.info('Calling save_to_netcdf4')
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
phase_file.a = self.a
phase_file.createDimension('v_dim', self.dim[0])
phase_file.createDimension('u_dim', self.dim[1])
phase_file.createDimension('v_dim', self.dim_uv[0])
phase_file.createDimension('u_dim', self.dim_uv[1])
phase = phase_file.createVariable('phase', 'f', ('v_dim', 'u_dim'))
phase[:] = self.phase
phase_file.close()
......@@ -334,8 +334,8 @@ class PhaseMap(object):
# Plot the phasemap:
im = axis.pcolormesh(phase, cmap=cmap, vmin=-limit, vmax=limit, norm=norm)
# Set the axes ticks and labels:
axis.set_xlim(0, self.dim[1])
axis.set_ylim(0, self.dim[0])
axis.set_xlim(0, self.dim_uv[1])
axis.set_ylim(0, self.dim_uv[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
......@@ -381,8 +381,8 @@ class PhaseMap(object):
fig = plt.figure()
axis = Axes3D(fig)
# Plot surface and contours:
u = range(self.dim[1])
v = range(self.dim[0])
u = range(self.dim_uv[1])
v = range(self.dim_uv[0])
uu, vv = np.meshgrid(u, v)
axis.plot_surface(uu, vv, phase, rstride=4, cstride=4, alpha=0.7, cmap=cmap,
linewidth=0, antialiased=False)
......@@ -458,8 +458,8 @@ class PhaseMap(object):
axis.set_title(title, fontsize=18)
axis.set_xlabel('x-axis [px]', fontsize=15)
axis.set_ylabel('y-axis [px]', fontsize=15)
axis.set_xlim(0, self.dim[1])
axis.set_ylim(0, self.dim[0])
axis.set_xlim(0, self.dim_uv[1])
axis.set_ylim(0, self.dim_uv[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
# Show Plot:
......@@ -467,9 +467,9 @@ class PhaseMap(object):
plt.show()
# Return plotting axis:
return axis
def display_combined(self, density=1, title='Combined Plot', interpolation='none',
grad_encode='dark'):
def display_combined(self, title='Combined Plot', cmap='RdBu', limit=None, norm=None,
density=1, interpolation='none', grad_encode='dark'):
'''Display the phase map and the resulting color coded holography image in one plot.
Parameters
......@@ -499,7 +499,7 @@ class PhaseMap(object):
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
self.display_phase(axis=phase_axis, show=False)
self.display_phase(cmap='RdBu', limit=limit, norm=norm, axis=phase_axis, show=False)
plt.show()
# Return the plotting axes:
return phase_axis, holo_axis
......@@ -537,6 +537,3 @@ class PhaseMap(object):
axis.xaxis.set_major_locator(NullLocator())
axis.yaxis.set_major_locator(NullLocator())
plt.show()
PhaseMap.make_color_wheel()
pyramid/phasemapper.py
+ 29
− 23
View file @ ec080940
......@@ -46,13 +46,13 @@ class PMAdapterFM(PhaseMapper):
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_2d, b_0, geometry))
self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, geometry), b_0)
def __call__(self, mag_data):
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
# TODO: test if mag_data fits in all aspects
phase_map = PhaseMap(self.a, np.zeros(self.projector.dim_2d))
phase_map = PhaseMap(self.a, np.zeros(self.projector.dim_uv))
phase_map.phase_vec = self.fwd_model(mag_data.mag_vec)
return phase_map
......@@ -94,8 +94,8 @@ class PMFourier(PhaseMapper):
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
# TODO: test if mag_data fits in all aspects (especially with projector)
v_dim, u_dim = self.projector.dim_2d
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_2d)
v_dim, u_dim = self.projector.dim_uv
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Create zero padded matrices:
u_pad = u_dim/2 * self.padding
v_pad = v_dim/2 * self.padding
......@@ -183,7 +183,7 @@ class PMElectric(PhaseMapper):
# Calculate mask:
mask = np.sqrt(np.sum(np.array(mag_data.magnitude)**2, axis=0)) > self.threshold
# Project and calculate phase:
projection = self.projector(mask.reshape(-1)).reshape(self.projector.dim)
projection = self.projector(mask.reshape(-1)).reshape(self.projector.dim_uv)
phase = v_0 * Ce * projection * self.a*1E-9
return PhaseMap(self.a, phase)
......@@ -200,14 +200,15 @@ class PMConvolve(PhaseMapper):
a : float
The grid spacing in nm.
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 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.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
kernel : :class:`~pyramid.kernel.Kernel`, optional
Specifies the kernel for the convolution with the magnetization data. If none is specified,
one will be created with `disc` as the default geometry.
Specifies the kernel for the convolution with the magnetization data. If none is
specified, one will be created with `disc` as the default geometry.
Returns
-------
......@@ -218,13 +219,14 @@ class PMConvolve(PhaseMapper):
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.b_0 = b_0
self.threshold = threshold
self.kernel = Kernel(a, projector.dim_2d, b_0, geometry)
self.kernel = Kernel(a, projector.dim_uv, geometry)
def __call__(self, mag_data):
# Docstring!
# Process input parameters:
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_2d)
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Fourier transform the projected magnetisation:
kernel = self.kernel
u_mag_fft = np.fft.rfftn(u_mag, kernel.dim_fft)
......@@ -233,7 +235,7 @@ class PMConvolve(PhaseMapper):
u_phase = np.fft.irfftn(u_mag_fft * kernel.u_fft, kernel.dim_fft)[kernel.slice_fft].copy()
v_phase = np.fft.irfftn(v_mag_fft * kernel.v_fft, kernel.dim_fft)[kernel.slice_fft].copy()
# Return the result:
return PhaseMap(self.a, u_phase - v_phase)
return PhaseMap(self.a, self.b_0*(u_phase-v_phase))
class PMReal(PhaseMapper):
......@@ -267,30 +269,34 @@ class PMReal(PhaseMapper):
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.b_0 = b_0
self.threshold = threshold
self.kernel = Kernel(a, projector.dim_2d, b_0, geometry)
self.kernel = Kernel(a, projector.dim_uv, geometry)
self.numcore = numcore
def __call__(self, mag_data):
# TODO: Docstring
# Process input parameters:
dim = self.projector.dim_2d
dim_uv = self.projector.dim_uv
threshold = self.threshold
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+dim)
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+dim_uv)
# Create kernel (lookup-tables for the phase of one pixel):
u_phi = self.kernel.u
v_phi = self.kernel.v
u_phi = self.b_0 * self.kernel.u
v_phi = self.b_0 * self.kernel.v
# Calculation of the phase:
phase = np.zeros(dim)
phase = np.zeros(dim_uv)
if self.numcore:
nc.phase_mag_real_core(dim[0], dim[1], v_phi, u_phi, v_mag, u_mag, phase, threshold)
nc.phase_mag_real_core(dim_uv[0], dim_uv[1], v_phi, u_phi,
v_mag, u_mag, phase, threshold)
else:
for j in range(dim[0]):
for i in range(dim[1]):
u_phase = u_phi[dim[0]-1-j:(2*dim[0]-1)-j, dim[1]-1-i:(2*dim[1]-1)-i]
for j in range(dim_uv[0]):
for i in range(dim_uv[1]):
u_phase = u_phi[dim_uv[0]-1-j:(2*dim_uv[0]-1)-j,
dim_uv[1]-1-i:(2*dim_uv[1]-1)-i]
if abs(u_mag[j, i]) > threshold:
phase += u_mag[j, i] * u_phase
v_phase = v_phi[dim[0]-1-j:(2*dim[0]-1)-j, dim[1]-1-i:(2*dim[1]-1)-i]
v_phase = v_phi[dim_uv[0]-1-j:(2*dim_uv[0]-1)-j,
dim_uv[1]-1-i:(2*dim_uv[1]-1)-i]
if abs(v_mag[j, i]) > threshold:
phase -= v_mag[j, i] * v_phase
# Return the phase:
......
pyramid/projector.py
+ 57
− 18
View file @ ec080940
......@@ -30,7 +30,7 @@ class Projector(object):
Attributes
----------
dim_2d : tuple (N=2)
dim_uv : tuple (N=2)
Dimensions (v, u) of the projected grid.
size_3d : int
Number of voxels of the 3-dimensional grid.
......@@ -48,26 +48,19 @@ class Projector(object):
__metaclass__ = abc.ABCMeta
@abc.abstractmethod
def __init__(self, (dim_v, dim_u), weight, coeff):
def __init__(self, dim_uv, weight, coeff):
self.log = logging.getLogger(__name__)
self.log.info('Calling __init__')
self.dim_2d = (dim_v, dim_u)
self.dim_uv = dim_uv
self.weight = weight
self.coeff = coeff
self.size_2d, self.size_3d = weight.shape
self.log.info('Created '+str(self))
def __call__(self, vector):
self.log.info('Calling as function')
if len(vector) == 3*self.size_3d: # mode == 'vector'
self.log.info('mode == vector')
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
self.log.info('mode == scalar')
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for ' \
'vector- or scalar-field-projection!')
self.log.info('Calling as function')
# print 'Projector - __call__:', len(vector)
return self.jac_dot(vector)
def _vector_field_projection(self, vector):
self.log.info('Calling _vector_field_projection')
......@@ -79,19 +72,42 @@ class Projector(object):
if self.coeff[0][1] != 0: # y to u
result[:size_2d] += self.coeff[0][1] * self.weight.dot(vector[size_3d:2*size_3d])
if self.coeff[0][2] != 0: # z to u
result[:size_2d] += self.coeff[0][2] * self.weight.dot(vector[2*size_3d])
result[:size_2d] += self.coeff[0][2] * self.weight.dot(vector[2*size_3d:])
if self.coeff[1][0] != 0: # x to v
result[size_2d:] += self.coeff[1][0] * self.weight.dot(vector[:size_3d])
if self.coeff[1][1] != 0: # y to v
result[size_2d:] += self.coeff[1][1] * self.weight.dot(vector[size_3d:2*size_3d])
if self.coeff[1][2] != 0: # z to v
result[size_2d:] += self.coeff[1][2] * self.weight.dot(vector[2*size_3d])
result[size_2d:] += self.coeff[1][2] * self.weight.dot(vector[2*size_3d:])
return result
def _vector_field_projection_T(self, vector):
self.log.info('Calling _vector_field_projection_T')
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(3*size_3d)
# Go over all possible component projections (z, y, x) to (u, v):
if self.coeff[0][0] != 0: # x to u
result[:size_3d] += self.coeff[0][0] * self.weight.T.dot(vector[:size_2d])
if self.coeff[0][1] != 0: # y to u
result[:size_3d] += self.coeff[0][1] * self.weight.T.dot(vector[size_2d:])
if self.coeff[0][2] != 0: # z to u
result[size_3d:2*size_3d] += self.coeff[0][2] * self.weight.T.dot(vector[:size_2d])
if self.coeff[1][0] != 0: # x to v
result[size_3d:2*size_3d] += self.coeff[1][0] * self.weight.T.dot(vector[size_2d:])
if self.coeff[1][1] != 0: # y to v
result[2*size_3d:] += self.coeff[1][1] * self.weight.T.dot(vector[:size_2d])
if self.coeff[1][2] != 0: # z to v
result[2*size_3d:] += self.coeff[1][2] * self.weight.T.dot(vector[size_2d:])
return result
def _scalar_field_projection(self, vector):
self.log.info('Calling _scalar_field_projection')
return np.array(self.weight.dot(vector))
def _scalar_field_projection_T(self, vector):
self.log.info('Calling _scalar_field_projection_T')
return np.array(self.weight.T.dot(vector))
def jac_dot(self, vector):
'''Multiply a `vector` with the jacobi matrix of this :class:`~.Projector` object.
......@@ -109,7 +125,30 @@ class Projector(object):
'''
self.log.info('Calling jac_dot')
return self(vector)
# print 'Projector - jac_dot:', len(vector)
if len(vector) == 3*self.size_3d: # mode == 'vector'
self.log.info('mode == vector')
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
self.log.info('mode == scalar')
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for ' \
'vector- or scalar-field-projection!')
def jac_T_dot(self, vector):
# TODO: Docstring!
self.log.info('Calling jac_T_dot')
# print 'Projector - jac_T_dot:', len(vector)
if len(vector) == 2*self.size_2d: # mode == 'vector'
self.log.info('mode == vector')
return self._vector_field_projection_T(vector)
elif len(vector) == self.size_2d: # mode == 'scalar'
self.log.info('mode == scalar')
return self._scalar_field_projection_T(vector)
else:
raise AssertionError('Vector size has to be suited either for ' \
'vector- or scalar-field-projection!')
......@@ -124,7 +163,7 @@ class YTiltProjector(Projector):
#
# Attributes
# ----------
# dim_2d : tuple (N=2)
# dim_uv : tuple (N=2)
# Dimensions (v, u) of the projected grid.
# size_3d : int
# Number of voxels of the 3-dimensional grid.
......@@ -239,7 +278,7 @@ class SimpleProjector(Projector):
#
# Attributes
# ----------
# dim_2d : tuple (N=2)
# dim_uv : tuple (N=2)
# Dimensions (v, u) of the projected grid.
# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
# The weight matrix containing the weighting coefficients which determine the influence of
......
scripts/edinburgh_test.py
+ 1
− 1
View file @ ec080940
......@@ -42,7 +42,7 @@ phasemapper = PMAdapterFM(mag_data.a, projector)
phase_map = phasemapper(mag_data)
phase_axis = phase_map.display_combined(density=20, interpolation='bilinear')[0]
phase_axis = phase_map.display_combined(density=20, interpolation='bilinear', grad_encode='bright')[0]
phase_axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:3.0f}'.format(x*mag_data.a)))
phase_axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:3.0f}'.format(x*mag_data.a)))
......
scripts/images_poster.py
+ 70
− 42
View file @ ec080940
......@@ -6,58 +6,86 @@ Created on Wed Feb 05 19:19:19 2014
"""
import os
from numpy import pi
import numpy as np
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.magdata import MagData
from pyramid.phasemapper import PMConvolve
from pyramid.projector import SimpleProjector
from pyramid.phasemap import PhaseMap
import matplotlib.pyplot as plt
from time import clock
directory = '../output/poster'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
a = 10.0 # nm
dim = (64, 128, 128)
# Slab:f
center = (32, 32, 32) # in px (z, y, x), index starts with 0!
width = (322, 48, 48) # in px (z, y, x)
mag_shape_slab = mc.Shapes.slab(dim, center, width)
# Disc:
center = (32, 32, 96) # in px (z, y, x), index starts with 0!
radius = 24 # in px
height = 24 # in px
mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
# Sphere:
center = (32, 96, 64) # in px (z, y, x), index starts with 0!
radius = 24 # in px
mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
# Create empty MagData object and add magnetized objects:
mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, pi/6))
mag_data += MagData(a, mc.create_mag_dist_vortex(mag_shape_disc, (32, 96)))
mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, -pi/4))
# Plot the magnetic distribution, phase map and holographic contour map:
mag_data_coarse = mag_data.copy()
mag_data_coarse.scale_down(2)
mag_data_coarse.quiver_plot()
plt.savefig(os.path.join(directory, 'mag.png'))
mag_data_coarse.quiver_plot3d()
projector = SimpleProjector(dim)
phase_map = PMConvolve(a, projector, b_0=0.1)(mag_data)
phase_map.display_phase()
plt.savefig(os.path.join(directory, 'phase.png'))
phase_map.display_holo(density=4, interpolation='bilinear')
plt.savefig(os.path.join(directory, 'holo.png'))
#directory = '../output/poster'
#if not os.path.exists(directory):
# os.makedirs(directory)
## Input parameters:
#a = 10.0 # nm
#dim = (64, 128, 128)
## Slab:f
#center = (32, 32, 32) # in px (z, y, x), index starts with 0!
#width = (322, 48, 48) # in px (z, y, x)
#mag_shape_slab = mc.Shapes.slab(dim, center, width)
## Disc:
#center = (32, 32, 96) # in px (z, y, x), index starts with 0!
#radius = 24 # in px
#height = 24 # in px
#mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
## Sphere:
#center = (32, 96, 64) # in px (z, y, x), index starts with 0!
#radius = 24 # in px
#mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
## Create empty MagData object and add magnetized objects:
#mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, pi/6))
#mag_data += MagData(a, mc.create_mag_dist_vortex(mag_shape_disc, (32, 96)))
#mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, -pi/4))
## Plot the magnetic distribution, phase map and holographic contour map:
#mag_data_coarse = mag_data.copy()
#mag_data_coarse.scale_down(2)
#mag_data_coarse.quiver_plot()
#plt.savefig(os.path.join(directory, 'mag.png'))
#mag_data_coarse.quiver_plot3d()
#projector = SimpleProjector(dim)
#phase_map = PMConvolve(a, projector, b_0=0.1)(mag_data)
#phase_map.display_phase()
#plt.savefig(os.path.join(directory, 'phase.png'))
#phase_map.display_holo(density=4, interpolation='bilinear')
#plt.savefig(os.path.join(directory, 'holo.png'))
#
#
#dim = (1, 9, 9)
#mag_shape = mc.Shapes.pixel(dim, (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)))
#mag_data_x = MagData(a, mc.create_mag_dist_homog(mag_shape, 0))
#mag_data_y = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/2))
#phasemapper = PMConvolve(a, SimpleProjector(dim))
#phasemapper(mag_data_x).display_phase()
#phasemapper(mag_data_y).display_phase()
dim = (1, 9, 9)
mag_shape = mc.Shapes.pixel(dim, (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)))
mag_data_x = MagData(a, mc.create_mag_dist_homog(mag_shape, 0))
mag_data_y = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/2))
a = 1.
dim = (128, 128, 128)
center = (dim[0]/2, dim[1]/2, dim[2]/2)
radius = dim[1]/4 # in px
mag_shape = mc.Shapes.sphere(dim, center, radius)
mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/4))
phasemapper = PMConvolve(a, SimpleProjector(dim))
phasemapper(mag_data_x).display_phase()
phasemapper(mag_data_y).display_phase()
\ No newline at end of file
start = clock()
phase_map = phasemapper(mag_data)
print 'TIME:', clock() - start
phase_ana = an.phase_mag_sphere(dim, a, pi/4, center, radius)
phase_diff = (phase_ana - phase_map).phase
print 'RMS%:', np.sqrt(np.mean((phase_diff)**2))/phase_ana.phase.max()*100
print 'max:', phase_ana.phase.max()
print 'RMS:', np.sqrt(np.mean((phase_diff)**2))
PhaseMap(a, phase_diff).display_phase()
mag_data.scale_down(2)
mag_data.quiver_plot()
mag_data.quiver_plot3d()
phase_map.display_phase()
scripts/logfile.log
+ 0
− 1
View file @ ec080940
2014-02-08 21:43:12: INFO @ <root>: Calling copy
scripts/paper 1/ch5-4-comparison_of_methods_new.py
+ 3
− 1
View file @ ec080940
......@@ -25,7 +25,7 @@ from pyramid.phasemapper import PMConvolve, PMFourier
import shelve
force_calculation = False
force_calculation = True
print '\nACCESS SHELVE'
......@@ -133,6 +133,8 @@ else:
print '------time (disc, real disc) =', data_disc_real_d[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
data_disc_real_d[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
print 'TIME:', data_disc_real_d[2, i]
print 'RMS%:', np.sqrt(np.mean(((phase_ana_disc-phase_num_disc).phase/phase_ana_disc.phase)**2))*100, '%'
print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
# Analytic solution:
......
scripts/paper 1/logfile.log
+ 84
− 85
View file @ ec080940
2014-02-08 20:46:32: INFO @ <pyramid.phasemap>: Calling make_color_wheel
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling __init__
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling __init__
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0A8519F0>
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0A8519F0>
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling as function
2014-02-08 20:46:43: INFO @ <pyramid.projector>: mode == vector
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling as function
2014-02-08 20:46:43: INFO @ <pyramid.projector>: mode == vector
2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __isub__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Adding an offset
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-08 20:46:44: INFO @ <pyramid.projector>: Calling as function
2014-02-08 20:46:44: INFO @ <pyramid.projector>: mode == vector
2014-02-08 20:46:44: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-08 20:46:45: INFO @ <pyramid.projector>: Calling as function
2014-02-08 20:46:45: INFO @ <pyramid.projector>: mode == vector
2014-02-08 20:46:45: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __isub__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Adding an offset
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-08 20:46:47: INFO @ <pyramid.projector>: Calling __init__
2014-02-08 20:46:47: INFO @ <pyramid.projector>: Calling __init__
2014-02-08 20:46:47: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0B117190>
2014-02-08 20:46:47: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0B117190>
2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __isub__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding an offset
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __neg__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __isub__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding an offset
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling display_phase
2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
scripts/simple_phasemapping.py
+ 1
− 1
View file @ ec080940
......@@ -15,7 +15,7 @@ from time import clock
#import cProfile
mag_data = MagData.load_from_netcdf4('../output/vtk data/rect_500x125x3.nc')
mag_data = MagData.load_from_netcdf4('../output/vtk data/tube_90x30x30.nc')
projector = SimpleProjector(mag_data.dim)
......
scripts/simple_reconstruction.py 0 → 100644
+ 171
− 0
View file @ ec080940
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 28 14:25:59 2014
@author: Jan
"""
import numpy as np
from numpy import pi
from pyramid.magdata import MagData
from pyramid.projector import YTiltProjector
from pyramid.phasemapper import PMConvolve
from pyramid.datacollection import DataCollection
import pyramid.optimizer as opt
from pyramid.kernel import Kernel
from pyramid.forwardmodel import ForwardModel
from pyramid.costfunction import Costfunction
a = 1.
b_0 = 1000.
dim = (3, 3, 3)
count = 1
###################################################################################################
print('--Generating input phase_maps')
magnitude = np.zeros((3,)+dim)
magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
mag_data = MagData(a, magnitude)
mag_data.quiver_plot3d()
tilts = np.linspace(0, 2*pi, num=count, endpoint=False)
projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
phase_maps = [pm(mag_data) for pm in phasemappers]
#[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
# for i, phase_map in enumerate(phase_maps)]
phase_maps[0].display_phase()
###################################################################################################
print('--Setting up data collection')
dim_uv = dim[1:3]
data_collection = DataCollection(a, dim_uv, b_0)
[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
###################################################################################################
print('--Test optimizer')
first_guess = MagData(a, np.zeros((3,)+dim))
first_guess.magnitude[1, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1
#first_guess.magnitude[0, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = -1
first_guess.quiver_plot3d()
phase_guess = PMConvolve(first_guess.a, projectors[0], b_0)(first_guess)
phase_guess.display_phase()
#mag_opt = opt.optimize_cg(data_collection, first_guess)
#
#mag_opt.quiver_plot3d()
#
#phase_opt = PMConvolve(mag_opt.a, projectors[0], b_0)(mag_opt)
#
#phase_opt.display_phase()
###################################################################################################
print('--Further testing')
data = data_collection
mag_0 = first_guess
x_0 = first_guess.mag_vec
y = data.phase_vec
kern = Kernel(data.a, data.dim_uv)
F = ForwardModel(data.projectors, kern, data.b_0)
C = Costfunction(y, F)
size_3d = np.prod(dim)
cost = C(first_guess.mag_vec)
cost_grad = C.jac(first_guess.mag_vec)
cost_grad_del_x = cost_grad.reshape((3, 3, 3, 3))[0, ...]
cost_grad_del_y = cost_grad.reshape((3, 3, 3, 3))[1, ...]
cost_grad_del_z = cost_grad.reshape((3, 3, 3, 3))[2, ...]
x_t = np.asmatrix(mag_data.mag_vec).T
y = np.asmatrix(y).T
K = np.array([F.jac_dot(x_0, np.eye(81)[:, i]) for i in range(81)]).T
K = np.asmatrix(K)
lam = 10. ** -10
KTK = K.T * K + lam * np.asmatrix(np.eye(81))
print lam,
#print pylab.cond(KTK),
x_f = KTK.I * K.T * y
print (np.asarray(K * x_f - y) ** 2).sum(),
print (np.asarray(K * x_t - y) ** 2).sum()
print x_f
x_rec = np.asarray(x_f).reshape(3,3,3,3)
print x_rec[0, ...]
#KTK = K.T * K
#u,s,v = pylab.svd(KTK, full_matrices=1)
#si = np.zeros_like(s)
#
#si[s>lam] = 1. / s[s>lam]
#si=np.asmatrix(np.diag(si))
#KTKI = np.asmatrix(v).T * si * np.asmatrix(u)
#
#
#x_f = KTKI * K.T * y
#print x_f
#
#
###################################################################################################
print('--Compliance test')
# test F
# test F.jac
# test F.jac_dot
#F_evaluate = F(mag_data.mag_vec)
#F_jacobi = F.jac_dot(None, mag_data.mag_vec)
#np.testing.assert_equal(F_evaluate, F_jacobi)
#
#F_reverse = F.jac_T_dot(None, F_jacobi)
#np.testing.assert_equal(F_reverse, mag_data.mag_vec)
from scipy.sparse.linalg import cg
K = np.asmatrix([F.jac_dot(x_0, np.eye(81)[:, i]) for i in range(81)]).T
lam = 10. ** -10
KTK = K.T * K + lam * np.asmatrix(np.eye(81))
A = KTK#np.array([F(np.eye(81)[:, i]) for i in range(81)])
b = F.jac_T_dot(None, y)
b_test = np.asarray((K.T * y).T)[0]
x_f = cg(A, b_test)[0]
mag_data_rec = MagData(a, np.zeros((3,)+dim))
mag_data_rec.mag_vec = x_f
mag_data_rec.quiver_plot3d()
scripts/test_projection.py 0 → 100644
+ 68
− 0
View file @ ec080940
# -*- coding: utf-8 -*-
"""
Created on Fri Jan 17 14:06:01 2014
@author: Jan
"""
import numpy as np
from numpy import pi
from pyramid.magdata import MagData
from pyramid.projector import YTiltProjector
from pyramid.phasemapper import PMConvolve
from time import clock
import matplotlib.pyplot as plt
#data = np.loadtxt('../output/data from Edinburgh/long_grain_remapped_0p0035.txt', delimiter=',')
#
#a = 1000 * (data[1, 2] - data[0, 2])
#
#dim = len(np.unique(data[:, 2])), len(np.unique(data[:, 1])), len(np.unique(data[:, 0]))
#
#mag_vec = np.concatenate([data[:, 3], data[:, 4], data[:, 5]])
#
#x_mag = np.reshape(data[:, 3], dim, order='F')
#y_mag = np.reshape(data[:, 4], dim, order='F')
#z_mag = np.reshape(data[:, 5], dim, order='F')
#
#magnitude = np.array((x_mag, y_mag, z_mag))
#
#mag_data = MagData(a, magnitude)
#
#mag_data.pad(30, 20, 0)
PATH = '../output/vtk data/tube_90x30x50_sat_edge_equil.gmr'
mag_data = MagData.load_from_netcdf4(PATH+'.nc')
#for tilt in np.linspace(0, 2*pi, num=20, endpoint=False):
# projector = YTiltProjector(mag_data.dim, tilt)
# phasemapper = PMConvolve(mag_data.a, projector)
# phase_map = phasemapper(mag_data)
# (-phase_map).display_combined(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilt/pi),
# limit=2., density=12, interpolation='bilinear',
# grad_encode='bright')
# plt.savefig(PATH+'_tilt_{:3.2f}pi.png'.format(tilt/pi))
mag_data.scale_down()
mag_data.quiver_plot3d()
#projectors = [YTiltProjector(mag_data.dim, i)
# for i in np.linspace(0, 2*pi, num=20, endpoint=False)]
#
#start = clock()
#phasemapper = PMConvolve[PMConvolve(mag_data.a, projector) for projector in projectors]
#print 'Overhead :', clock()-start
#
#start = clock()
#phase_maps = [pm(mag_data) for pm in phasemappers]
#print 'Phasemapping:', clock()-start
#
#[phase_map.display_combined(density=12, interpolation='bilinear', grad_encode='bright')
# for phase_map in phase_maps]
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment