From 275ecaaf9fb133c1efef3cafa600dd224c41c2cf Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Mon, 18 Aug 2014 00:36:10 +0200
Subject: [PATCH] further work on projectors, preparation for regularisation
class (next task!) projector: ALL projectors can now handle arbitrary
dimensions dataset: reverted to old strategy (dimensions must match)
reconstruction: least square method now also with regularization of order 1
scripts: cleanup and streamlining, added projector_test
---
pyramid/dataset.py | 9 +-
pyramid/projector.py | 96 ++--------
pyramid/reconstruction.py | 22 ++-
scripts/collaborations/rueffner_file.py | 97 +++-------
.../vtk_conversion_nanowire_full.py | 178 +++++-------------
scripts/collaborations/zi_an.py | 97 +++++-----
scripts/create distributions/create_sample.py | 4 +-
scripts/test methods/compare_pixel_fields.py | 43 ++---
scripts/test methods/projector_test.py | 29 +++
9 files changed, 213 insertions(+), 362 deletions(-)
create mode 100644 scripts/test methods/projector_test.py
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index bc91332..f51ddfc 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -48,7 +48,7 @@ class DataSet(object):
@property
def dim_uv(self):
- return self._dim_uv # TODO: delete
+ return self._dim_uv
@property
def phase_vec(self):
@@ -99,11 +99,8 @@ class DataSet(object):
self.LOG.debug('Calling append')
assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
'Argument has to be a tuple of a PhaseMap and a Projector object!'
-# assert 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!'
- # TODO: delete
- assert phase_map.dim_uv == projector.dim_uv, \
- 'PhaseMap and 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))
def display_phase(self, phase_maps=None, title='Phase Map', cmap='RdBu', limit=None, norm=None):
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 254a2c0..56d0a75 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -11,7 +11,6 @@ import abc
import itertools
from scipy.sparse import coo_matrix, csr_matrix
-import scipy.sparse as sp
import logging
@@ -246,6 +245,7 @@ class XTiltProjector(Projector):
if dim_uv is None:
dim_uv = (max(dim_perp, dim_proj), dim_rot) # x-y-plane
dim_v, dim_u = dim_uv # y, x
+ assert dim_v >= dim_perp and dim_u >= dim_rot, 'Projected dimensions are too small!'
size_2d = np.prod(dim_uv)
size_3d = np.prod(dim)
# Creating coordinate list of all voxels:
@@ -352,6 +352,7 @@ class YTiltProjector(Projector):
if dim_uv is None:
dim_uv = (dim_rot, max(dim_perp, dim_proj)) # x-y-plane
dim_v, dim_u = dim_uv # y, x
+ assert dim_v >= dim_rot and dim_u >= dim_perp, 'Projected dimensions are too small!'
size_2d = np.prod(dim_uv)
size_3d = np.prod(dim)
# Creating coordinate list of all voxels:
@@ -420,11 +421,13 @@ class SimpleProjector(Projector):
axis : {'z', 'y', 'x'}, optional
Main axis along which the magnetic distribution is projected (given as a string). Defaults
to the z-axis.
+ dim_uv : tuple (N=2), optional
+ Dimensions (v, u) of the projection. If not set it uses the 3D default dimensions.
'''
LOG = logging.getLogger(__name__+'.SimpleProjector')
- AXIS_DICT = {'z': (0, 1, 2), 'y': (1, 0, 2), 'x': (1, 2, 0)}
+ AXIS_DICT = {'z': (0, 1, 2), 'y': (1, 0, 2), 'x': (2, 1, 0)} # (0:z, 1:y, 2:x) -> (proj, v, u)
def __init__(self, dim, axis='z', dim_uv=None):
self.LOG.debug('Calling __init__')
@@ -432,12 +435,6 @@ class SimpleProjector(Projector):
proj, v, u = self.AXIS_DICT[axis]
dim_proj, dim_v, dim_u = dim[proj], dim[v], dim[u]
dim_z, dim_y, dim_x = dim
-
-# if dim_uv is None:
-# dim_uv = (dim_rot, max(dim_perp, dim_proj)) # x-y-plane
-# dim_v, dim_u = dim_uv # y, x
-# size_2d = np.prod(dim_uv)
-
size_2d = dim_u * dim_v
size_3d = np.prod(dim)
data = np.repeat(1, size_3d) # size_3d ones in the matrix (each voxel is projected)
@@ -450,87 +447,32 @@ class SimpleProjector(Projector):
elif axis == 'y':
self.LOG.debug('Projection along the y-axis')
coeff = [[1, 0, 0], [0, 0, 1]]
- indices = np.array([np.arange(row%dim_x, dim_x*dim_y, dim_x)
- + int(row/dim_x)*dim_x*dim_y
+ indices = np.array([np.arange(row%dim_x, dim_x*dim_y, dim_x)+int(row/dim_x)*dim_x*dim_y
for row in range(size_2d)]).reshape(-1)
elif axis == 'x':
self.LOG.debug('Projection along the x-axis')
- coeff = [[0, 1, 0], [0, 0, 1]]
- indices = np.array([np.arange(dim_proj) + row*dim_proj
+ coeff = [[0, 0, 1], [0, 1, 0]] # Caution, coordinate switch: u, v --> z, y (not y, z!)
+ # indices = np.array([np.arange(dim_proj) + row*dim_proj
+ # for row in range(size_2d)]).reshape(-1) # this is u, v --> y, z
+ indices = np.array([np.arange(dim_x) + (row%dim_z)*dim_x*dim_y + int(row/dim_z)*dim_x
for row in range(size_2d)]).reshape(-1)
if dim_uv is not None:
- indptr = indptr.tolist()
- d_v, d_u = int((dim_uv[0]-dim_v)/2), int((dim_uv[1]-dim_u)/2)
- print d_v, d_u
- print dim_v*dim_u, '-->', np.prod(dim_uv)
- print 'before:', indptr
- for i in range(d_v*dim_uv[1]): # last v_slices
- indptr.append(indptr[-1]) # append empty rows
- print 'last: ', indptr
+ indptr = indptr.tolist() # convert to use insert() and append()
+ d_v, d_u = int((dim_uv[0]-dim_v)/2), int((dim_uv[1]-dim_u)/2) # padding in u and v
+ indptr[-1:-1] = [indptr[-1]] * d_v*dim_uv[1] # append empty rows at the end
for i in np.arange(dim_v, 0, -1): # all slices in between
u, l = i*dim_u, (i-1)*dim_u+1 # upper / lower slice end
- print i, u, l
- print indptr[u], indptr[l]
indptr[u:u] = [indptr[u]] * d_u # end of the slice
indptr[l:l] = [indptr[l]] * d_u # start of the slice
- print 'middle:', indptr
- for i in range(d_v*dim_uv[1]): # first v_slices
- indptr.insert(0, 0) # append empty rows
+ indptr[0:0] = [0] * d_v*dim_uv[1] # insert empty rows at the beginning
size_2d = np.prod(dim_uv) # increase size_2d
- print 'after: ', indptr
- else:
+ # Make sure dim_uv is defined (used for the assertion)
+ if dim_uv is None:
dim_uv = dim_v, dim_u
-
-# TODO: assertions!
-
-# else: # dimensions of the projection differ from 3D-distribution
-# d_v, d_u = int((dim_uv[0]-dim_v)/2), int((dim_uv[1]-dim_u)/2)
-# filled_rows =
-# # else --> filled_rows = range(size_2d)
-# if axis == 'z':
-# self.LOG.debug('Projecting along the z-axis')
-# coeff = [[1, 0, 0], [0, 1, 0]]
-#
-# for i, row in enumerate()
-# indices = np.array([np.arange(row, size_3d, size_2d)
-# for row in range(size_2d)]).reshape(-1)
-# print indices.shape
-# elif axis == 'y':
-# self.LOG.debug('Projection along the y-axis')
-# coeff = [[1, 0, 0], [0, 0, 1]]
-# indices = np.array([np.arange(row%dim_x, dim_x*dim_y, dim_x)
-# + int(row/dim_x)*dim_x*dim_y
-# for row in range(size_2d)]).reshape(-1)
-# elif axis == 'x':
-# self.LOG.debug('Projection along the x-axis')
-# coeff = [[0, 1, 0], [0, 0, 1]]
-# indices = np.array([np.arange(dim_proj) + row*dim_proj
-# for row in range(size_2d)]).reshape(-1)
-
+ assert dim_uv[0] >= dim_v and dim_uv[1] >= dim_u, 'Projected dimensions are too small!'
+ # Create weight-matrix:
weight = csr_matrix((data, indices, indptr), shape=(size_2d, size_3d))
-# print len(data), data
-# print len(indices), indices
-# print len(indptr), indptr
-# print np.prod(dim_uv)
-# print size_3d
-# print size_2d
-# weight = csr_matrix((data, indices, indptr), shape=(size_2d+2*d_u, size_3d))
-# print dim_uv is not (dim_v, dim_u)
-# print (dim_uv is not (dim_v, dim_u))
-# print dim_uv, dim_v, dim_u
-# print (dim_uv is not None) and (dim_uv is not (dim_v, dim_u))
-# if (dim_uv is not None) and (dim_uv != (dim_v, dim_u)):
-# d_v, d_u = int((dim_uv[0]-dim_v)/2), int((dim_uv[1]-dim_u)/2)
-# print d_v, d_u
-# print weight.shape
-# print csr_matrix((d_u, size_2d)).shape
-# print csr_matrix((np.prod(dim_uv), d_v)).shape
-# weight = sp.hstack((csr_matrix((size_2d, d_u)), weight,
-# csr_matrix((size_2d, d_u))))
-# weight = sp.vstack((csr_matrix((d_v, np.prod(dim_uv))), weight,
-# csr_matrix((d_v, np.prod(dim_uv)))))
-# print weight.shape
- super(SimpleProjector, self).__init__(dim, dim_uv, weight, coeff) # TODO: dim_uv caution
+ super(SimpleProjector, self).__init__(dim, dim_uv, weight, coeff)
self.LOG.debug('Created '+str(self))
def get_info(self):
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index e61ba1d..e251d28 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -158,7 +158,7 @@ def optimize_cg(data, first_guess):
return mag_opt
-def optimize_simple_leastsq(phase_map, mask, b_0=1):
+def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
'''Reconstruct a magnetic distribution for a 2-D problem with known pixel locations.
Parameters
@@ -172,6 +172,11 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1):
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
+ lam : float, optional
+ The regularisation parameter. Defaults to 1E-4.
+ order : int {0, 1}, optional
+ order of the regularisation function. Default is 0 for a Tikhonov regularisation of order
+ zero. A first order regularisation, which uses the derivative is available with value 1.
Returns
-------
mag_data : :class:`~pyramid.magdata.MagData`
@@ -188,7 +193,6 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1):
a = phase_map.a # Grid spacing
dim = mask.shape # Dimensions of the mag. distr.
count = mask.sum() # Number of pixels with magnetization
- lam = 1e-4 # Regularisation parameter
# Create empty MagData object for the reconstruction:
mag_data_rec = MagData(a, np.zeros((3,)+dim))
@@ -198,16 +202,15 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1):
phase_map = PMConvolve(a, SimpleProjector(dim), b_0)(mag_data_rec)
return phase_map.phase_vec
- # TODO: cleanup
-
- # Cost function which should be minimized:
- def J(x_i):
+ # Cost function of order zero which should be minimized:
+ def J_0(x_i):
y_i = F(x_i)
term1 = (y_i - y_m)
term2 = lam * x_i
return np.concatenate([term1, term2])
- def J2(x_i):
+ # First order cost function which should be minimized:
+ def J_1(x_i):
y_i = F(x_i)
term1 = (y_i - y_m)
mag_data = mag_data_rec.magnitude
@@ -221,7 +224,10 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1):
term2 = lam * np.concatenate(term2)
return np.concatenate([term1, term2])
+ J_DICT = [J_0, J_1] # list of cost-functions with different regularisations
+ print J_DICT
+ print J_DICT[order]
# Reconstruct the magnetization components:
- x_rec, _ = leastsq(J2, np.zeros(3*count))
+ x_rec, _ = leastsq(J_DICT[order], np.zeros(3*count))
mag_data_rec.set_vector(mask, x_rec)
return mag_data_rec
diff --git a/scripts/collaborations/rueffner_file.py b/scripts/collaborations/rueffner_file.py
index 64f0af6..6b37fea 100644
--- a/scripts/collaborations/rueffner_file.py
+++ b/scripts/collaborations/rueffner_file.py
@@ -1,15 +1,10 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Jul 25 14:37:11 2014
-
-@author: Jan
-"""
+"""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
@@ -31,37 +26,26 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
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)
+dim = (16, 190, 220)
+dim_uv = (300, 500)
a = 1.
b_0 = 1.1
-projector_z = SimpleProjector(dim, axis='z')
-projector_x = SimpleProjector(dim, axis='x')
+dens_z = 4E3
+dens_x = 1E2
+lim_z = 22
+lim_x = 0.35
+###################################################################################################
+# Build projectors and phasemapper:
+projector_z = SimpleProjector(dim, axis='z', dim_uv=dim_uv)
+projector_x = SimpleProjector(dim, axis='x', dim_uv=dim_uv)
pm_z = PMConvolve(a, projector_z, b_0)
-pm_x = PMConvolve(a, projector_z, b_0)
-
+pm_x = PMConvolve(a, projector_x, b_0)
+# Read in hdf5-file and extract points:
h5file = tables.openFile(PATH+'dyn_0090mT_dyn.h5_dat.h5')
points = h5file.root.mesh.points.read()
-
+# Plot point distribution in 2D and 3D:
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),
@@ -76,56 +60,37 @@ iz_y = np.concatenate([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)
+xs = np.arange(-dim[2]/2, dim[2]/2)
+ys = np.arange(-dim[1]/2, dim[1]/2)
+zs = np.arange(-dim[0]/2, dim[0]/2)
xx, yy = np.meshgrid(xs, ys)
-
-#test = []
-
-for t in np.arange(0, 1001, 5):
+# Interpolation and phase calculation for all timesteps:
+for t in np.arange(865, 1001, 5):
print 't =', t
vectors = h5file.root.data.fields.m.read(field='m_CoFeb')[t, ...]
data = np.hstack((points, vectors))
-# if data_old is not None:
-# np.testing.assert_equal(data, data_old)
-# data_old = np.copy(data)
-
- zs_unique = np.unique(data[:, 2]) # TODO: not used
-
# Create empty magnitude:
magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
-
+ # Go over all z-slices:
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(data[:, 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[:, 3+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)
-
+ # Create magnetic distribution and phase maps:
mag_data = MagData(a, magnitude)
- mag_data.pad(20, 20, 0)
- phase_map = pm(mag_data)
- phase_map.unit = 'mrad'
- phase_map.display_combined(density=100, interpolation='bilinear', limit=None,
- 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()
+ phase_map_z = pm_z(mag_data)
+ phase_map_x = pm_x(mag_data)
+ phase_map_z.unit = 'mrad'
+ # Display phase maps and save them to png:
+ phase_map_z.display_combined(density=dens_z, interpolation='bilinear', limit=lim_z)
+ plt.savefig(PATH+'rueffner/phase_map_z_t_'+str(t)+'.png')
+ phase_map_x.display_combined(density=dens_x, interpolation='bilinear', limit=lim_x)
+ plt.savefig(PATH+'rueffner/phase_map_x_t_'+str(t)+'.png')
+ # Close all plots to avoid clutter:
+ plt.close('all')
diff --git a/scripts/collaborations/vtk_conversion_nanowire_full.py b/scripts/collaborations/vtk_conversion_nanowire_full.py
index a0ac4ae..a0f52f9 100644
--- a/scripts/collaborations/vtk_conversion_nanowire_full.py
+++ b/scripts/collaborations/vtk_conversion_nanowire_full.py
@@ -9,7 +9,6 @@ from numpy import pi
import matplotlib.pyplot as plt
from matplotlib._pylab_helpers import Gcf
from pylab import griddata
-from mayavi import mlab
import pickle
import vtk
@@ -18,7 +17,8 @@ from tqdm import tqdm
import pyramid
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector, YTiltProjector, XTiltProjector
-from pyramid.phasemapper import PMAdapterFM, PMConvolve
+from pyramid.phasemapper import PMConvolve
+from pyramid.phasemap import PhaseMap
from pyramid.dataset import DataSet
import logging
@@ -32,8 +32,11 @@ logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
PATH = '../../output/vtk data/longtube_withcap/CoFeB_tube_cap_4nm'
b_0 = 1.54
-gain = 1
force_calculation = False
+count = 16
+dim_uv_x = (500, 100)
+dim_uv_y = (100, 500)
+density = 8
###################################################################################################
# Load vtk-data:
if force_calculation or not os.path.exists(PATH+'.pickle'):
@@ -64,18 +67,10 @@ if force_calculation or not os.path.exists(PATH+'.pickle'):
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, 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'):
-
- # Determine the size of object:
+ # Determine the size of the 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()
@@ -106,15 +101,10 @@ if force_calculation or not os.path.exists(PATH+'.nc'):
xx, yy = np.meshgrid(xs, ys)
# Create empty magnitude:
magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
-
+ # Go over all slices:
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]
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!
if z <= 0:
grid_x = np.concatenate([z_slice[:, 0], iz_x])
@@ -126,131 +116,51 @@ if force_calculation or not os.path.exists(PATH+'.nc'):
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),
-# 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, :]
-# 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:
+ # Creating MagData object and save it as netcdf-file:
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(proj_axis='x')
###################################################################################################
-# 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))
-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()
-
-count = 16
-
-dim_uv_x = (500, 100)
-dim_uv_y = (100, 500)
-
-density = 8
-
+mag_data_scaled = mag_data.copy()
+mag_data_scaled.scale_down()
+mag_data_scaled.save_to_netcdf4(PATH+'_scaled.nc')
+# Create tilts and projectors:
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 = [YiltProjector(mag_data.dim, tilt, dim_uv=dim_uv_y) for tilt in tilts_miss]
-projectors_x = [XTiltProjector(mag_data.dim, tilt, dim_uv=dim_uv_x) 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.a, dim_uv_x, b_0)
-
-for i, pm in enumerate(phasemappers):
- data_set.append((pm(mag_data), projectors[i]))
-
-plt.close('all')
-
-data_set.display_combined(density=density, interpolation='bilinear')
-
+projectors_y = [YTiltProjector(mag_data_scaled.dim, tilt, dim_uv=dim_uv_y) for tilt in tilts_miss]
+projectors_x = [XTiltProjector(mag_data_scaled.dim, tilt, dim_uv=dim_uv_x) for tilt in tilts_miss]
+pm_y = [PMConvolve(mag_data_scaled.a, proj, b_0) for proj in projectors_y]
+pm_x = [PMConvolve(mag_data_scaled.a, proj, b_0) for proj in projectors_x]
+# Create data sets for x and y tilts:
+data_set_y = DataSet(mag_data_scaled.a, dim_uv_y, b_0)
+data_set_x = DataSet(mag_data_scaled.a, dim_uv_x, b_0)
+# Create and append phase maps and projectors:
+for i in range(len(pm_y)):
+ data_set_y.append((pm_y[i](mag_data_scaled), projectors_y[i]))
+ data_set_x.append((pm_x[i](mag_data_scaled), projectors_x[i]))
+# Display phase maps:
+data_set_y.display_combined(density=density, interpolation='bilinear')
+data_set_x.display_combined(density=density, interpolation='bilinear')
+# Save figures:
figures = [manager.canvas.figure for manager in Gcf.get_all_fig_managers()]
-
for i, figure in enumerate(figures):
figure.savefig(PATH+'_figure{}.png'.format(i))
+plt.close('all')
+###################################################################################################
+# Close ups:
+dim_uv = (600, 200)
+mag_data_tip = MagData(mag_data.a, mag_data.magnitude[:, 600:, ...])
+pm = PMConvolve(mag_data.a, SimpleProjector(mag_data_tip.dim, 'y', dim_uv))
+phase_map_tip = PhaseMap(mag_data.a, pm(mag_data_tip).phase[350:, :])
+phase_map_tip.display_combined('Phase Map Nanowire Tip', density=density,
+ interpolation='bilinear')
+plt.savefig(PATH+'_nanowire_tip.png')
+mag_data_bot = MagData(mag_data.a, mag_data.magnitude[:, :300, ...])
+pm = PMConvolve(mag_data.a, SimpleProjector(mag_data_bot.dim, 'y', dim_uv))
+phase_map_tip = PhaseMap(mag_data.a, pm(mag_data_bot).phase[:250, :])
+phase_map_tip.display_combined('Phase Map Nanowire Bottom', density=density,
+ interpolation='bilinear')
+plt.savefig(PATH+'_nanowire_bot.png')
diff --git a/scripts/collaborations/zi_an.py b/scripts/collaborations/zi_an.py
index 6fd95b1..e03e2d3 100644
--- a/scripts/collaborations/zi_an.py
+++ b/scripts/collaborations/zi_an.py
@@ -9,18 +9,18 @@ import os
import numpy as np
-from PIL import Image
-
from pyDM3reader import DM3lib as dm3
+import matplotlib.pyplot as plt
+
import pyramid
-from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PMConvolve
-import pyramid.magcreator as mc
import pyramid.reconstruction as rc
+from time import clock
+
import logging
import logging.config
@@ -31,62 +31,67 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
PATH = '../../output/zi-an/'
-threshold = 3
+threshold = 1
a = 1.0 # in nm
density = 5
b_0 = 1
inter = 'none'
dim_small = (64, 64)
+smoothed_pictures = True
+lam = 7.5E-5
+order = 1
###################################################################################################
-im_2_ele = Image.open(PATH+'18a_0102ele_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.jpg').convert('L')
-im_2_mag = Image.open(PATH+'18a_0102mag_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.jpg').convert('L')
-im_4_ele = Image.open(PATH+'07_0102ele60_q3_pha_01_sb280_sc512_vf3_med5.jpg').convert('L')
-im_4_mag = Image.open(PATH+'07_0102mag60_q3_pha_01_sb280_sc512_vf3_med5.jpg').convert('L')
-
+# Read in files:
+if smoothed_pictures:
+ dm3_2_mag = dm3.DM3(PATH+'Output333_hw512.dm3').image
+ dm3_4_mag = dm3.DM3(PATH+'Output334_hw512.dm3').image
+else:
+ dm3_2_mag = dm3.DM3(PATH+'07_0102mag60_q3_pha_01_sb280_sc512_vf3_med5.dm3').image
+ dm3_4_mag = dm3.DM3(PATH+'18a_0102mag_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
dm3_2_ele = dm3.DM3(PATH+'07_0102ele60_q3_pha_01_sb280_sc512_vf3_med5.dm3').image
-dm3_2_mag = dm3.DM3(PATH+'07_0102mag60_q3_pha_01_sb280_sc512_vf3_med5.dm3').image
-dm3_2_neg = dm3.DM3(PATH+'07_0102neg60_q3_pha_01_sb280_sc512_vf3.dm3').image
-dm3_2_pos = dm3.DM3(PATH+'07_0102pos60_q3_pha_01_sb280_sc512_vf3.dm3').image
dm3_4_ele = dm3.DM3(PATH+'18a_0102ele_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
-dm3_4_mag = dm3.DM3(PATH+'18a_0102mag_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
-dm3_4_neg = dm3.DM3(PATH+'18a_0102neg_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
-dm3_4_pos = dm3.DM3(PATH+'18a_0102pos_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
-
-#phase_map_2 = PhaseMap(a, np.array(im_2_mag.resize(dim_small))/255.-0.09230)
-phase_map_2 = PhaseMap(a, np.array(dm3_2_mag.resize(dim_small))-0.09147)
+# Construct phase maps and masks
+phase_map_2 = PhaseMap(a, np.array(dm3_2_mag.resize(dim_small))+0.101)
phase_map_2.display_combined(density=density, interpolation=inter)
-mask_2 = np.expand_dims(np.where(np.array(dm3_2_ele.resize(dim_small)) > threshold, True, False),
- axis=0)
-#mask_2 = np.expand_dims(np.where(np.array(im_2_ele.resize(dim_small)) > threshold, True, False),
-# axis=0)
-
-#phase_map_4 = PhaseMap(a, np.array(im_4_mag.resize(dim_small))/255.-0.32197)
-phase_map_4 = PhaseMap(a, np.array(dm3_4_mag.resize(dim_small))-0.22569)
+plt.savefig(PATH+'phase_map_2part.png')
+mask_2 = np.expand_dims(np.where(np.array(dm3_2_ele.resize(dim_small)) > threshold,
+ True, False), axis=0)
+phase_map_4 = PhaseMap(a, np.array(dm3_4_mag.resize(dim_small))-2.546)
phase_map_4.display_combined(density=density, interpolation=inter)
-mask_4 = np.expand_dims(np.where(np.array(dm3_4_ele.resize(dim_small)) > threshold, True, False),
- axis=0)
-#mask_4 = np.expand_dims(np.where(np.array(im_4_ele.resize(dim_small)) > threshold, True, False),
-# axis=0)
-
+plt.savefig(PATH+'phase_map_4part.png')
+mask_4 = np.expand_dims(np.where(np.array(dm3_4_ele.resize(dim_small)) > threshold,
+ True, False), axis=0)
# Reconstruct the magnetic distribution:
-mag_data_rec_2 = rc.optimize_simple_leastsq(phase_map_2, mask_2, b_0)
-mag_data_rec_4 = rc.optimize_simple_leastsq(phase_map_4, mask_4, b_0)
-
+tic = clock()
+mag_data_rec_2 = rc.optimize_simple_leastsq(phase_map_2, mask_2, b_0, lam=lam, order=order)
+print '2 particle reconstruction time:', clock() - tic
+tic = clock()
+mag_data_rec_4 = rc.optimize_simple_leastsq(phase_map_4, mask_4, b_0, lam=lam, order=order)
+print '4 particle reconstruction time:', clock() - tic
# Display the reconstructed phase map and holography image:
phase_map_rec_2 = PMConvolve(a, SimpleProjector(mag_data_rec_2.dim), b_0)(mag_data_rec_2)
phase_map_rec_2.display_combined('Reconstr. Distribution', density=density, interpolation=inter)
+plt.savefig(PATH+'phase_map_2part_rec.png')
phase_map_rec_4 = PMConvolve(a, SimpleProjector(mag_data_rec_4.dim), b_0)(mag_data_rec_4)
phase_map_rec_4.display_combined('Reconstr. Distribution', density=density, interpolation=inter)
-
-(mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot()
-(mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot()
-
-(phase_map_rec_2-phase_map_2).display_phase('Difference')
-(phase_map_rec_4-phase_map_4).display_phase('Difference')
-
-print 'Average difference (2 cubes):', np.average((phase_map_rec_2-phase_map_2).phase)
-print 'Average difference (4 cubes):', np.average((phase_map_rec_4-phase_map_4).phase)
-
-mag_data_test_2 = MagData(a, mc.create_mag_dist_homog(mask_4, -0.7+np.pi))
-PMConvolve(a, SimpleProjector(mag_data_test_2.dim))(mag_data_test_2).display_phase()
+plt.savefig(PATH+'phase_map_4part_rec.png')
+# Plot the magnetization:
+axis = (mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot()
+axis.set_xlim(20, 45)
+axis.set_ylim(20, 45)
+plt.savefig(PATH+'mag_data_2part.png')
+axis = (mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot()
+axis.set_xlim(20, 45)
+axis.set_ylim(20, 45)
+plt.savefig(PATH+'mag_data_4part.png')
+# Display the Phase:
+phase_diff_2 = phase_map_rec_2-phase_map_2
+phase_diff_2.display_phase('Difference')
+plt.savefig(PATH+'phase_map_2part_diff.png')
+phase_diff_4 = phase_map_rec_4-phase_map_4
+phase_diff_4.display_phase('Difference')
+plt.savefig(PATH+'phase_map_4part_diff.png')
+# Get the average difference from the experimental results:
+print 'Average difference (2 cubes):', np.average(phase_diff_2.phase)
+print 'Average difference (4 cubes):', np.average(phase_diff_4.phase)
diff --git a/scripts/create distributions/create_sample.py b/scripts/create distributions/create_sample.py
index e7def12..b69ad17 100644
--- a/scripts/create distributions/create_sample.py
+++ b/scripts/create distributions/create_sample.py
@@ -22,10 +22,10 @@ logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
+###################################################################################################
# Input parameters:
-#--------------------------------------------------------------------------------------------------
key = 'sphere'
-#--------------------------------------------------------------------------------------------------
+###################################################################################################
filename = directory + '/mag_dist_' + key + '.txt'
dim = (64, 64, 64) # in px (z, y, x)
a = 1.0 # in nm
diff --git a/scripts/test methods/compare_pixel_fields.py b/scripts/test methods/compare_pixel_fields.py
index a99e2be..bda9f0c 100644
--- a/scripts/test methods/compare_pixel_fields.py
+++ b/scripts/test methods/compare_pixel_fields.py
@@ -31,8 +31,8 @@ if not os.path.exists(directory):
# Input parameters:
a = 1.0 # in nm
phi = 0 # in rad
-dim = (5, 5, 5)
-pixel = (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2))
+dim = (1, 8, 8)
+pixel = (0, int(dim[1]/2), int(dim[2]/2))
limit = 0.35
@@ -52,32 +52,29 @@ def get_fourier_kernel():
# Create magnetic data and projector:
-mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, pixel), phi, theta=0))
+mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, pixel), phi))
mag_data.save_to_llg(directory + '/mag_dist_single_pixel.txt')
-projector_ref = SimpleProjector(dim, 'z', None) # (SimpleProjector(dim)
-test_ref = np.array(projector_ref.weight.todense())
-projector = SimpleProjector(dim, 'x', (15, 15)) # (SimpleProjector(dim)
-test = np.array(projector.weight.todense())
+projector = SimpleProjector(dim)
# Kernel of a disc in real space:
phase_map_disc = PMConvolve(a, projector, geometry='disc')(mag_data)
phase_map_disc.unit = 'mrad'
phase_map_disc.display_phase('Phase of one Pixel (Disc)', limit=limit)
-## Kernel of a slab in real space:
-#phase_map_slab = PMConvolve(a, projector, geometry='slab')(mag_data)
-#phase_map_slab.unit = 'mrad'
-#phase_map_slab.display_phase('Phase of one Pixel (Slab)', limit=limit)
-## Kernel of the Fourier method:
-#phase_map_fft = PMFourier(a, projector, padding=0)(mag_data)
-#phase_map_fft.unit = 'mrad'
-#phase_map_fft.display_phase('Phase of one Pixel (Fourier)', limit=limit)
-## Kernel of the Fourier method, calculated directly:
-#phase_map_fft_kernel = PhaseMap(a, get_fourier_kernel(), 'mrad')
-#phase_map_fft_kernel.display_phase('Phase of one Pixel (Fourier Kernel)', limit=limit)
-## Kernel differences:
-#print 'Fourier Kernel, direct and indirect method are identical:', \
-# np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
-#(phase_map_disc-phase_map_fft).display_phase('Phase difference of one Pixel (Disc - Fourier)')
-#
+# Kernel of a slab in real space:
+phase_map_slab = PMConvolve(a, projector, geometry='slab')(mag_data)
+phase_map_slab.unit = 'mrad'
+phase_map_slab.display_phase('Phase of one Pixel (Slab)', limit=limit)
+# Kernel of the Fourier method:
+phase_map_fft = PMFourier(a, projector, padding=0)(mag_data)
+phase_map_fft.unit = 'mrad'
+phase_map_fft.display_phase('Phase of one Pixel (Fourier)', limit=limit)
+# Kernel of the Fourier method, calculated directly:
+phase_map_fft_kernel = PhaseMap(a, get_fourier_kernel(), 'mrad')
+phase_map_fft_kernel.display_phase('Phase of one Pixel (Fourier Kernel)', limit=limit)
+# Kernel differences:
+print 'Fourier Kernel, direct and indirect method are identical:', \
+ np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
+(phase_map_disc-phase_map_fft).display_phase('Phase difference of one Pixel (Disc - Fourier)')
+
## Cross section plots of real space kernels:
#fig = plt.figure()
#axis = fig.add_subplot(1, 1, 1)
diff --git a/scripts/test methods/projector_test.py b/scripts/test methods/projector_test.py
new file mode 100644
index 0000000..616e69e
--- /dev/null
+++ b/scripts/test methods/projector_test.py
@@ -0,0 +1,29 @@
+# -*- coding: utf-8 -*-
+"""Created on Fri Aug 15 08:09:10 2014 @author: Jan"""
+
+
+import numpy as np
+
+import pyramid.magcreator as mc
+from pyramid.magdata import MagData
+from pyramid.projector import SimpleProjector, XTiltProjector, YTiltProjector
+from pyramid.phasemapper import PMConvolve
+
+
+dim = (16, 24, 32)
+center = (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2))
+width = (2, 8, 16)
+dim_uv = (32, 48)
+a = 2
+b_0 = 1.5
+
+shape = mc.Shapes.slab(dim, center, width)
+mag_data = MagData(a, mc.create_mag_dist_homog(shape, phi=np.pi/6, theta=np.pi/3))
+
+# PROJECTOR TESTING
+PMConvolve(a, SimpleProjector(dim, 'z', dim_uv), b_0)(mag_data).display_phase('z simple')
+PMConvolve(a, XTiltProjector(dim, 0, dim_uv), b_0)(mag_data).display_phase('z (x-tilt)')
+PMConvolve(a, SimpleProjector(dim, 'x', dim_uv), b_0)(mag_data).display_phase('x simple')
+PMConvolve(a, YTiltProjector(dim, np.pi/2, dim_uv), b_0)(mag_data).display_phase('x (y-tilt)')
+PMConvolve(a, XTiltProjector(dim, np.pi/2, dim_uv), b_0)(mag_data).display_phase('y (x-tilt)')
+PMConvolve(a, SimpleProjector(dim, 'y', dim_uv), b_0)(mag_data).display_phase('y simple')
--
GitLab