From 622c59d868685f7364f46548329c5cfac936a09d Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Thu, 16 Jul 2015 15:20:35 +0200
Subject: [PATCH] Phase offset and ramp can now be fitted! Many updates and
fixes to plotting! colormap: New module with colormap subclass which has
special functions to calculate directional rgb-colors, now used in
many plots! __init__: On star import now also has abbreviations for a few
modules. fwd_model: Now has flags for fitting an offset and/or ramp.
CAUTION: will be moved next update to new class! magdata: Removed
colorwheel stuff (all moved to 'colormap'). Overhaul of quiver plot
in 2D and 3D. Now shows mask and uses new colormap for angular color
coding. phasemap: Now uses the new colormap for holo plots. Removed color
wheel (moved to the new class). Added function to add ramps.
phasemapper: Now handles also offsets and ramps in the phase.
CAUTION: will be moved next update to new class! projector: Fixed bugs in
SimpleProjector. 1) leftmost pixels are no longer moved moved to
left corner of padded projection (when dim_uv is not None). 2) odd
padding when using dim_uv are now handled properly (one more pixel
on top and right if necessary). costfunction: n is now taken from the
fwd_model instead of the data_set. dataset: Introduced the count variable
which gives the number of images. magcreator: In the process of updating,
more on next commit. reconstruction: Also handles ramps and offsets now,
properly extracts them after reconstruction. regularisator:
New argument 'add_params' is used to cut off input (for offset
and ramp) which are not used in the regularisation. scripts: several small
changes...
---
pyramid/__init__.py | 7 +-
pyramid/colormap.py | 163 ++++++++++++++++++
pyramid/costfunction.py | 4 +-
pyramid/dataset.py | 4 +
pyramid/forwardmodel.py | 66 +++++--
pyramid/magcreator.py | 20 ++-
pyramid/magdata.py | 162 +++++++++--------
pyramid/phasemap.py | 162 ++++++++++-------
pyramid/phasemapper.py | 58 +++++--
pyramid/projector.py | 16 +-
pyramid/reconstruction.py | 15 +-
pyramid/regularisator.py | 37 ++--
.../magdata/magcreator/create_homog_slab.py | 4 +-
scripts/magdata/magdata_from_vtk.py | 73 ++++++--
scripts/phasemap/phasemap_from_dm3.py | 11 +-
scripts/phasemap/phasemap_from_image.py | 10 +-
.../reconstruction_2d_from_phasemap.py | 62 ++++---
.../reconstruction_3d_from_magdata.py | 60 +++++--
scripts/temp/custom_colormap.py | 2 +-
scripts/temp/test_dim_uv_projection.py | 24 +++
20 files changed, 715 insertions(+), 245 deletions(-)
create mode 100644 pyramid/colormap.py
create mode 100644 scripts/temp/test_dim_uv_projection.py
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index 1850baf..9b7f840 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -46,8 +46,11 @@ numcore
from . import analytic
+from . import analytic as an
from . import magcreator
+from . import magcreator as mc
from . import reconstruction
+from . import reconstruction as rc
from . import fft
from .costfunction import * # analysis:ignore
from .dataset import * # analysis:ignore
@@ -60,6 +63,7 @@ from .phasemapper import * # analysis:ignore
from .projector import * # analysis:ignore
from .regularisator import * # analysis:ignore
from .quaternion import * # analysis:ignore
+from .colormap import * # analysis:ignore
from .config import * # analysis:ignore
from .version import version as __version__
from .version import hg_revision as __hg_revision__
@@ -69,7 +73,7 @@ _log = logging.getLogger(__name__)
_log.info("Starting PYRAMID V{} HG{}".format(__version__, __hg_revision__))
del logging
-__all__ = ['analytic', 'magcreator', 'reconstruction', 'fft']
+__all__ = ['analytic', 'magcreator', 'reconstruction', 'fft', 'an', 'mc', 'rc']
__all__.extend(costfunction.__all__)
__all__.extend(dataset.__all__)
__all__.extend(diagnostics.__all__)
@@ -81,3 +85,4 @@ __all__.extend(phasemapper.__all__)
__all__.extend(projector.__all__)
__all__.extend(regularisator.__all__)
__all__.extend(quaternion.__all__)
+__all__.extend(colormap.__all__)
diff --git a/pyramid/colormap.py b/pyramid/colormap.py
new file mode 100644
index 0000000..94f3e49
--- /dev/null
+++ b/pyramid/colormap.py
@@ -0,0 +1,163 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jul 08 15:43:06 2015
+
+@author: Jan
+"""
+
+# TODO: Docstrings!
+
+
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from PIL import Image
+
+import logging
+
+
+__all__ = ['DirectionalColormap']
+
+
+class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
+
+ _log = logging.getLogger(__name__+'.DirectionalColormap')
+
+ CDICT = {'red': [(0.00, 0.0, 0.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 0.0, 0.0)],
+
+ 'green': [(0.00, 1.0, 1.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 1.0, 1.0)],
+
+ 'blue': [(0.00, 0.0, 0.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 0.0, 0.0)]}
+
+ CDICT_INV = {'red': [(0.00, 1.0, 1.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 1.0, 1.0)],
+
+ 'green': [(0.00, 0.0, 0.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 0.0, 0.0)],
+
+ 'blue': [(0.00, 1.0, 1.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 1.0, 1.0)]}
+
+ HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
+ HOLO_CMAP_INV = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT_INV, 256)
+
+
+ def __init__(self, inverted=False):
+ self._log.debug('Calling __create_directional_colormap')
+# r, g, b = [], [], [] # RGB lists
+# # Create saturation lists to encode up and down directions via black and white colors.
+# # example for 5 levels from black (down) to color (in-plane) to white:
+# # pos_sat: [ 0. 0.5 1. 1. 1. ]
+# # neg_sat: [ 0. 0. 0. 0.5 1. ]
+# center = levels//2
+# pos_sat = np.ones(levels)
+# pos_sat[0:center] = [i/center for i in range(center)]
+# neg_sat = np.zeros(levels)
+# neg_sat[center+1:] = [(i+1)/center for i in range(center)]
+#
+# # Iterate over all levels (the center level represents in-plane moments!):
+# for i in range(levels):
+# inter_points = np.linspace(i/levels, (i+1)/levels, 5) # interval points, current level
+# # Red:
+# r.append((inter_points[0], 0, neg_sat[i]))
+# r.append((inter_points[1], pos_sat[i], pos_sat[i]))
+# r.append((inter_points[2], pos_sat[i], pos_sat[i]))
+# r.append((inter_points[3], neg_sat[i], neg_sat[i]))
+# r.append((inter_points[4], neg_sat[i], 0))
+# # Green:
+# g.append((inter_points[0], 0, neg_sat[i]))
+# g.append((inter_points[1], neg_sat[i], neg_sat[i]))
+# g.append((inter_points[2], pos_sat[i], pos_sat[i]))
+# g.append((inter_points[3], pos_sat[i], pos_sat[i]))
+# g.append((inter_points[4], neg_sat[i], 0))
+# # Blue
+# b.append((inter_points[0], 0, pos_sat[i]))
+# b.append((inter_points[1], neg_sat[i], neg_sat[i]))
+# b.append((inter_points[2], neg_sat[i], neg_sat[i]))
+# b.append((inter_points[3], neg_sat[i], neg_sat[i]))
+# b.append((inter_points[4], pos_sat[i], 0))
+# # Combine to color dictionary and return:
+# cdict = {'red': r, 'green': g, 'blue': b}
+ cdict = self.CDICT_INV if inverted else self.CDICT
+ super(DirectionalColormap, self).__init__('directional_colormap', cdict, N=256)
+ self._log.debug('Created '+str(self))
+
+ @classmethod
+ def display_colorwheel(cls, mode='white_to_color'):
+ '''Display a color wheel to illustrate the color coding of the gradient direction.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ None
+
+ ''' # TODO: mode docstring!
+ cls._log.debug('Calling display_color_wheel')
+ yy, xx = np.indices((512, 512)) - 256
+ r = np.hypot(xx, yy)
+ # Create the wheel:
+ colorind = (1 + np.arctan2(yy, xx)/np.pi) / 2
+ saturation = r / 256 # 0 in center, 1 at borders of circle!
+ if mode == 'black_to_color':
+ pass
+ elif mode == 'color_to_black':
+ saturation = 1 - saturation
+ elif mode == 'white_to_color':
+ saturation = 2 - saturation
+ elif mode == 'white_to_color_to_black':
+ saturation = 2 - 2*saturation
+ else:
+ raise Exception('Invalid color wheel mode!')
+ saturation *= (r <= 256) # TODO: [r<=256]
+ rgb = cls.rgb_from_colorind_and_saturation(colorind, saturation)
+ color_wheel = Image.fromarray(rgb)
+ # Plot the color wheel:
+ fig = plt.figure(figsize=(4, 4))
+ axis = fig.add_subplot(1, 1, 1, aspect='equal')
+ axis.imshow(color_wheel, origin='lower')
+ plt.tick_params(axis='both', which='both', labelleft='off', labelbottom='off',
+ left='off', right='off', top='off', bottom='off')
+
+ @classmethod # TODO: Docstrings!
+ def rgb_from_colorind_and_saturation(cls, colorind, saturation):
+ c, s = colorind, saturation
+ rgb_norm = (np.minimum(s, np.ones_like(s)).T * cls.HOLO_CMAP(c)[..., :3].T).T
+ rgb_invs = (np.maximum(s-1, np.zeros_like(s)).T * cls.HOLO_CMAP_INV(c)[..., :3].T).T
+ return (255.999 * (rgb_norm + rgb_invs)).astype(np.uint8)
+
+ @classmethod
+ def rgb_from_angles(cls, phi, theta=np.pi/2):
+ colorind = (1 + phi/np.pi) / 2
+ saturation = 2 * (1 - theta/np.pi) # goes from 0 to 2 # TODO: explain!
+ return cls.rgb_from_colorind_and_saturation(colorind, saturation)
+
+ @classmethod
+ def rgb_from_direction(cls, x, y, z):
+ phi = np.arctan2(y, x)
+ r = np.sqrt(x**2 + y**2 + z**2)
+ theta = np.arccos(z / (r+1E-30))
+ return cls.rgb_from_angles(phi, theta)
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 41fb1b5..03ddcde 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -57,8 +57,8 @@ class Costfunction(object):
# Extract important information:
data_set = fwd_model.data_set
self.y = data_set.phase_vec
- self.n = data_set.n
- self.m = data_set.m
+ self.n = fwd_model.n
+ self.m = fwd_model.m
if data_set.Se_inv is None:
data_set.set_Se_inv_diag_with_conf()
self.Se_inv = data_set.Se_inv
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index bf73a19..a31802c 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -94,6 +94,10 @@ class DataSet(object):
def n(self):
return 3 * np.sum(self.mask)
+ @property
+ def count(self):
+ return len(self.projectors)
+
@property
def phase_vec(self):
return np.concatenate([p.phase_vec for p in self.phase_maps])
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 41a9e6e..24be8ad 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -46,12 +46,19 @@ class ForwardModel(object):
_log = logging.getLogger(__name__+'.ForwardModel')
- def __init__(self, data_set):
+ def __init__(self, data_set, fit_ramps=False, fit_offsets=False):
self._log.debug('Calling __init__')
self.data_set = data_set
+ if fit_ramps: # The ramps are not fitted without the offsets!
+ # TODO: fit über String -> eine Flag!
+ # TODO: immer mit self!
+ fit_offsets = True
+ self.fit_ramps = fit_ramps
+ self.fit_offsets = fit_offsets
self.phase_mappers = data_set.phase_mappers
self.m = data_set.m
- self.n = data_set.n
+ self.n = data_set.n + fit_offsets * data_set.count + fit_ramps * 2 * data_set.count
+ # TODO: bools nicht als integer verwenden!
self.shape = (self.m, self.n)
self.hook_points = data_set.hook_points
self.mag_data = MagData(data_set.a, np.zeros((3,)+data_set.dim))
@@ -82,13 +89,24 @@ class ForwardModel(object):
self._log.debug('Calling __str__')
return 'ForwardModel(data_set=%s)' % (self.data_set)
+# TODO: offset und ramp HIER handeln!
+
def __call__(self, x):
+ count = self.data_set.count
+ offsets = [None] * count
+ ramps = [None] * count
+ if self.fit_ramps:
+ x, offsets, u_ramps, v_ramps = np.split(x, [-3*count, -2*count, -count])
+ ramps = zip(u_ramps, v_ramps)
+ elif self.fit_offsets:
+ x, offsets = np.split(x, [-count])
self.mag_data.magnitude[...] = 0
self.mag_data.set_vector(x, self.data_set.mask)
result = np.zeros(self.m)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
- phase_map = self.phase_mappers[projector.dim_uv](projector(self.mag_data))
+ mapper = self.phase_mappers[projector.dim_uv]
+ phase_map = mapper(projector(self.mag_data), offsets[i], ramps[i])
result[hp[i]:hp[i+1]] = phase_map.phase_vec
return np.reshape(result, -1)
###################################################################################################
@@ -126,19 +144,27 @@ class ForwardModel(object):
`vector`.
'''
+ count = self.data_set.count
+ offsets = [None] * count
+ ramps = [None] * count
+ if self.fit_ramps:
+ vector, offsets, u_ramps, v_ramps = np.split(vector, [-3*count, -2*count, -count])
+ ramps = zip(u_ramps, v_ramps)
+ elif self.fit_offsets:
+ vector, offsets = np.split(vector, [-count])
self.mag_data.magnitude[...] = 0
self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.m)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
mag_vec = self.mag_data.mag_vec
- res = self.phase_mappers[projector.dim_uv].jac_dot(projector.jac_dot(mag_vec))
+ mapper = self.phase_mappers[projector.dim_uv]
+ res = mapper.jac_dot(projector.jac_dot(mag_vec), offsets[i], ramps[i])
result[hp[i]:hp[i+1]] = res
return result
def _jac_dot_element(self, mag_vec, projector, phasemapper):
- return phasemapper.jac_dot(projector.jac_dot(mag_vec))
-
+ return phasemapper.jac_dot(projector.jac_dot(mag_vec)) # TODO: ???
def jac_T_dot(self, x, vector):
''''Calculate the product of the transposed Jacobi matrix with a given `vector`.
@@ -159,10 +185,30 @@ class ForwardModel(object):
the input `vector`.
'''
- result = np.zeros(3*np.prod(self.data_set.dim))
+ offsets = []
+ u_ramps = []
+ v_ramps = []
+ proj_T_result = np.zeros(3*np.prod(self.data_set.dim))
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
vec = vector[hp[i]:hp[i+1]]
- result += projector.jac_T_dot(self.phase_mappers[projector.dim_uv].jac_T_dot(vec))
- self.mag_data.mag_vec = result
- return self.mag_data.get_vector(self.data_set.mask)
+ mapper = self.phase_mappers[projector.dim_uv]
+ map_T_result = mapper.jac_T_dot(vec, self.fit_ramps, self.fit_offsets)
+ if self.fit_ramps: # Extract offset and ramps (transposed):
+ map_T_result, add_params = np.split(map_T_result, [-3])
+ offsets.append(add_params[0])
+ u_ramps.append(add_params[1])
+ v_ramps.append(add_params[2])
+ elif self.fit_offsets: # Extract offset (transposed):
+ map_T_result, offset = np.split(map_T_result, [-1])
+ offsets.append(offset)
+ proj_T_result += projector.jac_T_dot(map_T_result)
+ self.mag_data.mag_vec = proj_T_result
+ result = self.mag_data.get_vector(self.data_set.mask)
+ if self.fit_ramps:
+ return np.concatenate((result, np.reshape(offsets, -1),
+ np.reshape(u_ramps, -1), np.reshape(v_ramps, -1)))
+ elif self.fit_offsets:
+ return np.concatenate((result, np.reshape(offsets, -1)))
+ else:
+ return result
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index 3083d45..a777709 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -16,6 +16,9 @@ center.
"""
+
+from __future__ import division
+
import numpy as np
from numpy import pi
@@ -65,12 +68,17 @@ class Shapes(object):
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a tuple of length 3!'
- mag_shape = np.array([[[abs(x - center[2]) <= width[2] / 2
- and abs(y - center[1]) <= width[1] / 2
- and abs(z - center[0]) <= width[0] / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ zz, yy, xx = np.indices(dim) + 0.5
+ xx_shape = np.where(abs(xx-center[2]) <= width[2]/2, True, False)
+ yy_shape = np.where(abs(yy-center[1]) <= width[1]/2, True, False)
+ zz_shape = np.where(abs(zz-center[0]) <= width[0]/2, True, False)
+ mag_shape = np.logical_and.reduce((xx_shape, yy_shape, zz_shape))
+# mag_shape = np.array([[[abs(x - center[2]) <= width[2] / 2
+# and abs(y - center[1]) <= width[1] / 2
+# and abs(z - center[0]) <= width[0] / 2
+# for x in range(dim[2])]
+# for y in range(dim[1])]
+# for z in range(dim[0])])
return mag_shape
@classmethod
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index af24642..04edc61 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -19,6 +19,7 @@ import matplotlib.cm as cmx
from matplotlib.ticker import MaxNLocator
from pyramid import fft
+from pyramid.colormap import DirectionalColormap
from numbers import Number
@@ -167,44 +168,44 @@ class MagData(object):
self._log.debug('Calling __imul__')
return self.__mul__(other)
- @classmethod
- def _create_directional_colormap(cls, levels=15, N=256):
- cls._log.debug('Calling __create_directional_colormap')
- r, g, b = [], [], [] # RGB lists
- # Create saturation lists to encode up and down directions via black and white colors.
- # example for 5 levels from black (down) to color (in-plane) to white:
- # pos_sat: [ 0. 0.5 1. 1. 1. ]
- # neg_sat: [ 0. 0. 0. 0.5 1. ]
- center = levels//2
- pos_sat = np.ones(levels)
- pos_sat[0:center] = [i/center for i in range(center)]
- neg_sat = np.zeros(levels)
- neg_sat[center+1:] = [(i+1)/center for i in range(center)]
-
- # Iterate over all levels (the center level represents in-plane moments!):
- for i in range(levels):
- inter_points = np.linspace(i/levels, (i+1)/levels, 5) # interval points, current level
- # Red:
- r.append((inter_points[0], 0, neg_sat[i]))
- r.append((inter_points[1], pos_sat[i], pos_sat[i]))
- r.append((inter_points[2], pos_sat[i], pos_sat[i]))
- r.append((inter_points[3], neg_sat[i], neg_sat[i]))
- r.append((inter_points[4], neg_sat[i], 0))
- # Green:
- g.append((inter_points[0], 0, neg_sat[i]))
- g.append((inter_points[1], neg_sat[i], neg_sat[i]))
- g.append((inter_points[2], pos_sat[i], pos_sat[i]))
- g.append((inter_points[3], pos_sat[i], pos_sat[i]))
- g.append((inter_points[4], neg_sat[i], 0))
- # Blue
- b.append((inter_points[0], 0, pos_sat[i]))
- b.append((inter_points[1], neg_sat[i], neg_sat[i]))
- b.append((inter_points[2], neg_sat[i], neg_sat[i]))
- b.append((inter_points[3], neg_sat[i], neg_sat[i]))
- b.append((inter_points[4], pos_sat[i], 0))
- # Combine to color dictionary and return:
- cdict = {'red': r, 'green': g, 'blue': b}
- return mpl.colors.LinearSegmentedColormap('directional_colormap', cdict, N)
+# @classmethod
+# def _create_directional_colormap(cls, levels=15, N=256):
+# cls._log.debug('Calling __create_directional_colormap')
+# r, g, b = [], [], [] # RGB lists
+# # Create saturation lists to encode up and down directions via black and white colors.
+# # example for 5 levels from black (down) to color (in-plane) to white:
+# # pos_sat: [ 0. 0.5 1. 1. 1. ]
+# # neg_sat: [ 0. 0. 0. 0.5 1. ]
+# center = levels//2
+# pos_sat = np.ones(levels)
+# pos_sat[0:center] = [i/center for i in range(center)]
+# neg_sat = np.zeros(levels)
+# neg_sat[center+1:] = [(i+1)/center for i in range(center)]
+#
+# # Iterate over all levels (the center level represents in-plane moments!):
+# for i in range(levels):
+# inter_points = np.linspace(i/levels, (i+1)/levels, 5) # interval points, current level
+# # Red:
+# r.append((inter_points[0], 0, neg_sat[i]))
+# r.append((inter_points[1], pos_sat[i], pos_sat[i]))
+# r.append((inter_points[2], pos_sat[i], pos_sat[i]))
+# r.append((inter_points[3], neg_sat[i], neg_sat[i]))
+# r.append((inter_points[4], neg_sat[i], 0))
+# # Green:
+# g.append((inter_points[0], 0, neg_sat[i]))
+# g.append((inter_points[1], neg_sat[i], neg_sat[i]))
+# g.append((inter_points[2], pos_sat[i], pos_sat[i]))
+# g.append((inter_points[3], pos_sat[i], pos_sat[i]))
+# g.append((inter_points[4], neg_sat[i], 0))
+# # Blue
+# b.append((inter_points[0], 0, pos_sat[i]))
+# b.append((inter_points[1], neg_sat[i], neg_sat[i]))
+# b.append((inter_points[2], neg_sat[i], neg_sat[i]))
+# b.append((inter_points[3], neg_sat[i], neg_sat[i]))
+# b.append((inter_points[4], pos_sat[i], 0))
+# # Combine to color dictionary and return:
+# cdict = {'red': r, 'green': g, 'blue': b}
+# return mpl.colors.LinearSegmentedColormap('directional_colormap', cdict, N)
def copy(self):
'''Returns a copy of the :class:`~.MagData` object
@@ -335,7 +336,7 @@ class MagData(object):
assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
cv[2*i:2*(i+1)] = values
cv *= np.resize([1, -1], len(cv))
- cv = np.where(crop_values == 0, None, cv)
+ cv = np.where(cv == 0, None, cv)
self.magnitude = self.magnitude[:, cv[0]:cv[1], cv[2]:cv[3], cv[4]:cv[5]]
def get_mask(self, threshold=0):
@@ -667,7 +668,8 @@ class MagData(object):
tree.write(filename, pretty_print=True)
def quiver_plot(self, title='Magnetization Distribution', axis=None, proj_axis='z',
- coloring='angle', ar_dens=1, ax_slice=None, log=False, scaled=True, scale=1.):
+ coloring='angle', ar_dens=1, ax_slice=None, log=False, scaled=True,
+ scale=1., show_mask=True):
'''Plot a slice of the magnetization as a quiver plot.
Parameters
@@ -700,6 +702,7 @@ class MagData(object):
The axis on which the graph is plotted.
'''
+ # TODO: include mask if available!
self._log.debug('Calling quiver_plot')
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
@@ -707,50 +710,53 @@ class MagData(object):
self._log.debug('proj_axis == z')
if ax_slice is None:
self._log.debug('ax_slice is None')
- ax_slice = int(self.dim[0]/2.)
+ ax_slice = self.dim[0] // 2
plot_u = np.copy(self.magnitude[0][ax_slice, ...]) # x-component
plot_v = np.copy(self.magnitude[1][ax_slice, ...]) # y-component
u_label = 'x [px]'
v_label = 'y [px]'
+ submask = self.get_mask()[ax_slice, ...]
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
self._log.debug('proj_axis == y')
if ax_slice is None:
self._log.debug('ax_slice is None')
- ax_slice = int(self.dim[1]/2.)
+ ax_slice = self.dim[1] // 2
plot_u = np.copy(self.magnitude[0][:, ax_slice, :]) # x-component
plot_v = np.copy(self.magnitude[2][:, ax_slice, :]) # z-component
u_label = 'x [px]'
v_label = 'z [px]'
+ submask = self.get_mask()[:, ax_slice, :]
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
self._log.debug('proj_axis == x')
if ax_slice is None:
self._log.debug('ax_slice is None')
- ax_slice = int(self.dim[2]/2.)
+ ax_slice = self.dim[2] // 2
plot_u = np.swapaxes(np.copy(self.magnitude[2][..., ax_slice]), 0, 1) # z-component
plot_v = np.swapaxes(np.copy(self.magnitude[1][..., ax_slice]), 0, 1) # y-component
u_label = 'z [px]'
v_label = 'y [px]'
+ submask = self.get_mask()[..., ax_slice]
# Prepare quiver (select only used arrows if ar_dens is specified):
dim_uv = plot_u.shape
- yy, xx = np.indices(dim_uv)
- xx = xx[::ar_dens, ::ar_dens]
- yy = yy[::ar_dens, ::ar_dens]
+ vv, uu = np.indices(dim_uv) + 0.5 # shift to center of pixel
+ uu = uu[::ar_dens, ::ar_dens]
+ vv = vv[::ar_dens, ::ar_dens]
plot_u = plot_u[::ar_dens, ::ar_dens]
plot_v = plot_v[::ar_dens, ::ar_dens]
amplitudes = np.hypot(plot_u, plot_v)
- angles = np.angle(plot_u+1j*plot_v, deg=True)
+ angles = np.angle(plot_u+1j*plot_v, deg=True).tolist()
# Calculate the arrow colors:
if coloring == 'angle':
self._log.debug('Encoding angles')
- colors = (1 - np.arctan2(plot_v, plot_u)/np.pi) / 2 # in-plane angle (rescaled: 0 - 1)
- cmap = self._create_directional_colormap(levels=1, N=256)
+ colorinds = (1 + np.arctan2(plot_v, plot_u)/np.pi) / 2 # in-plane color index (0 - 1)
+ cmap = DirectionalColormap()
elif coloring == 'amplitude':
self._log.debug('Encoding amplitude')
- colors = amplitudes / amplitudes.max()
+ colorinds = amplitudes / amplitudes.max()
cmap = 'jet'
elif coloring == 'uniform':
self._log.debug('No color encoding')
- colors = np.zeros_like(plot_u) # use black arrows!
+ colorinds = np.zeros_like(plot_u) # use black arrows!
cmap = 'gray'
else:
raise AttributeError("Invalid coloring mode! Use 'angles', 'amplitude' or 'uniform'!")
@@ -772,19 +778,27 @@ class MagData(object):
amplitudes = np.hypot(plot_u, plot_v) # Recalculate!
# Scale the magnitude of the arrows to the highest one (if specified):
if scaled:
- plot_u /= amplitudes.max()
- plot_v /= amplitudes.max()
- axis.quiver(xx, yy, plot_u, plot_v, colors, cmap=cmap, angles=angles,
+ plot_u /= amplitudes.max() + 1E-30
+ plot_v /= amplitudes.max() + 1E-30
+ axis.quiver(uu, vv, plot_u, plot_v, colorinds, cmap=cmap, clim=(0, 1), angles=angles,
pivot='middle', units='xy', scale_units='xy', scale=scale/ar_dens,
minlength=0.25, headwidth=6, headlength=7)
- axis.set_xlim(-1, dim_uv[1])
- axis.set_ylim(-1, dim_uv[0])
+ if show_mask and not np.all(submask): # Plot mask if desired and not trivial!
+ vv, uu = np.indices(dim_uv) + 0.5 # shift to center of pixel
+ axis.contour(uu, vv, submask, levels=[0.5], colors='k', linestyles='dotted')
+ axis.set_xlim(0, dim_uv[1])
+ axis.set_ylim(0, dim_uv[0])
axis.set_title(title, fontsize=18)
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
- axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
- axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ if dim_uv[0] >= dim_uv[1]:
+ u_bin, v_bin = np.max((2, np.floor(9*dim_uv[1]/dim_uv[0]))), 9
+ else:
+ u_bin, v_bin = 9, np.max((2, np.floor(9*dim_uv[0]/dim_uv[1])))
+ axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
+ axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
+ # Return plotting axis:
return axis
def quiver_plot3d(self, title='Magnetization Distribution', limit=None, cmap='jet',
@@ -819,6 +833,7 @@ class MagData(object):
'''
self._log.debug('Calling quiver_plot3D')
from mayavi import mlab
+ import warnings
a = self.a
dim = self.dim
if limit is None:
@@ -834,21 +849,28 @@ class MagData(object):
z_mag = self.magnitude[2][::ad, ::ad, ::ad].flatten()
# Plot them as vectors:
mlab.figure(size=(750, 700))
- plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode=mode, colormap=cmap)
+ with warnings.catch_warnings(record=True) as w: # Catch annoying warning
+ plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode=mode, colormap=cmap)
+ if len(w) != 0: # If a warning occurs make sure it is (only) the expected one!
+ assert len(w) == 1
+ assert issubclass(w[0].category, FutureWarning)
+ assert 'comparison' in str(w[0].message)
if coloring == 'angle': # Encodes the full angle via colorwheel and saturation
self._log.debug('Encoding full 3D angles')
from tvtk.api import tvtk
- phis = (1 - np.arctan2(y_mag, x_mag)/np.pi) / 2 # in-plane angle (rescaled: 0 to 1)
- thetas = np.arctan2(np.hypot(y_mag, x_mag), z_mag)/np.pi # vert. angle (also rescaled)
- levels = 15 # number of shades for up/down arrows (white is up, black is down)
- N = 256 # Overall number of colors for the colormap (only N/levels per level!)
- cmap = self._create_directional_colormap(levels, N)
- colors = []
- for i in range(len(x_mag)):
- level = np.floor((1-thetas[i]) * levels)
- lookup_value = (level + phis[i]) / levels
- rgba = tuple((np.asarray(cmap(lookup_value)) * (N-1)).astype(np.int))
- colors.append(rgba)
+ colorinds = DirectionalColormap.rgb_from_direction(x_mag, y_mag, z_mag)
+ colors = map(tuple, colorinds) # convert to list of tuples!
+# phis = (1 - np.arctan2(y_mag, x_mag)/np.pi) / 2 # in-plane angle (rescaled: 0 to 1)
+# thetas = np.arctan2(np.hypot(y_mag, x_mag), z_mag)/np.pi # vert. angle (also rescaled)
+# levels = 31 # number of shades for up/down arrows (white is up, black is down)
+# N = 2048 # Overall number of colors for the colormap (only N/levels per level!)
+# cmap = DirectionalColormap(levels, N)#self._create_directional_colormap(levels, N)
+# colors = []
+# for i in range(len(x_mag)):
+# level = np.floor((1-thetas[i]) * levels)
+# lookup_value = (level + phis[i]) / levels
+# rgba = tuple((np.asarray(cmap(lookup_value)) * 255.99999).astype(np.int))
+# colors.append(rgba)
sc = tvtk.UnsignedCharArray() # Used to hold the colors
sc.from_array(colors)
plot.mlab_source.dataset.point_data.scalars = sc
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index 01bdec7..033cf43 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -8,13 +8,14 @@
import os
import numpy as np
-from numpy import pi
from scipy.ndimage.interpolation import zoom
+from pyramid.colormap import DirectionalColormap
+
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
-from matplotlib.ticker import NullLocator, MaxNLocator, FuncFormatter
+from matplotlib.ticker import MaxNLocator, FuncFormatter
from PIL import Image
from numbers import Number
@@ -218,6 +219,7 @@ class PhaseMap(object):
else: # other is a Number
self._log.debug('Adding an offset')
return PhaseMap(self.a, self.phase+other, self.mask, self.confidence, self.unit)
+ # TODO: Add fitting numpy array! (useful for ramps)!
def __sub__(self, other): # self - other
self._log.debug('Calling __sub__')
@@ -254,6 +256,12 @@ class PhaseMap(object):
self._log.debug('Calling __imul__')
return self.__mul__(other)
+ def add_ramp(self, offset=0, ramp=(0, 0)):
+ # TODO: Docstring!
+ self._log.debug('Calling add_ramp')
+ vv, uu = self.a * np.indices(self.dim_uv)
+ self.phase += offset + ramp[0]*uu + ramp[1]*vv
+
def copy(self):
'''Returns a copy of the :class:`~.PhaseMap` object
@@ -399,7 +407,7 @@ class PhaseMap(object):
assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
cv[2*i:2*(i+1)] = values
cv *= np.resize([1, -1], len(cv))
- cv = np.where(crop_values == 0, None, cv)
+ cv = np.where(cv == 0, None, cv)
self.phase = self.phase[cv[0]:cv[1], cv[2]:cv[3]]
self.mask = self.mask[cv[0]:cv[1], cv[2]:cv[3]]
self.confidence = self.confidence[cv[0]:cv[1], cv[2]:cv[3]]
@@ -583,12 +591,17 @@ class PhaseMap(object):
# Plot the phasemap:
im = axis.pcolormesh(phase, cmap=cmap, vmin=-limit, vmax=limit, norm=norm)
if show_mask and not np.all(self.mask): # Plot mask if desired and not trivial!
- axis.contour(self.mask, levels=[0.999999], colors='g')
+ vv, uu = np.indices(self.dim_uv) + 0.5
+ axis.contour(uu, vv, self.mask, levels=[0.5], colors='k', linestyles='dotted')
# Set the axes ticks and labels:
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))
+ if self.dim_uv[0] >= self.dim_uv[1]:
+ u_bin, v_bin = np.max((2, np.floor(9*self.dim_uv[1]/self.dim_uv[0]))), 9
+ else:
+ u_bin, v_bin = 9, np.max((2, np.floor(9*self.dim_uv[0]/self.dim_uv[1])))
+ axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
+ axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
axis.tick_params(axis='both', which='major', labelsize=14)
@@ -673,38 +686,62 @@ class PhaseMap(object):
self._log.debug('Calling display_holo')
# Calculate gain if 'auto' is selected:
if gain == 'auto':
- gain = 4 * 2*pi/(np.abs(self.phase).max()+1E-30)
+ gain = 4 * 2*np.pi/(np.abs(self.phase).max()+1E-30)
# Set title if not set:
if title is None:
title = 'Contour Map (gain: %.2g)' % gain
# Calculate the holography image intensity:
- img_holo = (1 + np.cos(gain * self.phase)) / 2
+# img_holo = (1 + np.cos(gain * self.phase)) / 2
+ holo = (1 + np.cos(gain * self.phase)) / 2
# Calculate the phase gradients, expressed by magnitude and angle:
- phase_grad_y, phase_grad_x = np.gradient(self.phase, self.a, self.a)
- phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
- phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)
+ phase_grad_x, phase_grad_y = np.gradient(self.phase, self.a, self.a)
+# plt.figure()
+# plt.imshow(phase_grad_x)
+# plt.figure()
+# plt.imshow(phase_grad_y)
+ angles = (1 - np.arctan2(phase_grad_y, phase_grad_x)/np.pi) / 2
+ phase_grad_amp = np.hypot(phase_grad_y, phase_grad_x)
+# phase_angle = (1 + np.arctan2(phase_grad_y, phase_grad_x)/np.pi) / 2
+# phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)
# Calculate the color saturation:
- saturation = np.sin(phase_magnitude/(phase_magnitude.max()+1E-30) * pi / 2)
- phase_saturation = np.dstack((saturation,)*4)
- # Color code the angle and create the holography image:
+# saturation = np.sin(phase_magnitude/(phase_magnitude.max()+1E-30) * np.pi / 2)
+# phase_saturation = np.dstack((saturation,)*4)
+ saturations = np.sin(phase_grad_amp/(phase_grad_amp.max()+1E-30)*np.pi/2) # betw. 0 and 1
+# # Color code the angle and create the holography image:
+# if grad_encode == 'dark':
+# self._log.debug('gradient encoding: dark')
+# rgba = self.HOLO_CMAP(phase_angle)
+# rgb = (255.999 * img_holo.T * saturation.T * rgba[:, :, :3].T).T.astype(np.uint8)
+# elif grad_encode == 'bright_old':
+# self._log.debug('gradient encoding: bright')
+# rgba = self.HOLO_CMAP(phase_angle)+(1-phase_saturation)*self.HOLO_CMAP_INV(phase_angle)
+# rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
+# elif grad_encode == 'bright':
+# self._log.debug('gradient encoding: bright')
+# rgba = self.HOLO_CMAP(phase_angle)+(1-phase_saturation)*self.HOLO_CMAP_INV(phase_angle)
+# rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
+# elif grad_encode == 'color':
+# self._log.debug('gradient encoding: color')
+# rgba = self.HOLO_CMAP(phase_angle)
+# rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
+# elif grad_encode == 'none':
+# self._log.debug('gradient encoding: none')
+# rgba = self.HOLO_CMAP(phase_angle)+self.HOLO_CMAP_INV(phase_angle)
+# rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
+# else:
+# raise AssertionError('Gradient encoding not recognized!')
if grad_encode == 'dark':
- self._log.debug('gradient encoding: dark')
- rgba = self.HOLO_CMAP(phase_angle)
- rgb = (255.999 * img_holo.T * saturation.T * rgba[:, :, :3].T).T.astype(np.uint8)
+ pass
elif grad_encode == 'bright':
- self._log.debug('gradient encoding: bright')
- rgba = self.HOLO_CMAP(phase_angle)+(1-phase_saturation)*self.HOLO_CMAP_INV(phase_angle)
- rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
+ saturations = 2 - saturations
elif grad_encode == 'color':
- self._log.debug('gradient encoding: color')
- rgba = self.HOLO_CMAP(phase_angle)
- rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
+ saturations = np.ones_like(saturations)
elif grad_encode == 'none':
- self._log.debug('gradient encoding: none')
- rgba = self.HOLO_CMAP(phase_angle)+self.HOLO_CMAP_INV(phase_angle)
- rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
+ saturations = 2*np.ones_like(saturations)
else:
raise AssertionError('Gradient encoding not recognized!')
+ rgb = DirectionalColormap.rgb_from_colorind_and_saturation(angles, saturations)
+ rgb = (holo.T * rgb.T).T.astype(np.uint8)
holo_image = Image.fromarray(rgb)
# If no axis is specified, a new figure is created:
if axis is None:
@@ -712,17 +749,20 @@ class PhaseMap(object):
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
# Plot the image and set axes:
- axis.imshow(holo_image, origin='lower', interpolation=interpolation)
+ axis.imshow(holo_image, origin='lower', interpolation=interpolation,
+ extent=(0, self.dim_uv[1], 0, self.dim_uv[0]))
# Set the title and the axes labels:
axis.set_title(title)
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('u-axis [px]', fontsize=15)
axis.set_ylabel('v-axis [px]', fontsize=15)
- 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))
+ if self.dim_uv[0] >= self.dim_uv[1]:
+ u_bin, v_bin = np.max((2, np.floor(9*self.dim_uv[1]/self.dim_uv[0]))), 9
+ else:
+ u_bin, v_bin = 9, np.max((2, np.floor(9*self.dim_uv[0]/self.dim_uv[1])))
+ axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
+ axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
# Return plotting axis:
return axis
@@ -776,33 +816,33 @@ class PhaseMap(object):
# Return the plotting axes:
return phase_axis, holo_axis
- @classmethod
- def make_color_wheel(cls):
- '''Display a color wheel to illustrate the color coding of the gradient direction.
-
- Parameters
- ----------
- None
-
- Returns
- -------
- None
-
- '''
- cls._log.debug('Calling make_color_wheel')
- yy, xx = np.indices((512, 512)) - 256
- r = np.sqrt(xx**2 + yy**2)
- # Create the wheel:
- color_wheel_magnitude = (1 - np.cos(r * pi/360)) / 2
- color_wheel_magnitude *= 0 * (r > 256) + 1 * (r <= 256)
- color_wheel_angle = (1 - np.arctan2(xx, -yy)/pi) / 2
- # Color code the angle and create the holography image:
- rgba = cls.HOLO_CMAP(color_wheel_angle)
- rgb = (255.999 * color_wheel_magnitude.T * rgba[:, :, :3].T).T.astype(np.uint8)
- color_wheel = Image.fromarray(rgb)
- # Plot the color wheel:
- fig = plt.figure(figsize=(4, 4))
- axis = fig.add_subplot(1, 1, 1, aspect='equal')
- axis.imshow(color_wheel, origin='lower')
- axis.xaxis.set_major_locator(NullLocator())
- axis.yaxis.set_major_locator(NullLocator())
+# @classmethod
+# def make_color_wheel(cls):
+# '''Display a color wheel to illustrate the color coding of the gradient direction.
+#
+# Parameters
+# ----------
+# None
+#
+# Returns
+# -------
+# None
+#
+# '''
+# cls._log.debug('Calling make_color_wheel')
+# yy, xx = np.indices((512, 512)) - 256
+# r = np.sqrt(xx**2 + yy**2)
+# # Create the wheel:
+# color_wheel_magnitude = (1 - np.cos(r * pi/360)) / 2
+# color_wheel_magnitude *= 0 * (r > 256) + 1 * (r <= 256)
+# color_wheel_angle = (1 - np.arctan2(xx, -yy)/pi) / 2
+# # Color code the angle and create the holography image:
+# rgba = cls.HOLO_CMAP(color_wheel_angle)
+# rgb = (255.999 * color_wheel_magnitude.T * rgba[:, :, :3].T).T.astype(np.uint8)
+# color_wheel = Image.fromarray(rgb)
+# # Plot the color wheel:
+# fig = plt.figure(figsize=(4, 4))
+# axis = fig.add_subplot(1, 1, 1, aspect='equal')
+# axis.imshow(color_wheel, origin='lower')
+# axis.xaxis.set_major_locator(NullLocator())
+# axis.yaxis.set_major_locator(NullLocator())
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 3fe0345..81f5c54 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -52,15 +52,15 @@ class PhaseMapper(object):
_log = logging.getLogger(__name__+'.PhaseMapper')
@abc.abstractmethod
- def __call__(self, mag_data):
+ def __call__(self, mag_data, offset=None):
raise NotImplementedError()
@abc.abstractmethod
- def jac_dot(self, vector):
+ def jac_dot(self, vector, offset=None):
raise NotImplementedError()
@abc.abstractmethod
- def jac_T_dot(self, vector):
+ def jac_T_dot(self, vector, offseth=None):
raise NotImplementedError()
@@ -113,7 +113,7 @@ class PhaseMapperRDFC(PhaseMapper):
self._log.debug('Calling __str__')
return 'PhaseMapperRDFC(kernel=%s)' % (self.kernel)
- def __call__(self, mag_data):
+ def __call__(self, mag_data, offset=None, ramp=None):
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.kernel.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
@@ -121,7 +121,19 @@ class PhaseMapperRDFC(PhaseMapper):
# Process input parameters:
self.u_mag[self.kernel.slice_mag] = mag_data.magnitude[0, 0, ...] # u-component
self.v_mag[self.kernel.slice_mag] = mag_data.magnitude[1, 0, ...] # v-component
- return PhaseMap(mag_data.a, self._convolve())
+ # Handle offset:
+ if offset is not None:
+ off_phase = offset * np.ones(self.kernel.dim_uv)
+ else:
+ off_phase = 0
+ # Handle ramp:
+ if ramp is not None:
+ u_ramp, v_ramp = ramp # in rad / nm
+ vv, uu = np.indices(self.kernel.dim_uv) * self.kernel.a
+ ramp_phase = u_ramp * uu + v_ramp * vv
+ else:
+ ramp_phase = 0
+ return PhaseMap(mag_data.a, self._convolve() + off_phase + ramp_phase)
def _convolve(self):
# Fourier transform the projected magnetisation:
@@ -132,7 +144,7 @@ class PhaseMapperRDFC(PhaseMapper):
# Return the result:
return fft.irfftn(self.phase_fft)[self.kernel.slice_phase]
- def jac_dot(self, vector):
+ def jac_dot(self, vector, offset=None, ramp=None):
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
@@ -152,10 +164,21 @@ class PhaseMapperRDFC(PhaseMapper):
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
self.u_mag[self.kernel.slice_mag], self.v_mag[self.kernel.slice_mag] = \
np.reshape(vector, (2,)+self.kernel.dim_uv)
- result = self._convolve().flatten()
- return result
+ # Handle offset:
+ if offset is not None:
+ off_phase = offset * np.ones(self.kernel.dim_uv)
+ else:
+ off_phase = 0
+ # Handle ramp:
+ if ramp is not None:
+ u_ramp, v_ramp = ramp # in rad / nm
+ vv, uu = np.indices(self.kernel.dim_uv) * self.kernel.a
+ ramp_phase = u_ramp * uu + v_ramp * vv
+ else:
+ ramp_phase = 0
+ return np.ravel(self._convolve() + off_phase + ramp_phase)
- def jac_T_dot(self, vector):
+ def jac_T_dot(self, vector, return_ramp=False, return_offset=False):
'''Calculate the product of the transposed Jacobi matrix with a given `vector`.
Parameters
@@ -179,8 +202,23 @@ class PhaseMapperRDFC(PhaseMapper):
v_phase_adj_fft = mag_adj_fft * -self.kernel.v_fft
u_phase_adj = fft.rfftn_adj(u_phase_adj_fft)[self.kernel.slice_phase]
v_phase_adj = fft.rfftn_adj(v_phase_adj_fft)[self.kernel.slice_phase]
- return -np.concatenate((u_phase_adj.flatten(), v_phase_adj.flatten())) # TODO: why minus?
+ # TODO: offset subtract?
+ if return_ramp:
+ offset = [vector.sum()]
+ vv, uu = np.indices(self.kernel.dim_uv) * self.kernel.a
+ u_ramp = [np.sum(vector * uu.flatten())]
+ v_ramp = [np.sum(vector * vv.flatten())]
+ result = np.concatenate((-u_phase_adj.flatten(), -v_phase_adj.flatten(),
+ offset, u_ramp, v_ramp))
+ elif return_offset:
+ offset = [vector.sum()]
+ result = np.concatenate((-u_phase_adj.flatten(), -v_phase_adj.flatten(), offset))
+ else:
+ result = np.concatenate((-u_phase_adj.flatten(), -v_phase_adj.flatten()))
+ # TODO: Why minus?
+ return result
+# TODO: ALL Docstrings with offset!!!
class PhaseMapperRDRC(PhaseMapper):
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 2e7a3ad..7438372 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -634,16 +634,16 @@ class SimpleProjector(Projector):
for row in range(size_2d)]).reshape(-1)
if dim_uv is not None:
indptr = indptr.tolist() # convert to use insert() and append()
- d_v, d_u = (dim_uv[0]-dim_v)//2, (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
+ d_v = np.floor((dim_uv[0]-dim_v)/2), np.ceil((dim_uv[0]-dim_v)/2) # padding in v
+ d_u = np.floor((dim_uv[1]-dim_u)/2), np.ceil((dim_uv[1]-dim_u)/2) # padding in u
+ indptr.extend([indptr[-1]]*d_v[1]*dim_uv[1]) # add empty lines at the end
for i in np.arange(dim_v, 0, -1): # all slices in between
- u, l = i*dim_u, (i-1)*dim_u+1 # upper / lower slice end
- indptr[u:u] = [indptr[u]] * d_u # end of the slice
- indptr[l:l] = [indptr[l]] * d_u # start of the slice
- indptr[0:0] = [0] * d_v*dim_uv[1] # insert empty rows at the beginning
+ up, lo = i*dim_u, (i-1)*dim_u # upper / lower slice end
+ indptr[up:up] = [indptr[up]] * d_u[1] # end of the slice
+ indptr[lo:lo] = [indptr[lo]] * d_u[0] # start of the slice
+ indptr = [0]*d_v[0]*dim_uv[1] + indptr # insert empty rows at the beginning
size_2d = np.prod(dim_uv) # increase size_2d
- # Make sure dim_uv is defined (used for the assertion)
- if dim_uv is None:
+ else: # Make sure dim_uv is defined (used for the assertion)
dim_uv = dim_v, dim_u
assert dim_uv[0] >= dim_v and dim_uv[1] >= dim_u, 'Projected dimensions are too small!'
# Create weight-matrix:
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 96664e8..3d2cf4a 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -48,11 +48,22 @@ def optimize_linear(costfunction, max_iter=None):
_log.info('Cost before optimization: {}'.format(costfunction(np.zeros(costfunction.n))))
x_opt = jcg.conj_grad_minimize(costfunction, max_iter=max_iter).x
_log.info('Cost after optimization: {}'.format(costfunction(x_opt)))
- # Create and return fitting MagData object:
+ # Extract additional info if necessary:
+ add_info = []
data_set = costfunction.fwd_model.data_set
+ count = data_set.count
+ if costfunction.fwd_model.fit_ramps:
+ x_opt, offsets, u_ramps, v_ramps = np.split(x_opt, [-3*count, -2*count, -count])
+ ramps = zip(u_ramps, v_ramps)
+ add_info.append(offsets)
+ add_info.append(ramps)
+ elif costfunction.fwd_model.fit_offsets:
+ x_opt, offsets = np.split(x_opt, [-count])
+ add_info.append(offsets)
+ # Create and return fitting MagData object:
mag_opt = MagData(data_set.a, np.zeros((3,) + data_set.dim))
mag_opt.set_vector(x_opt, data_set.mask)
- return mag_opt
+ return mag_opt, add_info # TODO: Docstring
def optimize_nonlin(costfunction, first_guess=None):
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index a6fac95..5d638ae 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -13,7 +13,6 @@ from scipy import sparse
import jutil.norms as jnorm
import jutil.diff as jdiff
-import jutil.operator as joperator
import logging
@@ -37,24 +36,32 @@ class Regularisator(object):
lam: float
Regularisation parameter determining the weighting between measurements and regularisation.
- '''
+ ''' # TODO: ALL dostrings with add_params of input!
__metaclass__ = abc.ABCMeta
_log = logging.getLogger(__name__+'.Regularisator')
@abc.abstractmethod
- def __init__(self, norm, lam):
+ def __init__(self, norm, lam, add_params=None):
+ # TODO: Cut off additional parameters at the end (offset and stuff)!!!
+ # OR fix x in init and cut off the rest so you don't have to care what else is needed for F
self._log.debug('Calling __init__')
self.norm = norm
self.lam = lam
+ if add_params is not None:
+ self.add_params = -add_params # Workaround! For indexing -None does NOT exist!
+ #TODO: define slice!
+ else:
+ self.add_params = None
self._log.debug('Created '+str(self))
def __call__(self, x):
self._log.debug('Calling __call__')
- return self.lam * self.norm(x)
+ # TODO: assert len(x) - self.add_params korrekt
+ return self.lam * self.norm(x[:self.add_params])
def __repr__(self):
- self._log.debug('Calling __repr__')
+ self._log.debug('Calling __repr__') # TODO: add_params
return '%s(norm=%r, lam=%r)' % (self.__class__, self.norm, self.lam)
def __str__(self):
@@ -75,7 +82,9 @@ class Regularisator(object):
Jacobi vector which represents the cost derivative of all voxels of the magnetization.
'''
- return self.lam * self.norm.jac(x)
+ result = np.zeros_like(x)
+ result[:self.add_params] = self.lam * self.norm.jac(x[:self.add_params])
+ return result
def hess_dot(self, x, vector):
'''Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
@@ -96,7 +105,9 @@ class Regularisator(object):
Product of the input `vector` with the Hessian matrix.
'''
- return self.lam * self.norm.hess_dot(x, vector)
+ result = np.zeros_like(x)
+ result[:self.add_params] = self.lam * self.norm.hess_dot(x, vector[:self.add_params])
+ return result
def hess_diag(self, x):
''' Return the diagonal of the Hessian.
@@ -116,7 +127,9 @@ class Regularisator(object):
'''
self._log.debug('Calling hess_diag')
- return self.lam * self.norm.hess_diag(x)
+ result = np.zeros_like(x)
+ result[:self.add_params] = self.lam * self.norm.hess_diag(x[:self.add_params])
+ return result
class ComboRegularisator(Regularisator):
@@ -313,14 +326,14 @@ class ZeroOrderRegularisator(Regularisator):
_log = logging.getLogger(__name__+'.ZeroOrderRegularisator')
- def __init__(self, _=None, lam=1E-4, p=2):
+ def __init__(self, _=None, lam=1E-4, p=2, add_params=None):
self._log.debug('Calling __init__')
self.p = p
if p == 2:
norm = jnorm.L2Square()
else:
norm = jnorm.LPPow(p, 1e-12)
- super(ZeroOrderRegularisator, self).__init__(norm, lam)
+ super(ZeroOrderRegularisator, self).__init__(norm, lam, add_params)
self._log.debug('Created '+str(self))
@@ -343,7 +356,7 @@ class FirstOrderRegularisator(Regularisator):
'''
- def __init__(self, mask, lam=1E-4, p=2):
+ def __init__(self, mask, lam=1E-4, p=2, add_params=None):
self.p = p
D0 = jdiff.get_diff_operator(mask, 0, 3)
D1 = jdiff.get_diff_operator(mask, 1, 3)
@@ -353,5 +366,5 @@ class FirstOrderRegularisator(Regularisator):
norm = jnorm.WeightedL2Square(D)
else:
norm = jnorm.WeightedTV(jnorm.LPPow(p, 1e-12), D, [D0.shape[0], D.shape[0]])
- super(FirstOrderRegularisator, self).__init__(norm, lam)
+ super(FirstOrderRegularisator, self).__init__(norm, lam, add_params)
self._log.debug('Created '+str(self))
diff --git a/scripts/magdata/magcreator/create_homog_slab.py b/scripts/magdata/magcreator/create_homog_slab.py
index 9fc206d..5d1394c 100644
--- a/scripts/magdata/magcreator/create_homog_slab.py
+++ b/scripts/magdata/magcreator/create_homog_slab.py
@@ -13,8 +13,8 @@ logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
# Parameters:
dim = (32, 32, 32)
a = 1.0 # in nm
-phi = np.pi/4 # in rad
-theta = 0*np.pi/2 # in rad
+phi = np.pi/2 # in rad
+theta = np.pi/2 # in rad
magnitude = 1
filename = 'magdata_mc_homog_slab.nc'
diff --git a/scripts/magdata/magdata_from_vtk.py b/scripts/magdata/magdata_from_vtk.py
index 478319d..0897ee7 100644
--- a/scripts/magdata/magdata_from_vtk.py
+++ b/scripts/magdata/magdata_from_vtk.py
@@ -6,6 +6,7 @@ import os
import imp
import numpy as np
from pylab import griddata
+from scipy.spatial import cKDTree as KDTree
import matplotlib.pyplot as plt
import vtk
from tqdm import tqdm
@@ -16,8 +17,8 @@ import logging.config
logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
###################################################################################################
-filename = 'tube_CoFeB_with_cap_3nm.vtk'
-b_0 = 1.54
+filename = 'gyroid_homx.vtk'
+b_0 = 1.
###################################################################################################
filepath = os.path.join(py.DIR_FILES, 'vtk', filename)
@@ -77,23 +78,71 @@ y = np.arange(y_cent-y_diff, y_cent+y_diff, a)
z = np.arange(z_min, z_max, a)
xx, yy = np.meshgrid(x, y)
+def excluding_mesh(x, y, nx=30, ny=30):
+ """
+ Construct a grid of points, that are some distance away from points (x,
+ """
+
+ dx = x.ptp() / nx
+ dy = y.ptp() / ny
+
+ xp, yp = np.mgrid[x.min()-2*dx:x.max()+2*dx:(nx+2)*1j,
+ y.min()-2*dy:y.max()+2*dy:(ny+2)*1j]
+ xp = xp.ravel()
+ yp = yp.ravel()
+
+ # Use KDTree to answer the question: "which point of set (x,y) is the
+ # nearest neighbors of those in (xp, yp)"
+ tree = KDTree(np.c_[x, y])
+ dist, j = tree.query(np.c_[xp, yp], k=1)
+
+ # Select points sufficiently far away
+ m = (dist > np.hypot(dx, dy))
+ return xp[m], yp[m]
+
+## Prepare fake data points
+#xp, yp = excluding_mesh(x, y, nx=35, ny=35)
+#zp = np.nan + np.zeros_like(xp)
+#
+## Grid the data plus fake data points
+#xi, yi = np.ogrid[-3:3:350j, -3:3:350j]
+#zi = griddata((np.r_[x,xp], np.r_[y,yp]), np.r_[z, zp], (xi, yi),
+# method='linear')
+#plt.imshow(zi)
+#plt.show()
+
+
# Create empty magnitude:
magnitude = np.zeros((3, len(z), len(y), len(x)), dtype=py.fft.FLOAT)
print 'Mag Dimensions:', magnitude.shape[1:]
+
+
+#import pdb; pdb.set_trace()
# Fill magnitude slice per slice:
for i, zi in tqdm(enumerate(z), total=len(z)):
- # Take all points that lie in one voxel of the new regular grid into account (use weights!):
+ # Take all points that lie in one z-voxel of the new regular grid into account (use weights!):
z_slice = data[np.abs(data[:, 2]-zi) <= a/2., :]
weights = 1 - np.abs(z_slice[:, 2]-zi)*2/a # If z is regular everywhere, weights are always 1!
+
+ # Prepare fake data points
+ x_nan, y_nan = excluding_mesh(z_slice[:, 0], z_slice[:, 1], nx=len(x)//10, ny=len(y)//10)
+ z_nan = np.nan + np.zeros_like(x_nan)
+
+
+ if i == 10:
+ axis = plt.figure().add_subplot(1, 1, 1)
+ axis.scatter(x_nan, y_nan)
+ plt.show()
+
+ grid_x = np.r_[z_slice[:, 0], x_nan]
+ grid_y = np.r_[z_slice[:, 1], y_nan]
for j in range(3): # For all 3 components!
- if condition(zi): # use insert point if condition is True!
- grid_x = np.concatenate([z_slice[:, 0], insert_x])
- grid_y = np.concatenate([z_slice[:, 1], insert_y])
- grid_z = np.concatenate([weights*z_slice[:, 3+j], np.zeros(len(insert_x))])
- else:
- grid_x = z_slice[:, 0]
- grid_y = z_slice[:, 1]
- grid_z = weights*z_slice[:, 3+j]
+# if condition(zi): # use insert point if condition is True!
+# grid_x = np.concatenate([z_slice[:, 0], insert_x])
+# grid_y = np.concatenate([z_slice[:, 1], insert_y])
+# grid_z = np.concatenate([weights*z_slice[:, 3+j], np.zeros(len(insert_x))])
+# else:
+ grid_z = np.r_[weights*z_slice[:, 3+j], z_nan]
gridded_subdata = griddata(grid_x, grid_y, grid_z, xx, yy)
magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0).astype(py.fft.FLOAT)
# Convert a to nm:
@@ -102,3 +151,5 @@ a *= 10
print 'CREATE AND SAVE MAGDATA OBJECT!'
mag_data = py.MagData(a, magnitude)
mag_data.save_to_netcdf4('magdata_vtk_{}'.format(filename.replace('.vtk', '.nc')))
+
+py.pm(mag_data).display_combined()
\ No newline at end of file
diff --git a/scripts/phasemap/phasemap_from_dm3.py b/scripts/phasemap/phasemap_from_dm3.py
index 604f3d8..f76b8ac 100644
--- a/scripts/phasemap/phasemap_from_dm3.py
+++ b/scripts/phasemap/phasemap_from_dm3.py
@@ -12,13 +12,12 @@ import logging.config
logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
###################################################################################################
-path_mag = 'zi_an_magnetite_2_particles_magnetic.dm3'
-path_ele = 'zi_an_magnetite_2_particles_electric.dm3'
-filename = 'phasemap_dm3_zi_an_magnetite_2_particles.nc'
+path_mag = 'zi_an_elongated_nanorod.dm3'
+path_ele = 'zi_an_elongated_nanorod_mip.dm3'
+filename = 'phasemap_dm3_zi_an_elongated_nanorod.nc'
a = 1.
dim_uv = None
-threshold = 2
-offset = 0.23280109
+threshold = 4.5
###################################################################################################
# Load images:
@@ -28,7 +27,7 @@ if dim_uv is not None:
im_mag = im_mag.resize(dim_uv)
im_ele = im_ele.resize(dim_uv)
# Calculate phase and mask:
-phase = np.asarray(im_mag) - offset
+phase = np.asarray(im_mag)
mask = np.where(np.asarray(im_ele) >= threshold, True, False)
# Create and save PhaseMap object:
diff --git a/scripts/phasemap/phasemap_from_image.py b/scripts/phasemap/phasemap_from_image.py
index 54c6dce..d0d1dd6 100644
--- a/scripts/phasemap/phasemap_from_image.py
+++ b/scripts/phasemap/phasemap_from_image.py
@@ -12,14 +12,14 @@ import logging.config
logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
###################################################################################################
-path_mag = 'kiyou_skyrmion.tif'
-path_ele = 'kiyou_skyrmion_mask.tif'
-filename = 'phasemap_tif_kiyou_skyrmion.nc'
-a = 1.
+path_mag = 'Arnaud_M.tif'
+path_ele = 'Arnaud_MIP_mask.tif'
+filename = 'phasemap_tif_martial_skyrmion.nc'
+a = 2
dim_uv = None
max_phase = 1
threshold = 0.5
-offset = 0.182
+offset = 0
###################################################################################################
# Load images:
diff --git a/scripts/reconstruction/reconstruction_2d_from_phasemap.py b/scripts/reconstruction/reconstruction_2d_from_phasemap.py
index c656a7b..4e17081 100644
--- a/scripts/reconstruction/reconstruction_2d_from_phasemap.py
+++ b/scripts/reconstruction/reconstruction_2d_from_phasemap.py
@@ -3,51 +3,67 @@
import numpy as np
-import pyramid as py
+import pyramid as pr
from jutil.taketime import TakeTime
import logging.config
-logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
+logging.config.fileConfig(pr.LOGGING_CONFIG, disable_existing_loggers=False)
###################################################################################################
-phase_name = 'phasemap_bmp_trevor_magnetite_1'
-b_0 = 1
+phase_name = 'phasemap_dm3_zi_an_elongated_nanorod'
+b_0 = 0.6 # in T
lam = 1E-3
max_iter = 100
-buffer_pixel = 20
+buffer_pixel = 0
+fit_ramps = True
+fit_offsets = True
###################################################################################################
# Load phase map:
-phase_map = py.PhaseMap.load_from_netcdf4(phase_name+'.nc')
-phase_map.pad(buffer_pixel, buffer_pixel)
+phase_map = pr.PhaseMap.load_from_netcdf4(phase_name+'.nc')
+phase_map.pad((buffer_pixel, buffer_pixel))
dim = (1,) + phase_map.dim_uv
-# Construct data set and regularisator:
-data = py.DataSet(phase_map.a, dim, b_0)
-data.append(phase_map, py.SimpleProjector(dim))
+# Construct regularisator, forward model and costfunction:
+data = pr.DataSet(phase_map.a, dim, b_0)
+data.append(phase_map, pr.SimpleProjector(dim))
data.set_3d_mask()
-reg = py.FirstOrderRegularisator(data.mask, lam)
+if fit_ramps:
+ add_param_count = 3
+elif fit_offsets:
+ add_param_count = 1
+else:
+ add_param_count = None
+reg = pr.FirstOrderRegularisator(data.mask, lam, add_params=add_param_count)
+fwd_model = pr.ForwardModel(data, fit_offsets=fit_offsets, fit_ramps=fit_ramps)
+cost = pr.Costfunction(fwd_model, reg)
# Reconstruct and save:
with TakeTime('reconstruction time'):
- mag_data_rec, cost = py.reconstruction.optimize_linear(data, reg, max_iter=max_iter)
+ mag_data_rec, add_info = pr.reconstruction.optimize_linear(cost, max_iter=max_iter)
+if fit_ramps:
+ offset, ramp = add_info[0][0], add_info[1][0]
+elif fit_offsets:
+ offset = add_info[0][0]
mag_data_buffer = mag_data_rec.copy()
-mag_data_rec.crop([0]*2+[buffer_pixel]*4)
+mag_data_rec.crop((0, buffer_pixel, buffer_pixel))
mag_name = '{}_lam={}'.format(phase_name.replace('phasemap', 'magdata_rec'), lam)
mag_data_rec.save_to_netcdf4(mag_name+'.nc')
# Plot stuff:
-mag_data_rec.quiver_plot('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim)/64.))
-phase_map.crop([buffer_pixel]*4)
+mag_data_rec.quiver_plot('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim)/128.))
+phase_map.crop((buffer_pixel, buffer_pixel))
phase_map.display_combined('Input Phase')
-phase_map_rec = py.pm(mag_data_rec)
+phase_map_rec = pr.pm(mag_data_rec)
phase_map_rec.mask = phase_map.mask
-phase_map_rec.display_combined('Reconstructed Phase')
+title = 'Reconstructed Phase'
+if fit_offsets:
+ print 'offset:', offset
+ title += ', fitted Offset: {:.2g} [rad]'.format(offset)
+if fit_ramps:
+ print 'ramp:', ramp
+ title += ', (Fitted Ramp: (u:{:.2g}, v:{:.2g}) [rad/nm]'.format(*ramp)
+phase_map_rec.display_combined(title)
difference = (phase_map_rec.phase-phase_map.phase).mean()
-(phase_map_rec-phase_map).display_phase('Difference (mean: {:.3g})'.format(difference))
-
-bp = buffer_pixel
-mag_data_buffer.magnitude[:, :, bp:-bp, bp:-bp] = 0
-mag_data_buffer.quiver_plot(ar_dens=np.ceil(np.max(dim)/64.))
-py.pm(mag_data_buffer).display_phase()
+(phase_map_rec-phase_map).display_phase('Difference (mean: {:.2g})'.format(difference))
diff --git a/scripts/reconstruction/reconstruction_3d_from_magdata.py b/scripts/reconstruction/reconstruction_3d_from_magdata.py
index f5cfd37..8d1a763 100644
--- a/scripts/reconstruction/reconstruction_3d_from_magdata.py
+++ b/scripts/reconstruction/reconstruction_3d_from_magdata.py
@@ -3,64 +3,92 @@
import numpy as np
-import pyramid as py
+import pyramid as pr
from jutil.taketime import TakeTime
import logging.config
-logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
+logging.config.fileConfig(pr.LOGGING_CONFIG, disable_existing_loggers=False)
###################################################################################################
-mag_name = 'magdata_mc_array_sphere_disc_slab'
+mag_name = 'magdata_mc_vortex_disc'
dim_uv = None
angles = np.linspace(-90, 90, num=7)
axes = {'x': True, 'y': True}
b_0 = 1
-noise = 0
-lam = 1E-4
+noise = 0.1
+lam = 1E-1
max_iter = 100
-save_images = False
-use_internal_mask = False
+use_internal_mask = True
+offset_max = 2
+ramp_max = 0.01
+fit_ramps = True
+fit_offsets = True
###################################################################################################
# Load magnetization distribution:
-mag_data = py.MagData.load_from_netcdf4(mag_name+'.nc')
+mag_data = pr.MagData.load_from_netcdf4(mag_name+'.nc')
dim = mag_data.dim
# Construct data set and regularisator:
-data = py.DataSet(mag_data.a, mag_data.dim, b_0)
+data = pr.DataSet(mag_data.a, mag_data.dim, b_0)
# Construct projectors:
projectors = []
for angle in angles:
angle_rad = angle*np.pi/180
if axes['x']:
- projectors.append(py.XTiltProjector(mag_data.dim, angle_rad, dim_uv))
+ projectors.append(pr.XTiltProjector(mag_data.dim, angle_rad, dim_uv))
if axes['y']:
- projectors.append(py.YTiltProjector(mag_data.dim, angle_rad, dim_uv))
+ projectors.append(pr.YTiltProjector(mag_data.dim, angle_rad, dim_uv))
# Add projectors and construct according phase maps:
data.projectors = projectors
data.phase_maps = data.create_phase_maps(mag_data)
+for i, phase_map in enumerate(data.phase_maps):
+ offset = np.random.uniform(-offset_max, offset_max)
+ ramp = np.random.uniform(-ramp_max, ramp_max), np.random.uniform(-ramp_max, ramp_max)
+ phase_map.add_ramp(offset, ramp)
+
# Add noise if necessary:
if noise != 0:
for i, phase_map in enumerate(data.phase_maps):
- phase_map += py.PhaseMap(mag_data.a, np.random.normal(0, noise, dim_uv))
+ phase_map.phase += np.random.normal(0, noise, phase_map.dim_uv)
+
data.phase_maps[i] = phase_map
+# Plot input:
+for i, phase_map in enumerate(data.phase_maps):
+ phase_map.display_phase()
+
# Construct mask:
if use_internal_mask:
data.mask = mag_data.get_mask() # Use perfect mask from mag_data!
else:
data.set_3d_mask() # Construct mask from 2D phase masks!
-# Construct regularisator:
-reg = py.FirstOrderRegularisator(data.mask, lam)
+# Construct regularisator, forward model and costfunction:
+if fit_ramps:
+ add_param_count = 3 * data.count
+elif fit_offsets:
+ add_param_count = data.count
+else:
+ add_param_count = None
+reg = pr.FirstOrderRegularisator(data.mask, lam, add_params=add_param_count)
+fwd_model = pr.ForwardModel(data, fit_offsets=fit_offsets, fit_ramps=fit_ramps)
+cost = pr.Costfunction(fwd_model, reg)
# Reconstruct and save:
with TakeTime('reconstruction time'):
- mag_data_rec, cost = py.reconstruction.optimize_linear(data, reg, max_iter=max_iter)
+ mag_data_rec, add_info = pr.reconstruction.optimize_linear(cost, max_iter=max_iter)
+if fit_ramps:
+ offset, ramp = add_info[0], add_info[1]
+ print 'offsets:', offset
+ print 'ramps:', ramp
+elif fit_offsets:
+ offset = add_info[0]
+ print 'offset:', offset
mag_name_rec = mag_name.replace('magdata', 'magdata_rec')
mag_data_rec.save_to_netcdf4(mag_name_rec+'.nc')
@@ -68,3 +96,5 @@ mag_data_rec.save_to_netcdf4(mag_name_rec+'.nc')
data.display_mask(ar_dens=np.ceil(np.max(dim)/64.))
mag_data.quiver_plot3d('Original Distribution', ar_dens=np.ceil(np.max(dim)/64.))
mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim)/64.))
+mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=np.ceil(np.max(dim)/64.),
+ coloring='amplitude')
diff --git a/scripts/temp/custom_colormap.py b/scripts/temp/custom_colormap.py
index 2ed480a..d3f32a7 100644
--- a/scripts/temp/custom_colormap.py
+++ b/scripts/temp/custom_colormap.py
@@ -170,5 +170,5 @@ Z = np.cos(X) * np.sin(Y) * 10
fig = plt.figure(figsize=(7, 7))
axis = fig.add_subplot(1, 1, 1)
-im = axis.imshow(Z, cmap=create_custom_colormap(levels=3, N=3))
+im = axis.imshow(Z, cmap=create_custom_colormap(levels=3, N=156))
fig.colorbar(im)
diff --git a/scripts/temp/test_dim_uv_projection.py b/scripts/temp/test_dim_uv_projection.py
new file mode 100644
index 0000000..b231245
--- /dev/null
+++ b/scripts/temp/test_dim_uv_projection.py
@@ -0,0 +1,24 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Jul 14 19:11:33 2015
+
+@author: Jan
+"""
+
+import numpy as np
+import pyramid as pr
+
+
+dim = (1, 2, 2)
+a = 5
+dim_uv = (a, a)
+width = (1, 1, 1)
+center = (0.5, 1, 1)
+mag_shape = pr.magcreator.Shapes.slab(dim, center, width)
+mag_data = pr.MagData(10, pr.magcreator.create_mag_dist_homog(mag_shape, np.pi/4))
+mag_data.quiver_plot()
+projector = pr.SimpleProjector(dim, axis='z', dim_uv=dim_uv)
+mag_proj = projector(mag_data)
+mag_proj.quiver_plot()
+pr.pm(mag_proj).display_combined()
+weight = np.array(projector.weight.todense())
--
GitLab