From d136112fabf661c9d0c159771e13fce640ade8b5 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Tue, 21 Jul 2015 17:56:56 +0200
Subject: [PATCH] Dedicated class for polynomial ramps of arbitrary order!
Cleanup and documentation! ramp: New class for polynomial ramps. Instantiated
in the ForwardModel. forwardmodel: Now uses new ramp class. phasemapper: Ramp
functionality moved to ramp class. reconstruction: Ramp order not used, ramp
can be now retrieved from the ramp class after
reconstruction. analytic, magcreator: Better and consistent definition of
pixel centers and their influence on the phase.
testing: Minor corrections for new syntax.
---
pyramid/__init__.py | 11 +-
pyramid/analytic.py | 28 ++-
pyramid/colormap.py | 153 +++++++++-----
pyramid/config.py | 1 +
pyramid/costfunction.py | 14 +-
pyramid/dataset.py | 34 +++-
pyramid/diagnostics.py | 2 +-
pyramid/fft.py | 3 +-
pyramid/forwardmodel.py | 117 +++++------
pyramid/kernel.py | 2 +-
pyramid/magcreator.py | 107 ++++------
pyramid/magdata.py | 70 +------
pyramid/phasemap.py | 85 +-------
pyramid/phasemapper.py | 54 +----
pyramid/projector.py | 5 +-
pyramid/quaternion.py | 3 +
pyramid/ramp.py | 190 ++++++++++++++++++
pyramid/reconstruction.py | 17 +-
pyramid/regularisator.py | 47 +++--
pyramid/version.py | 2 +-
scripts/phasemap/phasemap_from_image.py | 19 +-
.../reconstruction_2d_from_phasemap.py | 34 ++--
.../reconstruction_3d_from_magdata.py | 30 +--
scripts/temp/test_dim_uv_projection.py | 10 +-
scripts/temp/warning_control.py | 24 +++
tests/test_analytic.py | 5 +-
tests/test_magcreator.py | 28 +--
tests/test_magcreator/ref_ellipsoid.npy | Bin 584 -> 584 bytes
tests/test_regularisator.py | 2 +-
29 files changed, 586 insertions(+), 511 deletions(-)
create mode 100644 pyramid/ramp.py
create mode 100644 scripts/temp/warning_control.py
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index 9b7f840..6765e79 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -30,10 +30,14 @@ reconstruction
Reconstruct magnetic distributions from given phasemaps.
regularisator
Class to instantiate different regularisation strategies.
+ramp
+ Class which is used to add polynomial ramps to phasemaps.
diagnostics
Class to calculate diagnostics
quaternion
Class which is used for easy rotations in the Projector classes.
+colormap
+ Class which implements a custom direction encoding colormap.
fft
Class for custom FFT functions using numpy or FFTW.
@@ -46,11 +50,8 @@ 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
@@ -62,6 +63,7 @@ from .phasemap import * # analysis:ignore
from .phasemapper import * # analysis:ignore
from .projector import * # analysis:ignore
from .regularisator import * # analysis:ignore
+from .ramp import * # analysis:ignore
from .quaternion import * # analysis:ignore
from .colormap import * # analysis:ignore
from .config import * # analysis:ignore
@@ -73,7 +75,7 @@ _log = logging.getLogger(__name__)
_log.info("Starting PYRAMID V{} HG{}".format(__version__, __hg_revision__))
del logging
-__all__ = ['analytic', 'magcreator', 'reconstruction', 'fft', 'an', 'mc', 'rc']
+__all__ = ['analytic', 'magcreator', 'reconstruction', 'fft']
__all__.extend(costfunction.__all__)
__all__.extend(dataset.__all__)
__all__.extend(diagnostics.__all__)
@@ -84,5 +86,6 @@ __all__.extend(phasemap.__all__)
__all__.extend(phasemapper.__all__)
__all__.extend(projector.__all__)
__all__.extend(regularisator.__all__)
+__all__.extend(ramp.__all__)
__all__.extend(quaternion.__all__)
__all__.extend(colormap.__all__)
diff --git a/pyramid/analytic.py b/pyramid/analytic.py
index 3585c89..e7ee296 100644
--- a/pyramid/analytic.py
+++ b/pyramid/analytic.py
@@ -68,8 +68,8 @@ def phase_mag_slab(dim, a, phi, center, width, b_0=1):
+ F_0(y-y0+Ly/2, x-x0+Lx/2)))
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = a * (center[1] + 0.5) # y0, x0 define the center of a pixel,
- x0 = a * (center[2] + 0.5) # hence: (cellindex + 0.5) * grid spacing
+ y0 = a * center[1] # y0, x0 define the center of a pixel,
+ x0 = a * center[2] # hence: (cellindex + 0.5) * grid spacing
Lz, Ly, Lx = a * width[0], a * width[1], a * width[2]
coeff = - b_0 / (4*PHI_0) # Minus because of negative z-direction
# Create grid:
@@ -111,16 +111,15 @@ def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
# Function for the phase:
def phi_mag(x, y):
- r = np.hypot(x - x0, y - y0)
- r[center[1], center[2]] = 1E-30
- result = coeff * Lz * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
- result *= np.where(r <= R, 1, (R / r) ** 2)
+ r = np.hypot(x-x0, y-y0)
+ result = coeff * Lz * ((y-y0) * np.cos(phi) - (x-x0) * np.sin(phi))
+ result *= np.where(r <= R, 1, (R / (r+1E-30)) ** 2)
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
- x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
+ y0 = a * center[1]
+ x0 = a * center[2]
Lz = a * height
R = a * radius
coeff = pi * b_0 / (2*PHI_0) # Minus is gone because of negative z-direction
@@ -161,16 +160,15 @@ def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
# Function for the phase:
def phi_mag(x, y):
- r = np.hypot(x - x0, y - y0)
- r[center[1], center[2]] = 1E-30
- result = coeff * R ** 3 / r ** 2 * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
+ r = np.hypot(x-x0, y-y0)
+ result = coeff * R ** 3 / (r+1E-30) ** 2 * ((y-y0) * np.cos(phi) - (x-x0) * np.sin(phi))
result *= np.where(r > R, 1, (1 - (1 - (r / R) ** 2) ** (3. / 2.)))
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
- x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
+ y0 = a * center[1]
+ x0 = a * center[2]
R = a * radius
coeff = 2./3. * pi * b_0 / PHI_0 # Minus is gone because of negative z-direction
# Create grid:
@@ -216,8 +214,8 @@ def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
- x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
+ y0 = a * center[1]
+ x0 = a * center[2]
Lz = a * height
R = a * radius
coeff = - pi * b_0 * Lz / PHI_0 # Minus because of negative z-direction
diff --git a/pyramid/colormap.py b/pyramid/colormap.py
index 94f3e49..4f9e7b4 100644
--- a/pyramid/colormap.py
+++ b/pyramid/colormap.py
@@ -1,11 +1,9 @@
# -*- coding: utf-8 -*-
-"""
-Created on Wed Jul 08 15:43:06 2015
-
-@author: Jan
-"""
-
-# TODO: Docstrings!
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides a custom :class:`~.DirectionalColormap` colormap class which has a few
+additional functions and can encode three-dimensional directions."""
import numpy as np
@@ -21,6 +19,28 @@ __all__ = ['DirectionalColormap']
class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
+ '''Colormap subclass for encoding 3D-directions with colors..
+
+ This class is a subclass of the :class:`~matplotlib.pyplot.colors.LinearSegmentedColormap`
+ class with a few classmethods which can be used for convenience. The
+ :method:`.display_colorwheel` method can be used to display a colorhweel which is used for
+ the directional encoding and three `rgb_from_...` methods are used to calculate rgb tuples
+ from 3D-directions, angles or directly from a colorindex and a saturation value. This is
+ useful for quiverplots where the arrow colors can be set individually by providing said rgb-
+ tuples. Arrows in plane will be encoded with full color saturation and arrow pointing down
+ will gradually lose saturation until they are black when pointing down. Arrows pointing up
+ will 'oversaturate' (the saturation value will go from 1 up to 2), which in consequence will
+ partially add the inverse colormap to the arrows turning them white if they point up (rgb-
+ tuple: 255, 255, 255).
+
+ Attributes
+ ----------
+ inverted: boolean, optional
+ Flag which is used to invert the internal colormap (default is False).
+ Just influences the classical use as a normal colormap, not the classmethods!
+
+ '''
+
_log = logging.getLogger(__name__+'.DirectionalColormap')
CDICT = {'red': [(0.00, 0.0, 0.0),
@@ -62,43 +82,8 @@ class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
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))
@@ -109,13 +94,15 @@ class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
Parameters
----------
- None
+ mode : {'white_to_color', 'color_to_black', 'black_to_color', 'white_to_color_to_black'}
+ Optional string for determining the color scheme of the color wheel. Describes the
+ order of colors from the center to the outline.
Returns
-------
None
- ''' # TODO: mode docstring!
+ '''
cls._log.debug('Calling display_color_wheel')
yy, xx = np.indices((512, 512)) - 256
r = np.hypot(xx, yy)
@@ -123,7 +110,7 @@ class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
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
+ pass # no further modification necessary!
elif mode == 'color_to_black':
saturation = 1 - saturation
elif mode == 'white_to_color':
@@ -132,7 +119,7 @@ class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
saturation = 2 - 2*saturation
else:
raise Exception('Invalid color wheel mode!')
- saturation *= (r <= 256) # TODO: [r<=256]
+ saturation *= np.where(r <= 256, 1, 0) # Cut out the wheel!
rgb = cls.rgb_from_colorind_and_saturation(colorind, saturation)
color_wheel = Image.fromarray(rgb)
# Plot the color wheel:
@@ -142,21 +129,89 @@ class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
plt.tick_params(axis='both', which='both', labelleft='off', labelbottom='off',
left='off', right='off', top='off', bottom='off')
- @classmethod # TODO: Docstrings!
+ @classmethod
def rgb_from_colorind_and_saturation(cls, colorind, saturation):
- c, s = colorind, saturation
+ '''Construct a rgb tuple from a colorindex and a saturation value.
+
+ Parameters
+ ----------
+ colorind : float, [0, 1)
+ Colorindex specifying the directional color according to the CDICT colormap.
+ The colormap is periodic so that a value of 1 is equivalent to 0 again.
+ saturation : [0, 2]float, optional
+ Saturation value for the color. The lower hemisphere uses values from 0 to 1 in a
+ traditional sense of saturation with no saturation for directions pointing down (black)
+ and full saturation in plane (full colors). Higher values (between 1 and 2) add the
+ inverse colormap described in CDICT_INV to gradually increase the complementary colors
+ so that arrows pointing up appear white.
+ Azimuthal angle of the desired direction to encode (in rad). Default: pi/2 (in-plane).
+
+ Returns
+ -------
+ rgb : tuple (N=3)
+ rgb tuple with the encoded direction color.
+
+ Notes
+ -----
+ Also works with numpy arrays as input! Always returns array of shape (..., 3)!
+
+ '''
+ cls._log.debug('Calling rgb_from_colorind_and_saturation')
+ c, s = np.ravel(colorind), np.ravel(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):
+ '''Construct a rgb tuple from two angles representing a 3D direction.
+
+ Parameters
+ ----------
+ phi : float
+ Polar angle of the desired direction to encode (in rad).
+ theta : float, optional
+ Azimuthal angle of the desired direction to encode (in rad). Default: pi/2 (in-plane).
+
+ Returns
+ -------
+ rgb : tuple (N=3)
+ rgb tuple with the encoded direction color.
+
+ Notes
+ -----
+ Also works with numpy arrays as input!
+
+ '''
+ cls._log.debug('Calling rgb_from_angles')
colorind = (1 + phi/np.pi) / 2
- saturation = 2 * (1 - theta/np.pi) # goes from 0 to 2 # TODO: explain!
+ saturation = 2 * (1 - theta/np.pi) # goes from 0 to 2!
return cls.rgb_from_colorind_and_saturation(colorind, saturation)
@classmethod
def rgb_from_direction(cls, x, y, z):
+ '''Construct a rgb tuple from three coordinates representing a 3D direction.
+
+ Parameters
+ ----------
+ x : float
+ x-coordinate of the desired direction to encode.
+ y : float
+ y-coordinate of the desired direction to encode.
+ z : float
+ z-coordinate of the desired direction to encode.
+
+ Returns
+ -------
+ rgb : tuple (N=3)
+ rgb tuple with the encoded direction color.
+
+ Notes
+ -----
+ Also works with numpy arrays as input!
+
+ '''
+ cls._log.debug('Calling rgb_from_direction')
phi = np.arctan2(y, x)
r = np.sqrt(x**2 + y**2 + z**2)
theta = np.arccos(z / (r+1E-30))
diff --git a/pyramid/config.py b/pyramid/config.py
index 32122b8..3b49a4c 100644
--- a/pyramid/config.py
+++ b/pyramid/config.py
@@ -4,6 +4,7 @@
#
"""Config file for the pyramid package."""
+
import os
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 03ddcde..e60f50a 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -17,6 +17,7 @@ __all__ = ['Costfunction']
class Costfunction(object):
+
'''Class for calculating the cost of a 3D magnetic distributions in relation to 2D phase maps.
Represents a strategy for the calculation of the `cost` of a 3D magnetic distribution in
@@ -55,13 +56,12 @@ class Costfunction(object):
if self.regularisator is None:
self.regularisator = NoneRegularisator()
# Extract important information:
- data_set = fwd_model.data_set
- self.y = data_set.phase_vec
- 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
+ self.y = self.fwd_model.data_set.phase_vec
+ self.n = self.fwd_model.n
+ self.m = self.fwd_model.m
+ if self.fwd_model.data_set.Se_inv is None:
+ self.fwd_model.data_set.set_Se_inv_diag_with_conf()
+ self.Se_inv = self.fwd_model.data_set.Se_inv
self._log.debug('Created '+str(self))
def __repr__(self):
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index a31802c..4262212 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -8,9 +8,7 @@ and additional data like corresponding projectors."""
import numpy as np
from numbers import Number
-
from scipy import sparse
-
import matplotlib.pyplot as plt
from pyramid.phasemap import PhaseMap
@@ -265,7 +263,7 @@ class DataSet(object):
'''
self._log.debug('Calling display_mask')
if self.mask is not None:
- from mayavi import mlab
+ from mayavi import mlab # TODO: Supress annoying warning from traits!
zz, yy, xx = np.indices(self.dim)
ad = ar_dens
zz = zz[::ad, ::ad, ::ad].flatten()
@@ -358,7 +356,31 @@ class DataSet(object):
phase_maps = self.create_phase_maps(mag_data)
else:
phase_maps = self.phase_maps
- [phase_map.display_combined('{} ({})'.format(title, self.projectors[i].get_info()),
- cmap, limit, norm, gain, interpolation, grad_encode)
- for (i, phase_map) in enumerate(phase_maps)]
+ for (i, phase_map) in enumerate(phase_maps):
+ phase_map.display_combined('{} ({})'.format(title, self.projectors[i].get_info()),
+ cmap, limit, norm, gain, interpolation, grad_encode)
plt.show()
+
+
+# TODO: Multiprocessing! ##########################################################################
+# TODO: Use proxy objects? (https://docs.python.org/2/library/multiprocessing.html#proxy-objects)
+# class DistributedDataSet(DataSet):
+#
+# @property
+# def count(self):
+# return np.sum([len(data_set.projectors) for data_set in self.data_sets])
+#
+# def __init__(self, a, dim, b_0=1, mask=None, Se_inv=None):
+# # TODO: get number of processes!
+# self.nprocs = 4
+# self.data_sets = []
+# for proc_ind in range(self.nprocs):
+# # TODO: Create processes and let DataSet live there!
+# self.data_sets.append(DataSet(a, dim, b_0, mask, Se_inv))
+#
+# def __get_item__(self, key):
+# return self.data_sets[key]
+#
+# def append(self, phase_map, projector):
+# self.data_sets[self.count % self.nprocs].append(phase_map, projector)
+###################################################################################################
diff --git a/pyramid/diagnostics.py b/pyramid/diagnostics.py
index 5930052..7aa2b01 100644
--- a/pyramid/diagnostics.py
+++ b/pyramid/diagnostics.py
@@ -10,7 +10,6 @@ import numpy as np
import matplotlib.pyplot as plt
import jutil
-
from pyramid import fft
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
@@ -199,6 +198,7 @@ class Diagnostics(object):
x = mag_x.sum(axis=(0, 1))
y = mag_y.sum(axis=(0, 2))
z = mag_z.sum(axis=(1, 2))
+ # Plot helpful stuff:
plt.figure()
plt.axhline(y=0, ls='-', color='k')
plt.axhline(y=1, ls='-', color='k')
diff --git a/pyramid/fft.py b/pyramid/fft.py
index f65df71..d7d10f2 100644
--- a/pyramid/fft.py
+++ b/pyramid/fft.py
@@ -16,9 +16,10 @@ import numpy as np
import cPickle as pickle
import os
+from pyramid.config import NTHREADS
+
import logging
-from pyramid.config import NTHREADS
# pyFFTW depends on this
try:
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 24be8ad..da4b139 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -9,6 +9,7 @@ threedimensional magnetization distribution onto a two-dimensional phase map."""
import numpy as np
from pyramid.magdata import MagData
+from pyramid.ramp import Ramp
import logging
@@ -21,22 +22,20 @@ class ForwardModel(object):
'''Class for mapping 3D magnetic distributions to 2D phase maps.
Represents a strategy for the mapping of a 3D magnetic distribution to two-dimensional
- phase maps. Can handle a list of `projectors` of :class:`~.Projector` objects, which describe
- different projection angles, so many phase_maps can be created from one magnetic distribution.
- All required data should be given in a :class:`~DataSet` object.
+ phase maps. A :class:`~.DataSet` object is given which is used as input for the model
+ (projectors, phase_mappers, etc.). A `ramp_order` can be specified to add polynomial ramps
+ to the constructed phase maps (which can also be reconstructed!). A :class:`~.Ramp` class
+ object will be constructed accordingly, which also holds all info about the ramps after a
+ reconstruction.
Attributes
----------
data_set: :class:`~dataset.DataSet`
:class:`~dataset.DataSet` object, which stores all required information calculation.
- projectors : list of :class:`~.Projector`
- A list of all :class:`~.Projector` objects representing the projection directions.
- kernel : :class:`~.Kernel`
- A kernel which describes the phasemapping of the 2D projected magnetization distribution.
- a : float
- The grid spacing in nm.
- dim : tuple (N=3)
- Dimensions of the 3D magnetic distribution.
+ ramp_order : int or None (default)
+ Polynomial order of the additional phase ramp which will be added to the phase maps.
+ All ramp parameters have to be at the end of the input vector and are split automatically.
+ Default is None (no ramps are added).
m: int
Size of the image space. Number of pixels of the 2-dimensional projected grid.
n: int
@@ -46,22 +45,19 @@ class ForwardModel(object):
_log = logging.getLogger(__name__+'.ForwardModel')
- def __init__(self, data_set, fit_ramps=False, fit_offsets=False):
+ def __init__(self, data_set, ramp_order=None):
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 + fit_offsets * data_set.count + fit_ramps * 2 * data_set.count
- # TODO: bools nicht als integer verwenden!
+ self.ramp_order = ramp_order
+ self.phase_mappers = self.data_set.phase_mappers
+ self.ramp = Ramp(self.data_set, self.ramp_order)
+ self.m = self.data_set.m
+ self.n = self.data_set.n
+ if self.ramp.n is not None: # Additional parameters have to be fitted!
+ self.n += self.ramp.n
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))
+ self.mag_data = MagData(self.data_set.a, np.zeros((3,)+self.data_set.dim))
self._log.debug('Creating '+str(self))
# TODO: Multiprocessing! ##########################################################################
# nprocs = 4
@@ -89,27 +85,22 @@ 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])
+ # Extract ramp parameters if necessary (x will be shortened!):
+ x = self.ramp.extract_ramp_params(x)
+ # Reset mag_data and fill with vector:
self.mag_data.magnitude[...] = 0
self.mag_data.set_vector(x, self.data_set.mask)
+ # Simulate all phase maps and create result vector:
result = np.zeros(self.m)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
mapper = self.phase_mappers[projector.dim_uv]
- phase_map = mapper(projector(self.mag_data), offsets[i], ramps[i])
+ phase_map = mapper(projector(self.mag_data))
+ phase_map += self.ramp(i) # add ramp!
result[hp[i]:hp[i+1]] = phase_map.phase_vec
return np.reshape(result, -1)
-###################################################################################################
+# TODO: Multiprocessing! ##########################################################################
# nprocs = 4
# # Set up processes:
# for i, projector in enumerate(self.data_set.projectors):
@@ -135,7 +126,8 @@ class ForwardModel(object):
implemented for the case that in the future nonlinear problems have to be solved.
vector : :class:`~numpy.ndarray` (N=1)
Vectorized form of the 3D magnetization distribution. First the `x`, then the `y` and
- lastly the `z` components are listed.
+ lastly the `z` components are listed. Ramp parameters are also added at the end if
+ necessary.
Returns
-------
@@ -144,22 +136,19 @@ 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])
+ # Extract ramp parameters if necessary (vector will be shortened!):
+ vector = self.ramp.extract_ramp_params(vector)
+ # Reset mag_data and fill with vector:
self.mag_data.magnitude[...] = 0
self.mag_data.set_vector(vector, self.data_set.mask)
+ # Simulate all phase maps and create result vector:
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
mapper = self.phase_mappers[projector.dim_uv]
- res = mapper.jac_dot(projector.jac_dot(mag_vec), offsets[i], ramps[i])
+ res = mapper.jac_dot(projector.jac_dot(mag_vec))
+ res += self.ramp.jac_dot(i) # add ramp!
result[hp[i]:hp[i+1]] = res
return result
@@ -182,33 +171,27 @@ class ForwardModel(object):
-------
result_vector : :class:`~numpy.ndarray` (N=1)
Product of the transposed Jacobi matrix (which is not explicitely calculated) with
- the input `vector`.
+ the input `vector`. If necessary, transposed ramp parameters are concatenated.
'''
- 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]]
+ sub_vec = vector[hp[i]:hp[i+1]]
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)
+ proj_T_result += projector.jac_T_dot(mapper.jac_T_dot(sub_vec))
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
+ ramp_params = self.ramp.jac_T_dot(vector) # calculate ramp_params separately!
+ return np.concatenate((result, ramp_params))
+
+# TODO: Multiprocessing! ##########################################################################
+#class DistributedForwardModel(ForwardModel):
+#
+# def __init__(self, distributed_data_set, ramp_order=None):
+# self.nprocs = distributed_data_set.nprocs
+# self.fwd_models = []
+# for proc_ind in range(self.nprocs):
+# data_set = distributed_data_set[proc_ind]
+# self.fwd_models.append(ForwardModel(data_set))
+###################################################################################################
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 724514f..317ce38 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -123,7 +123,7 @@ class Kernel(object):
def _get_elementary_phase(self, geometry, n, m, a):
self._log.debug('Calling _get_elementary_phase')
if geometry == 'disc':
- in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
+ in_or_out = ~ np.logical_and(n == 0, m == 0)
return m / (n**2 + m**2 + 1E-30) * in_or_out
elif geometry == 'slab':
def F_a(n, m):
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index a777709..e42d097 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -21,7 +21,6 @@ from __future__ import division
import numpy as np
from numpy import pi
-
import abc
import logging
@@ -72,14 +71,7 @@ class Shapes(object):
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
+ return np.logical_and.reduce((xx_shape, yy_shape, zz_shape))
@classmethod
def disc(cls, dim, center, radius, height, axis='z'):
@@ -110,25 +102,18 @@ class Shapes(object):
assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
+ zz, yy, xx = np.indices(dim) + 0.5
+ xx -= center[2]
+ yy -= center[1]
+ zz -= center[0]
if axis == 'z':
- mag_shape = np.array([[[np.hypot(x - center[2], y - center[1]) <= radius
- and abs(z - center[0]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ uu, vv, ww = xx, yy, zz
elif axis == 'y':
- mag_shape = np.array([[[np.hypot(x - center[2], z - center[0]) <= radius
- and abs(y - center[1]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ uu, vv, ww = zz, xx, yy
elif axis == 'x':
- mag_shape = np.array([[[np.hypot(y - center[1], z - center[0]) <= radius
- and abs(x - center[2]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
- return mag_shape
+ uu, vv, ww = yy, zz, xx
+ return np.logical_and(np.where(np.hypot(uu, vv) <= radius, True, False),
+ np.where(abs(ww) <= height/2, True, False))
@classmethod
def ellipse(cls, dim, center, width, height, axis='z'):
@@ -159,28 +144,19 @@ class Shapes(object):
assert np.shape(width) == (2,), 'Parameter width has to be a a tuple of length 2!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
+ zz, yy, xx = np.indices(dim) + 0.5
+ xx -= center[2]
+ yy -= center[1]
+ zz -= center[0]
if axis == 'z':
- mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
- (y-center[1])/(width[0]/2.)) <= 1
- and abs(z - center[0]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ uu, vv, ww = xx, yy, zz
elif axis == 'y':
- mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
- (z-center[0])/(width[0]/2.)) <= 1
- and abs(y-center[1]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ uu, vv, ww = xx, zz, yy
elif axis == 'x':
- mag_shape = np.array([[[np.hypot((y-center[1])/(width[1]/2.),
- (z-center[0])/(width[0]/2.)) <= 1
- and abs(x - center[2]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
- return mag_shape
+ uu, vv, ww = yy, zz, xx
+ distance = np.hypot(uu/(width[1]/2), vv/(width[0]/2))
+ return np.logical_and(np.where(distance <= 1, True, False),
+ np.where(abs(ww) <= height/2, True, False))
@classmethod
def sphere(cls, dim, center, radius):
@@ -205,13 +181,9 @@ class Shapes(object):
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
- mag_shape = np.array([[[np.sqrt((x-center[2])**2
- + (y-center[1])**2
- + (z-center[0])**2) <= radius
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
- return mag_shape
+ zz, yy, xx = np.indices(dim) + 0.5
+ distance = np.sqrt((xx-center[2])**2 + (yy-center[1])**2 + (zz-center[0])**2)
+ return np.where(distance <= radius, True, False)
@classmethod
def ellipsoid(cls, dim, center, width):
@@ -236,13 +208,11 @@ class Shapes(object):
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a a tuple of length 3!'
- mag_shape = np.array([[[(x-center[2])**2/(width[2]/2)**2
- + (y-center[1])**2/(width[1]/2)**2
- + (z-center[0])**2/(width[0]/2)**2 <= 1
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
- return mag_shape
+ zz, yy, xx = np.indices(dim) + 0.5
+ distance = np.sqrt(((xx-center[2]) / (width[2]/2))**2
+ + ((yy-center[1]) / (width[1]/2))**2
+ + ((zz-center[0]) / (width[0]/2))**2)
+ return np.where(distance <= 1, True, False)
@classmethod
def filament(cls, dim, pos, axis='y'):
@@ -350,12 +320,19 @@ def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
the vortex direction (left-handed instead of right-handed).
axis : {'z', 'y', 'x'}, optional
The orientation of the vortex axis. The default is 'z'.
+
Returns
-------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The magnetic distribution as a tuple of the 3 components in
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
+ Notes
+ -----
+ To avoid singularities, the vortex center should lie between the pixel centers (which
+ reside at coordinates with _.5 at the end), i.e. integer values should be used as center
+ coordinates (e.g. coordinate 1 lies between the first and the second pixel).
+
'''
_log.debug('Calling create_mag_dist_vortex')
dim = np.shape(mag_shape)
@@ -364,12 +341,12 @@ def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
assert center is None or len(center) in {2, 3}, \
'Vortex center has to be defined in 3D or 2D or not at all!'
if center is None:
- center = (int(dim[1]/2)-0.5, int(dim[2]/2)-0.5) # center has to be between (!) two pixels
+ center = (dim[1]/2, dim[2]/2)
if axis == 'z':
if len(center) == 3: # if a 3D-center is given, just take the x and y components
center = (center[1], center[2])
- u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
- v = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
+ u = np.linspace(-center[1], dim[2]-1-center[1], dim[2]) + 0.5 # pixel center!
+ v = np.linspace(-center[0], dim[1]-1-center[0], dim[1]) + 0.5 # pixel center!
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=0)
phi = np.tile(phi, (dim[0], 1, 1))
@@ -379,8 +356,8 @@ def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
elif axis == 'y':
if len(center) == 3: # if a 3D-center is given, just take the x and z components
center = (center[0], center[2])
- u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
- v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
+ u = np.linspace(-center[1], dim[2]-1-center[1], dim[2]) + 0.5 # pixel center!
+ v = np.linspace(-center[0], dim[0]-1-center[0], dim[0]) + 0.5 # pixel center!
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=1)
phi = np.tile(phi, (1, dim[1], 1))
@@ -390,8 +367,8 @@ def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
elif axis == 'x':
if len(center) == 3: # if a 3D-center is given, just take the z and y components
center = (center[0], center[1])
- u = np.linspace(-center[1], dim[1]-1-center[1], dim[1])
- v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
+ u = np.linspace(-center[1], dim[1]-1-center[1], dim[1]) + 0.5 # pixel center!
+ v = np.linspace(-center[0], dim[0]-1-center[0], dim[0]) + 0.5 # pixel center!
uu, vv = np.meshgrid(u, v)
phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=2)
phi = np.tile(phi, (1, 1, dim[2]))
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 04edc61..88473d8 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -8,12 +8,11 @@
from __future__ import division
import os
-
import numpy as np
from numpy.linalg import norm
from scipy.ndimage.interpolation import zoom
-
-import matplotlib as mpl
+from numbers import Number
+import netCDF4
import matplotlib.pyplot as plt
import matplotlib.cm as cmx
from matplotlib.ticker import MaxNLocator
@@ -21,10 +20,6 @@ from matplotlib.ticker import MaxNLocator
from pyramid import fft
from pyramid.colormap import DirectionalColormap
-from numbers import Number
-
-import netCDF4
-
import logging
@@ -168,45 +163,6 @@ 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)
-
def copy(self):
'''Returns a copy of the :class:`~.MagData` object
@@ -695,6 +651,8 @@ class MagData(object):
scale: float, optional
Additional multiplicative factor scaling the arrow length. Default is 1
(no further scaling).
+ show_mask: boolean
+ Default is True. Shows the outlines of the mask slice if available.
Returns
-------
@@ -702,7 +660,6 @@ 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).'
@@ -832,8 +789,7 @@ class MagData(object):
'''
self._log.debug('Calling quiver_plot3D')
- from mayavi import mlab
- import warnings
+ from mayavi import mlab # TODO: Supress annoying warning from traits!
a = self.a
dim = self.dim
if limit is None:
@@ -849,6 +805,7 @@ class MagData(object):
z_mag = self.magnitude[2][::ad, ::ad, ::ad].flatten()
# Plot them as vectors:
mlab.figure(size=(750, 700))
+ import warnings
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!
@@ -858,19 +815,8 @@ class MagData(object):
if coloring == 'angle': # Encodes the full angle via colorwheel and saturation
self._log.debug('Encoding full 3D angles')
from tvtk.api import tvtk
- 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)
+ rgb = DirectionalColormap.rgb_from_direction(x_mag, y_mag, z_mag)
+ colors = map(tuple, rgb) # convert to list of tuples!
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 033cf43..4547745 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -6,21 +6,17 @@
import os
-
import numpy as np
from scipy.ndimage.interpolation import zoom
-
-from pyramid.colormap import DirectionalColormap
-
+from numbers import Number
+import netCDF4
+from PIL import Image
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import MaxNLocator, FuncFormatter
-from PIL import Image
-from numbers import Number
-
-import netCDF4
+from pyramid.colormap import DirectionalColormap
import logging
@@ -219,7 +215,6 @@ 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__')
@@ -256,12 +251,6 @@ 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
@@ -691,45 +680,13 @@ class PhaseMap(object):
if title is None:
title = 'Contour Map (gain: %.2g)' % gain
# Calculate the holography image intensity:
-# 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_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) * 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!')
+ # Calculate color encoding:
if grad_encode == 'dark':
pass
elif grad_encode == 'bright':
@@ -740,6 +697,7 @@ class PhaseMap(object):
saturations = 2*np.ones_like(saturations)
else:
raise AssertionError('Gradient encoding not recognized!')
+ # Calculate colored holo image:
rgb = DirectionalColormap.rgb_from_colorind_and_saturation(angles, saturations)
rgb = (holo.T * rgb.T).T.astype(np.uint8)
holo_image = Image.fromarray(rgb)
@@ -815,34 +773,3 @@ class PhaseMap(object):
self.display_phase(cmap='RdBu', limit=limit, norm=norm, axis=phase_axis)
# 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())
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 81f5c54..a4c63d5 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -5,14 +5,11 @@
"""This module executes several forward models to calculate the magnetic or electric phase map from
a given projection of a 3-dimensional magnetic distribution (see :mod:`~pyramid.projector`).
For the magnetic phase map, an approach using real space and one using Fourier space is provided.
-The electrostatic contribution is calculated by using the assumption of a mean inner potential.
-
-"""
+The electrostatic contribution is calculated by using the assumption of a mean inner potential."""
import numpy as np
from numpy import pi
-
import abc
import pyramid.numcore.phasemapper_core as nc
@@ -113,7 +110,7 @@ class PhaseMapperRDFC(PhaseMapper):
self._log.debug('Calling __str__')
return 'PhaseMapperRDFC(kernel=%s)' % (self.kernel)
- def __call__(self, mag_data, offset=None, ramp=None):
+ def __call__(self, mag_data):
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,19 +118,7 @@ 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
- # 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)
+ return PhaseMap(mag_data.a, self._convolve())
def _convolve(self):
# Fourier transform the projected magnetisation:
@@ -144,7 +129,7 @@ class PhaseMapperRDFC(PhaseMapper):
# Return the result:
return fft.irfftn(self.phase_fft)[self.kernel.slice_phase]
- def jac_dot(self, vector, offset=None, ramp=None):
+ def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
@@ -164,21 +149,9 @@ 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)
- # 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)
+ return np.ravel(self._convolve())
- def jac_T_dot(self, vector, return_ramp=False, return_offset=False):
+ def jac_T_dot(self, vector):
'''Calculate the product of the transposed Jacobi matrix with a given `vector`.
Parameters
@@ -202,23 +175,10 @@ 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]
- # 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()))
+ 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 7438372..937fc0a 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -10,12 +10,9 @@ from __future__ import division
import numpy as np
from numpy import pi
-
-import abc
-
import itertools
-
from scipy.sparse import coo_matrix, csr_matrix
+import abc
from pyramid.magdata import MagData
from pyramid.quaternion import Quaternion
diff --git a/pyramid/quaternion.py b/pyramid/quaternion.py
index 9a9a431..f23291c 100644
--- a/pyramid/quaternion.py
+++ b/pyramid/quaternion.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""This module provides the :class:`~.Quaternion` class which can be used for rotations."""
diff --git a/pyramid/ramp.py b/pyramid/ramp.py
new file mode 100644
index 0000000..c12b0a0
--- /dev/null
+++ b/pyramid/ramp.py
@@ -0,0 +1,190 @@
+# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides the :class:`~.Ramp` class which implements polynomial phase ramps."""
+
+
+import numpy as np
+
+from pyramid.phasemap import PhaseMap
+
+
+__all__ = ['Ramp']
+
+
+class Ramp(object):
+
+ '''Class representing a polynomial phase ramp.
+
+ Sometimes additional phase ramps occur in phase maps which do not stem from a magnetization
+ distribution inside the FOV. This class allows the construction (and via the derivative
+ functions also the reconstruction) of a polynomial ramp. This class is generally constructed
+ within the ForwardModel and can be retrieved as its attribute if ramp information should be
+ accessed.
+
+ Attributes
+ ----------
+ data_set : :class:`~dataset.DataSet`
+ :class:`~dataset.DataSet` object, which stores all required information calculation.
+ order : int or None (default)
+ Polynomial order of the additional phase ramp which will be added to the phase maps.
+ All ramp parameters have to be at the end of the input vector and are split automatically.
+ Default is None (no ramps are added).
+ deg_of_freedom : int
+ Number of degrees of freedom. This is calculated to ``1 + 2 * order``. There is just one
+ degree of freedom for a ramp of order zero (offset), every higher order contributes two
+ degrees of freedom.
+ param_cache : :class:`numpy.ndarray` (N=2)
+ Parameter cache which is used to store the polynomial coefficients. Higher coefficients
+ (one for each degree of freedom) are saved along the first axis, values for different
+ images along the second axis.
+ n : int
+ Size of the input space. Coincides with the numer of entries in `param_cache` and
+ calculates to ``deg_of_freedom * data_set.count``.
+
+ Notes
+ -----
+ After a reconstruction the relevant polynomial ramp information is stored in the
+ `param_cache`. If a phasemap with index `i` in the DataSet should be corrected use:
+
+ >>> phase_map -= ramp(i, dof_list)
+
+ The optional parameter `dof_list` can be used to specify a list of degrees of freedom which
+ should be used for the ramp (e.g. `[0]` will just apply the offset, `[0, 1, 2]` will apply
+ the offset and linear ramps in both directions).
+
+ Fitting polynoms of higher orders than `order = 1` is possible but not recommended, because
+ features which stem from the magnetization could be covered by the polynom, decreasing the
+ phase contribution of the magnetization distribution, leading to a false retrieval.
+
+ '''
+
+ def __init__(self, data_set, order=None):
+ self.data_set = data_set
+ assert order is None or (isinstance(order, int) and order >= 0), \
+ 'Order has to be None or a positive integer!'
+ self.order = order
+ self.deg_of_freedom = (1 + 2 * self.order) if self.order is not None else 0
+ self.param_cache = np.zeros((self.deg_of_freedom, self.data_set.count))
+ self.n = self.deg_of_freedom * self.data_set.count
+
+ def __call__(self, index, dof_list=None):
+ if self.order is None: # Do nothing if order is None!
+ return 0
+ else:
+ if dof_list is None: # if no specific list is supplied!
+ dof_list = range(self.deg_of_freedom) # use all available degrees of freedom
+ dim_uv = self.data_set.projectors[index].dim_uv
+ phase_ramp = np.zeros(dim_uv)
+ # Iterate over all degrees of freedom:
+ for dof in dof_list:
+ # Add the contribution of the current degree of freedom:
+ phase_ramp += self.param_cache[dof][index] * self.create_poly_mesh(dof, dim_uv)
+ return PhaseMap(self.data_set.a, phase_ramp, mask=np.zeros(dim_uv, dtype=np.bool))
+
+ def jac_dot(self, index):
+ '''Calculate the product of the Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the 3D magnetization distribution. First the `x`, then the `y` and
+ lastly the `z` components are listed. Ramp parameters are also added at the end if
+ necessary.
+
+ Returns
+ -------
+ result_vector : :class:`~numpy.ndarray` (N=1)
+ Product of the Jacobi matrix (which is not explicitely calculated) with the input
+ `vector`. Just the ramp contribution is calculated!
+
+ '''
+ if self.order is None: # Do nothing if order is None!
+ return 0
+ else:
+ dim_uv = self.data_set.projectors[index].dim_uv
+ phase_ramp = np.zeros(dim_uv)
+ # Iterate over all degrees of freedom:
+ for dof in range(self.deg_of_freedom):
+ # Add the contribution of the current degree of freedom:
+ phase_ramp += self.param_cache[dof][index] * self.create_poly_mesh(dof, dim_uv)
+ return np.ravel(phase_ramp)
+
+ def jac_T_dot(self, vector):
+ ''''Calculate the transposed ramp parameters from a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of all 2D phase maps one after another in one vector.
+
+ Returns
+ -------
+ result_vector : :class:`~numpy.ndarray` (N=1)
+ Transposed ramp parameters.
+
+ '''
+ result = []
+ hp = self.data_set.hook_points
+ # Iterate over all degrees of freedom:
+ for dof in range(self.deg_of_freedom):
+ # Iterate over all projectors:
+ for i, projector in enumerate(self.data_set.projectors):
+ sub_vec = vector[hp[i]:hp[i+1]]
+ poly_mesh = self.create_poly_mesh(dof, projector.dim_uv)
+ # Transposed ramp parameters: summed product of the vector with the poly-mesh:
+ result.append(np.sum(sub_vec * np.ravel(poly_mesh)))
+ return result
+
+ def create_poly_mesh(self, deg_of_freedom, dim_uv):
+ '''Create a polynomial mesh for the ramp calculation for a specific degree of freedom.
+
+ Parameters
+ ----------
+ deg_of_freedom : int
+ Current degree of freedom for which the mesh should be created. 0 corresponds to a
+ simple offset, 1 corresponds to a linear ramp in u-direction, 2 to a linear ramp in
+ v-direction and so on.
+ dim_uv : tuple (N=2)
+ Dimensions of the 2D mesh that should be created.
+
+ Returns
+ -------
+ result_mesh : :class:`~numpy.ndarray` (N=2)
+ Polynomial mesh that was created and can be used for further calculations.
+
+ '''
+ # Determine if u-direction (u_or_v == 1) or v-direction (u_or_v == 0)!
+ u_or_v = (deg_of_freedom - 1) % 2
+ # Determine polynomial order:
+ order = (deg_of_freedom + 1) // 2
+ # Return polynomial mesh:
+ return (np.indices(dim_uv)[u_or_v] * self.data_set.a) ** order
+
+ def extract_ramp_params(self, x):
+ '''Extract the ramp parameters of an input vector and return the rest.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Input vector which consists of the vectorised magnetization distribution and the ramp
+ parameters at the end which will be extracted.
+
+ Returns
+ -------
+ result_vector : :class:`~numpy.ndarray` (N=1)
+ Inpput vector without the extracted ramp parameters.
+
+ Notes
+ -----
+ This method should always be used before a vector `x` is processed if it is known that
+ ramp parameters are present so that other functions do not have to bother with them
+ and the :class:`.~ramp` already knows all important parameters for its own functions.
+
+ '''
+ if self.order is not None: # Do nothing if order is None!
+ # Split off ramp parameters and fill cache:
+ x, ramp_params = np.split(x, [-self.n])
+ self.param_cache = ramp_params.reshape((self.deg_of_freedom, self.data_set.count))
+ return x
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 3d2cf4a..80271a1 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -48,22 +48,13 @@ 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)))
- # 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)
+ # Cut ramp parameters if necessary (this also saves the final parameters in the ramp class!):
+ x_opt = costfunction.fwd_model.ramp.extract_ramp_params(x_opt)
# Create and return fitting MagData object:
+ data_set = costfunction.fwd_model.data_set
mag_opt = MagData(data_set.a, np.zeros((3,) + data_set.dim))
mag_opt.set_vector(x_opt, data_set.mask)
- return mag_opt, add_info # TODO: Docstring
+ return mag_opt
def optimize_nonlin(costfunction, first_guess=None):
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index 5d638ae..73cc7de 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -6,10 +6,9 @@
which adds additional constraints to a costfunction to minimize."""
-import abc
-
import numpy as np
from scipy import sparse
+import abc
import jutil.norms as jnorm
import jutil.diff as jdiff
@@ -27,46 +26,50 @@ class Regularisator(object):
Represents a certain constraint for the 3D magnetization distribution whose cost is to minimize
in addition to the derivation from the 2D phase maps. Important is the used `norm` and the
regularisation parameter `lam` (lambda) which determines the weighting between the two cost
- parts (measurements and regularisation).
+ parts (measurements and regularisation). Additional parameters at the end of the input
+ vector, which are not relevant for the regularisation can be discarded by specifying the
+ number in `add_params`.
Attributes
----------
- norm: :class:`~jutil.norm.WeightedNorm`
+ norm : :class:`~jutil.norm.WeightedNorm`
Norm, which is used to determine the cost of the regularisation term.
- lam: float
+ lam : float
Regularisation parameter determining the weighting between measurements and regularisation.
+ add_params : int
+ Number of additional parameters which are not used in the regularisation. Used to cut
+ the input vector into the appropriate size.
- ''' # TODO: ALL dostrings with add_params of input!
+ '''
__metaclass__ = abc.ABCMeta
_log = logging.getLogger(__name__+'.Regularisator')
@abc.abstractmethod
- 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
+ def __init__(self, norm, lam, add_params=0):
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!
+ self.add_params = add_params
+ if self.add_params > 0:
+ self.slice = slice(-add_params)
else:
- self.add_params = None
+ self.slice = slice(None)
self._log.debug('Created '+str(self))
def __call__(self, x):
self._log.debug('Calling __call__')
- # TODO: assert len(x) - self.add_params korrekt
- return self.lam * self.norm(x[:self.add_params])
+ return self.lam * self.norm(x[self.slice])
def __repr__(self):
- self._log.debug('Calling __repr__') # TODO: add_params
- return '%s(norm=%r, lam=%r)' % (self.__class__, self.norm, self.lam)
+ self._log.debug('Calling __repr__')
+ return '%s(norm=%r, lam=%r, add_params=%r)' % (self.__class__, self.norm, self.lam,
+ self.add_params)
def __str__(self):
self._log.debug('Calling __str__')
- return 'Regularisator(norm=%s, lam=%s)' % (self.norm, self.lam)
+ return 'Regularisator(norm=%s, lam=%s, add_params=%s)' % (self.norm, self.lam,
+ self.add_params)
def jac(self, x):
'''Calculate the derivative of the regularisation term for a given magnetic distribution.
@@ -83,7 +86,7 @@ class Regularisator(object):
'''
result = np.zeros_like(x)
- result[:self.add_params] = self.lam * self.norm.jac(x[:self.add_params])
+ result[self.slice] = self.lam * self.norm.jac(x[self.slice])
return result
def hess_dot(self, x, vector):
@@ -105,8 +108,8 @@ class Regularisator(object):
Product of the input `vector` with the Hessian matrix.
'''
- result = np.zeros_like(x)
- result[:self.add_params] = self.lam * self.norm.hess_dot(x, vector[:self.add_params])
+ result = np.zeros_like(vector)
+ result[self.slice] = self.lam * self.norm.hess_dot(x, vector[self.slice])
return result
def hess_diag(self, x):
@@ -128,7 +131,7 @@ class Regularisator(object):
'''
self._log.debug('Calling hess_diag')
result = np.zeros_like(x)
- result[:self.add_params] = self.lam * self.norm.hess_diag(x[:self.add_params])
+ result[self.slice] = self.lam * self.norm.hess_diag(x[self.slice])
return result
diff --git a/pyramid/version.py b/pyramid/version.py
index bf68b51..9d9f9c9 100644
--- a/pyramid/version.py
+++ b/pyramid/version.py
@@ -1,3 +1,3 @@
# THIS FILE IS GENERATED BY THE PYRAMID SETUP.PY
version = "0.1.0-dev"
-hg_revision = "e2e912d70bc2+"
+hg_revision = "711e66d3213b+"
diff --git a/scripts/phasemap/phasemap_from_image.py b/scripts/phasemap/phasemap_from_image.py
index d0d1dd6..35f2175 100644
--- a/scripts/phasemap/phasemap_from_image.py
+++ b/scripts/phasemap/phasemap_from_image.py
@@ -12,25 +12,30 @@ import logging.config
logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
###################################################################################################
-path_mag = 'Arnaud_M.tif'
-path_ele = 'Arnaud_MIP_mask.tif'
-filename = 'phasemap_tif_martial_skyrmion.nc'
-a = 2
+path_mag = 'trevor_magnetite_m20.bmp'
+path_mask = 'trevor_magnetite_mask20.bmp'
+filename = 'phasemap_bmp_trevor_magnetite_m20.nc'
+a = 0.4 # nm
dim_uv = None
max_phase = 1
threshold = 0.5
offset = 0
+flip_up_down = True
###################################################################################################
# Load images:
im_mag = Image.open(os.path.join(py.DIR_FILES, 'images', path_mag)).convert('P')
-im_ele = Image.open(os.path.join(py.DIR_FILES, 'images', path_ele)).convert('P')
+if flip_up_down:
+ im_mag = im_mag.transpose(Image.FLIP_TOP_BOTTOM)
+im_mask = Image.open(os.path.join(py.DIR_FILES, 'images', path_mask)).convert('P')
+if flip_up_down:
+ im_mask = im_mask.transpose(Image.FLIP_TOP_BOTTOM)
if dim_uv is not None:
im_mag = im_mag.resize(dim_uv)
- im_ele = im_ele.resize(dim_uv)
+ im_mask = im_mask.resize(dim_uv)
# Calculate phase and mask:
phase = np.asarray(im_mag)/255.*max_phase - offset
-mask = np.where(np.asarray(im_ele)/255. >= threshold, True, False)
+mask = np.where(np.asarray(im_mask)/255. >= threshold, True, False)
# Create and save PhaseMap object:
phase_map = py.PhaseMap(a, phase, mask, confidence=None, unit='rad')
diff --git a/scripts/reconstruction/reconstruction_2d_from_phasemap.py b/scripts/reconstruction/reconstruction_2d_from_phasemap.py
index 4e17081..1efb602 100644
--- a/scripts/reconstruction/reconstruction_2d_from_phasemap.py
+++ b/scripts/reconstruction/reconstruction_2d_from_phasemap.py
@@ -11,13 +11,12 @@ import logging.config
logging.config.fileConfig(pr.LOGGING_CONFIG, disable_existing_loggers=False)
###################################################################################################
-phase_name = 'phasemap_dm3_zi_an_elongated_nanorod'
+phase_name = 'phasemap_bmp_trevor_magnetite_m20'
b_0 = 0.6 # in T
lam = 1E-3
max_iter = 100
buffer_pixel = 0
-fit_ramps = True
-fit_offsets = True
+order = 1
###################################################################################################
# Load phase map:
@@ -29,22 +28,18 @@ dim = (1,) + phase_map.dim_uv
data = pr.DataSet(phase_map.a, dim, b_0)
data.append(phase_map, pr.SimpleProjector(dim))
data.set_3d_mask()
-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)
+
+fwd_model = pr.ForwardModel(data, order)
+reg = pr.FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
cost = pr.Costfunction(fwd_model, reg)
# Reconstruct and save:
with TakeTime('reconstruction time'):
- 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:
+ mag_data_rec = pr.reconstruction.optimize_linear(cost, max_iter=max_iter)
+ add_info = cost.fwd_model.ramp.param_cache
+if order >= 1:
+ offset, ramp = add_info[0][0], (add_info[1][0], add_info[2][0])
+elif order == 0:
offset = add_info[0][0]
mag_data_buffer = mag_data_rec.copy()
mag_data_rec.crop((0, buffer_pixel, buffer_pixel))
@@ -55,15 +50,18 @@ mag_data_rec.save_to_netcdf4(mag_name+'.nc')
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 -= fwd_model.ramp(index=0)
+phase_map.display_combined('Input Phase (ramp corrected)')
phase_map_rec = pr.pm(mag_data_rec)
-phase_map_rec.mask = phase_map.mask
title = 'Reconstructed Phase'
-if fit_offsets:
+if order >= 0:
print 'offset:', offset
title += ', fitted Offset: {:.2g} [rad]'.format(offset)
-if fit_ramps:
+if order >= 1:
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: {:.2g})'.format(difference))
+if order is not None:
+ fwd_model.ramp(0).display_combined('Fitted Ramp')
diff --git a/scripts/reconstruction/reconstruction_3d_from_magdata.py b/scripts/reconstruction/reconstruction_3d_from_magdata.py
index 8d1a763..2fa182e 100644
--- a/scripts/reconstruction/reconstruction_3d_from_magdata.py
+++ b/scripts/reconstruction/reconstruction_3d_from_magdata.py
@@ -22,8 +22,8 @@ max_iter = 100
use_internal_mask = True
offset_max = 2
ramp_max = 0.01
-fit_ramps = True
-fit_offsets = True
+order = 1
+show_input = True
###################################################################################################
# Load magnetization distribution:
@@ -59,8 +59,9 @@ if noise != 0:
data.phase_maps[i] = phase_map
# Plot input:
-for i, phase_map in enumerate(data.phase_maps):
- phase_map.display_phase()
+if show_input:
+ for i, phase_map in enumerate(data.phase_maps):
+ phase_map.display_phase()
# Construct mask:
if use_internal_mask:
@@ -69,30 +70,19 @@ else:
data.set_3d_mask() # Construct mask from 2D phase masks!
# 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)
+fwd_model = pr.ForwardModel(data, ramp_order=order)
+reg = pr.FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
cost = pr.Costfunction(fwd_model, reg)
# Reconstruct and save:
with TakeTime('reconstruction time'):
- 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_data_rec = pr.reconstruction.optimize_linear(cost, max_iter=max_iter)
+ param_cache = cost.fwd_model.ramp.param_cache
mag_name_rec = mag_name.replace('magdata', 'magdata_rec')
mag_data_rec.save_to_netcdf4(mag_name_rec+'.nc')
# Plot stuff:
+
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.))
diff --git a/scripts/temp/test_dim_uv_projection.py b/scripts/temp/test_dim_uv_projection.py
index b231245..2634862 100644
--- a/scripts/temp/test_dim_uv_projection.py
+++ b/scripts/temp/test_dim_uv_projection.py
@@ -9,11 +9,10 @@ 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)
+dim = (1, 5, 5)
+dim_uv = None
+width = (1, 2, 2)
+center = (0.5, 2.5, 2.5)
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()
@@ -22,3 +21,4 @@ mag_proj = projector(mag_data)
mag_proj.quiver_plot()
pr.pm(mag_proj).display_combined()
weight = np.array(projector.weight.todense())
+pr.an.phase_mag_vortex(dim, 10, center, 12.5, 1).display_combined()
diff --git a/scripts/temp/warning_control.py b/scripts/temp/warning_control.py
new file mode 100644
index 0000000..324fadb
--- /dev/null
+++ b/scripts/temp/warning_control.py
@@ -0,0 +1,24 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Jul 19 13:03:30 2015
+
+@author: Jan
+"""
+
+
+import pyramid as pr
+import warnings
+import logging.config
+
+
+logging.config.fileConfig(pr.LOGGING_CONFIG, disable_existing_loggers=False)
+
+#with warnings.catch_warnings():
+# warnings.simplefilter("ignore")
+logging.disable = True
+warnings.warn('TEST')
+from mayavi import mlab
+warnings.warn('TEST')
+logging.disable = False
+
+print mlab
diff --git a/tests/test_analytic.py b/tests/test_analytic.py
index dbdd9d6..7b7637f 100644
--- a/tests/test_analytic.py
+++ b/tests/test_analytic.py
@@ -19,7 +19,7 @@ class TestCaseAnalytic(unittest.TestCase):
dim = (4, 4, 4)
a = 10.0
phi = pi/4
- center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2-0.5)
+ center = (dim[0]/2, dim[1]/2, dim[2]/2)
radius = dim[2]/4
def test_phase_mag_slab(self):
@@ -42,7 +42,8 @@ class TestCaseAnalytic(unittest.TestCase):
radius = self.dim[2]/4
phase = an.phase_mag_sphere(self.dim, self.a, self.phi, self.center, radius).phase
reference = np.load(os.path.join(self.path, 'ref_phase_sphere.npy'))
- assert_allclose(phase, reference, err_msg='Unexpected behavior in phase_mag_sphere()')
+ assert_allclose(phase, reference, err_msg='Unexpected behavior in phase_mag_sphere()',
+ atol=1E-10)
def test_phase_mag_vortex(self):
'''Test of the phase_mag_vortex method.'''
diff --git a/tests/test_magcreator.py b/tests/test_magcreator.py
index 6ba7371..980d631 100644
--- a/tests/test_magcreator.py
+++ b/tests/test_magcreator.py
@@ -17,14 +17,14 @@ class TestCaseMagCreator(unittest.TestCase):
path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_magcreator')
def test_shape_slab(self):
- test_slab = mc.Shapes.slab((5, 6, 7), (2, 3, 4), (1, 3, 5))
+ test_slab = mc.Shapes.slab((5, 6, 7), (2.5, 3.5, 4.5), (1, 3, 5))
assert_allclose(test_slab, np.load(os.path.join(self.path, 'ref_slab.npy')),
err_msg='Created slab does not match expectation!')
def test_shape_disc(self):
- test_disc_z = mc.Shapes.disc((5, 6, 7), (2, 3, 4), 2, 3, 'z')
- test_disc_y = mc.Shapes.disc((5, 6, 7), (2, 3, 4), 2, 3, 'y')
- test_disc_x = mc.Shapes.disc((5, 6, 7), (2, 3, 4), 2, 3, 'x')
+ test_disc_z = mc.Shapes.disc((5, 6, 7), (2.5, 3.5, 4.5), 2, 3, 'z')
+ test_disc_y = mc.Shapes.disc((5, 6, 7), (2.5, 3.5, 4.5), 2, 3, 'y')
+ test_disc_x = mc.Shapes.disc((5, 6, 7), (2.5, 3.5, 4.5), 2, 3, 'x')
assert_allclose(test_disc_z, np.load(os.path.join(self.path, 'ref_disc_z.npy')),
err_msg='Created disc in z-direction does not match expectation!')
assert_allclose(test_disc_y, np.load(os.path.join(self.path, 'ref_disc_y.npy')),
@@ -33,23 +33,23 @@ class TestCaseMagCreator(unittest.TestCase):
err_msg='Created disc in x-direction does not match expectation!')
def test_shape_ellipse(self):
- test_ellipse_z = mc.Shapes.ellipse((7, 8, 9), (3, 4, 5), (3, 5), 1, 'z')
- test_ellipse_y = mc.Shapes.ellipse((7, 8, 9), (3, 4, 5), (3, 5), 1, 'y')
- test_ellipse_x = mc.Shapes.ellipse((7, 8, 9), (3, 4, 5), (3, 5), 1, 'x')
+ test_ellipse_z = mc.Shapes.ellipse((7, 8, 9), (3.5, 4.5, 5.5), (3, 5), 1, 'z')
+ test_ellipse_y = mc.Shapes.ellipse((7, 8, 9), (3.5, 4.5, 5.5), (3, 5), 1, 'y')
+ test_ellipse_x = mc.Shapes.ellipse((7, 8, 9), (3.5, 4.5, 5.5), (3, 5), 1, 'x')
assert_allclose(test_ellipse_z, np.load(os.path.join(self.path, 'ref_ellipse_z.npy')),
- err_msg='Created ellipse does not match expectation!')
+ err_msg='Created ellipse does not match expectation (z)!')
assert_allclose(test_ellipse_y, np.load(os.path.join(self.path, 'ref_ellipse_y.npy')),
- err_msg='Created ellipse does not match expectation!')
+ err_msg='Created ellipse does not match expectation (y)!')
assert_allclose(test_ellipse_x, np.load(os.path.join(self.path, 'ref_ellipse_x.npy')),
- err_msg='Created ellipse does not match expectation!')
+ err_msg='Created ellipse does not match expectation (x)!')
def test_shape_sphere(self):
- test_sphere = mc.Shapes.sphere((5, 6, 7), (2, 3, 4), 2)
+ test_sphere = mc.Shapes.sphere((5, 6, 7), (2.5, 3.5, 4.5), 2)
assert_allclose(test_sphere, np.load(os.path.join(self.path, 'ref_sphere.npy')),
err_msg='Created sphere does not match expectation!')
def test_shape_ellipsoid(self):
- test_ellipsoid = mc.Shapes.ellipsoid((7, 8, 9), (3, 4, 4), (3, 5, 7))
+ test_ellipsoid = mc.Shapes.ellipsoid((7, 8, 9), (3.5, 4.5, 4.5), (3, 5, 7))
assert_allclose(test_ellipsoid, np.load(os.path.join(self.path, 'ref_ellipsoid.npy')),
err_msg='Created ellipsoid does not match expectation!')
@@ -70,13 +70,13 @@ class TestCaseMagCreator(unittest.TestCase):
err_msg='Created pixel does not match expectation!')
def test_create_mag_dist_homog(self):
- mag_shape = mc.Shapes.disc((1, 10, 10), (0, 4.5, 4.5), 3, 1)
+ mag_shape = mc.Shapes.disc((1, 10, 10), (0, 5, 5), 3, 1)
magnitude = mc.create_mag_dist_homog(mag_shape, pi/4)
assert_allclose(magnitude, np.load(os.path.join(self.path, 'ref_mag_disc.npy')),
err_msg='Created homog. magnetic distribution does not match expectation')
def test_create_mag_dist_vortex(self):
- mag_shape = mc.Shapes.disc((1, 10, 10), (0, 4.5, 4.5), 3, 1)
+ mag_shape = mc.Shapes.disc((1, 10, 10), (0, 5, 5), 3, 1)
magnitude = mc.create_mag_dist_vortex(mag_shape)
assert_allclose(magnitude, np.load(os.path.join(self.path, 'ref_mag_vort.npy')),
err_msg='Created vortex magnetic distribution does not match expectation')
diff --git a/tests/test_magcreator/ref_ellipsoid.npy b/tests/test_magcreator/ref_ellipsoid.npy
index 8abca7c48658e9f9df018d7869ba2dc4461b1d9a..f040dfd61b0d6ab846cce0e0a9f9defb108a4eb4 100644
GIT binary patch
delta 39
acmX@Xa)M=pGov&EBLpyDGA73{z5)Pd@&uUx
delta 39
WcmX@Xa)M=pGov&kZZJ8H@f84PUId>2
diff --git a/tests/test_regularisator.py b/tests/test_regularisator.py
index 39ef31c..656ff94 100644
--- a/tests/test_regularisator.py
+++ b/tests/test_regularisator.py
@@ -128,7 +128,7 @@ class TestCaseFirstOrderRegularisator(unittest.TestCase):
err_msg='Unexpected behaviour in hess_dot()!')
def test_hess_diag(self):
- hess_diag = self.reg.hess_diag(np.ones(self.n))
+ hess_diag = self.reg.hess_diag(np.ones(3*self.n)) # derivatives in all directions!
hess_diag_ref = np.zeros(3*self.n) # derivatives in all directions!
first_order_jac_ref = np.load(os.path.join(self.path, 'first_order_jac_ref.npy'))
hess_diag_ref[0:self.n] = np.diag(first_order_jac_ref)
--
GitLab