utils: New subpackage with the GUIs, pm and shortcut reconstructions for 2D/3D.

pyramid/__init__.py

@@ -53,7 +53,7 @@ from . import magcreator
 from . import reconstruction
 from . import fft
 from . import fieldconverter
-from . import gui
+from . import utils
 from .costfunction import *
 from .dataset import *
 from .diagnostics import *
@@ -76,7 +76,7 @@ _log = logging.getLogger(__name__)
 _log.info("Starting Pyramid V{} HG{}".format(__version__, __hg_revision__))
 del logging
 
-__all__ = ['analytic', 'magcreator', 'reconstruction', 'fft', 'fieldconverter', 'gui']
+__all__ = ['analytic', 'magcreator', 'reconstruction', 'fft', 'fieldconverter', 'utils']
 __all__.extend(costfunction.__all__)
 __all__.extend(dataset.__all__)
 __all__.extend(diagnostics.__all__)

pyramid/fielddata.py

@@ -472,8 +472,7 @@ class VectorData(FieldData):
         for i, values in enumerate(pad_values):
             assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
             pv[2 * i:2 * (i + 1)] = values
-        self.field = np.pad(self.field,
-                            ((0, 0), (pv[0], pv[1]), (pv[2], pv[3]), (pv[4], pv[5])),
+        self.field = np.pad(self.field, ((0, 0), (pv[0], pv[1]), (pv[2], pv[3]), (pv[4], pv[5])),
                             mode='constant')
 
     def crop(self, crop_values):

pyramid/gui/__init__.py

-# -*- coding: utf-8 -*-
-"""Package containing Pyramid GUI functions."""
-
-from .phasemap_creator import phasemap_creator
-from .mag_slicer import mag_slicer
-
-__all__ = ['phasemap_creator', 'mag_slicer']

pyramid/magcreator/__init__.py

 # -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
 
-from .magcreator import *
 from . import shapes
 from . import examples
+from .magcreator import *
 
 import logging
 _log = logging.getLogger(__name__)
 del logging
 
+
 __all__ = ['shapes', 'examples']
 __all__.extend(magcreator.__all__)

pyramid/magcreator/examples.py

 # -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
 
 import logging
 
@@ -10,6 +13,12 @@ from . import magcreator as mc
 from . import shapes
 from ..fielddata import VectorData
 
+
+__all__ = ['pyramid_logo', 'singularity', 'homog_pixel', 'homog_slab', 'homog_disc',
+           'homog_sphere', 'homog_filament', 'homog_alternating_filament',
+           'homog_array_sphere_disc_slab', 'homog_random_pixels', 'homog_random_slabs',
+           'vortex_slab', 'vortex_disc', 'vortex_alternating_discs', 'vortex_sphere',
+           'vortex_horseshoe', 'core_shell_disc', 'core_shell_sphere']
 _log = logging.getLogger(__name__)
 
 

pyramid/phasemapper.py

@@ -13,13 +13,11 @@ import logging
 import numpy as np
 from numpy import pi
 
-from pyramid import fft
-from pyramid.kernel import Kernel
-from pyramid.fielddata import VectorData, ScalarData
-from pyramid.phasemap import PhaseMap
-from pyramid.projector import RotTiltProjector, XTiltProjector, YTiltProjector, SimpleProjector
+from . import fft
+from .fielddata import VectorData, ScalarData
+from .phasemap import PhaseMap
 
-__all__ = ['PhaseMapperRDFC', 'PhaseMapperFDFC', 'PhaseMapperMIP', 'pm']
+__all__ = ['PhaseMapperRDFC', 'PhaseMapperFDFC', 'PhaseMapperMIP']
 _log = logging.getLogger(__name__)
 
 PHI_0 = 2067.83  # magnetic flux in T*nm²
@@ -409,48 +407,3 @@ class PhaseMapperMIP(PhaseMapper):
 
         """
         raise NotImplementedError()  # TODO: Implement right!
-
-
-def pm(mag_data, mode='z', b_0=1, **kwargs):
-    """Convenience function for fast magnetic phase mapping.
-
-    Parameters
-    ----------
-    mag_data : :class:`~.VectorData`
-        A :class:`~.VectorData` object, from which the projected phase map should be calculated.
-    mode: {'z', 'y', 'x', 'x-tilt', 'y-tilt', 'rot-tilt'}, optional
-        Projection mode which determines the :class:`~.pyramid.projector.Projector` subclass, which
-        is used for the projection. Default is a simple projection along the `z`-direction.
-    b_0 : float, optional
-        Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
-    **kwargs : additional arguments
-        Additional arguments like `dim_uv`, 'tilt' or 'rotation', which are passed to the
-        projector-constructor, defined by the `mode`.
-
-    Returns
-    -------
-    phase_map : :class:`~pyramid.phasemap.PhaseMap`
-        The calculated phase map as a :class:`~.PhaseMap` object.
-
-    """
-    _log.debug('Calling pm')
-    # Determine projection mode:
-    if mode == 'rot-tilt':
-        projector = RotTiltProjector(mag_data.dim, **kwargs)
-    elif mode == 'x-tilt':
-        projector = XTiltProjector(mag_data.dim, **kwargs)
-    elif mode == 'y-tilt':
-        projector = YTiltProjector(mag_data.dim, **kwargs)
-    elif mode in ['x', 'y', 'z']:
-        projector = SimpleProjector(mag_data.dim, axis=mode, **kwargs)
-    else:
-        raise ValueError("Invalid mode (use 'x', 'y', 'z', 'x-tilt', 'y-tilt' or 'rot-tilt')")
-    # Project:
-    mag_proj = projector(mag_data)
-    # Set up phasemapper and map phase:
-    phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv, b_0=b_0))
-    phase_map = phasemapper(mag_proj)
-    # Get mask from magdata:
-    phase_map.mask = mag_proj.get_mask()[0, ...]
-    # Return phase:
-    return phase_map

pyramid/tests/test_phasemapper.py

@@ -10,7 +10,7 @@ from numpy.testing import assert_allclose
 from pyramid.kernel import Kernel
 from pyramid.fielddata import VectorData, ScalarData
 from pyramid.phasemap import PhaseMap
-from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC, PhaseMapperMIP, pm
+from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC, PhaseMapperMIP
 
 
 class TestCasePhaseMapperRDFC(unittest.TestCase):
@@ -178,24 +178,6 @@ class TestCasePhaseMapperMIP(unittest.TestCase):
         self.assertRaises(NotImplementedError, self.mapper.jac_T_dot, None)
 
 
-class TestCasePM(unittest.TestCase):
-    def setUp(self):
-        self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_phasemapper')
-        self.mag_proj = VectorData.load_from_hdf5(os.path.join(self.path, 'mag_proj.hdf5'))
-
-    def tearDown(self):
-        self.path = None
-        self.mag_proj = None
-        self.mapper = None
-
-    def test_pm(self):
-        phase_ref = PhaseMap.load_from_hdf5(os.path.join(self.path, 'phase_map.hdf5'))
-        phase_map = pm(self.mag_proj)
-        assert_allclose(phase_map.phase, phase_ref.phase, atol=1E-7,
-                        err_msg='Unexpected behavior in pm()!')
-        assert_allclose(phase_map.a, phase_ref.a, err_msg='Unexpected behavior in pm()!')
-
-
 if __name__ == '__main__':
     import nose
     nose.run(defaultTest=__name__)

pyramid/tests/test_pm.py

+# -*- coding: utf-8 -*-
+"""Testcase for the pm function."""
+
+import os
+import unittest
+
+from numpy.testing import assert_allclose
+
+from pyramid.fielddata import VectorData
+from pyramid.phasemap import PhaseMap
+from pyramid.utils import pm
+
+
+class TestCasePM(unittest.TestCase):
+    def setUp(self):
+        self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_phasemapper')
+        self.mag_proj = VectorData.load_from_hdf5(os.path.join(self.path, 'mag_proj.hdf5'))
+
+    def tearDown(self):
+        self.path = None
+        self.mag_proj = None
+        self.mapper = None
+
+    def test_pm(self):
+        phase_ref = PhaseMap.load_from_hdf5(os.path.join(self.path, 'phase_map.hdf5'))
+        phase_map = pm(self.mag_proj)
+        assert_allclose(phase_map.phase, phase_ref.phase, atol=1E-7,
+                        err_msg='Unexpected behavior in pm()!')
+        assert_allclose(phase_map.a, phase_ref.a, err_msg='Unexpected behavior in pm()!')
+
+
+if __name__ == '__main__':
+    import nose
+    nose.run(defaultTest=__name__)

pyramid/utils/__init__.py

+# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""Package containing Pyramid utility functions."""
+
+from .pm import pm
+from .reconstruction_2d_from_phasemap import reconstruction_2d_from_phasemap
+from .reconstruction_3d_from_magdata import reconstruction_3d_from_magdata
+from .phasemap_creator import gui_phasemap_creator
+from .mag_slicer import gui_mag_slicer
+
+import logging
+_log = logging.getLogger(__name__)
+del logging
+
+
+__all__ = ['pm', 'reconstruction_2d_from_phasemap', 'reconstruction_3d_from_magdata',
+           'gui_phasemap_creator', 'gui_mag_slicer']

pyramid/gui/mag_slicer.py → pyramid/utils/mag_slicer.py

 # -*- coding: utf-8 -*-
-# Form implementation generated from reading ui file 'mag_slicer2.ui'
+# Form implementation generated from reading ui file 'mag_slicer.ui'
 #
 # Created: Sun Aug 31 20:39:52 2014
 #      by: PyQt4 UI code generator 4.9.6
@@ -17,7 +17,10 @@ from matplotlib.figure import Figure
 from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
 from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar
 
-import pyramid as pr
+from ..projector import SimpleProjector
+from ..kernel import Kernel
+from ..phasemapper import PhaseMapperRDFC
+from ..fielddata import VectorData
 
 
 ui_location = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'mag_slicer.ui')
@@ -87,13 +90,13 @@ class Main(QMainWindow, UI_MainWindow):
                 gain = 'auto'
             else:
                 gain = self.spinBoxGain.value()
-            self.projector = pr.SimpleProjector(self.mag_data.dim, axis=self.mode)
+            self.projector = SimpleProjector(self.mag_data.dim, axis=self.mode)
             self.spinBoxSlice.setMaximum(length)
             self.scrollBarSlice.setMaximum(length)
             self.spinBoxSlice.setValue(int(length / 2.))
             self.update_slice()
-            kernel = pr.Kernel(self.mag_data.a, self.projector.dim_uv)
-            self.phase_mapper = pr.PhaseMapperRDFC(kernel)
+            kernel = Kernel(self.mag_data.a, self.projector.dim_uv)
+            self.phase_mapper = PhaseMapperRDFC(kernel)
             self.phase_map = self.phase_mapper(self.projector(self.mag_data))
             self.canvasPhase.figure.axes[0].clear()
             self.phase_map.display_phase(axis=self.canvasPhase.figure.axes[0], cbar=False)
@@ -124,7 +127,7 @@ class Main(QMainWindow, UI_MainWindow):
             return  # Abort if no conf_path is selected!
         import hyperspy.api as hs
         print(hs.load(mag_file))
-        self.mag_data = pr.VectorData.load_from_hdf5(mag_file)
+        self.mag_data = VectorData.load_from_hdf5(mag_file)
         if not self.is_mag_data_loaded:
             self.addmpl()
         self.is_mag_data_loaded = True
@@ -132,7 +135,7 @@ class Main(QMainWindow, UI_MainWindow):
         self.update_phase()
 
 
-def mag_slicer():
+def gui_mag_slicer():
     app = QtGui.QApplication(sys.argv)
     main = Main()
     main.show()

pyramid/gui/phasemap_creator.py → pyramid/utils/phasemap_creator.py

@@ -167,7 +167,7 @@ class Main(QMainWindow, UI_MainWindow):
             return  # Abort if no export_path is selected or self.phase_map doesn't exist yet!
 
 
-def phasemap_creator():
+def gui_phasemap_creator():
     app = QtGui.QApplication(sys.argv)
     main = Main()
     main.show()

pyramid/utils/pm.py

+# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+import logging
+
+from ..kernel import Kernel
+from ..phasemapper import PhaseMapperRDFC
+from ..projector import RotTiltProjector, XTiltProjector, YTiltProjector, SimpleProjector
+
+
+__all__ = ['pm']
+_log = logging.getLogger(__name__)
+
+
+def pm(mag_data, mode='z', b_0=1, **kwargs):
+    """Convenience function for fast magnetic phase mapping.
+
+    Parameters
+    ----------
+    mag_data : :class:`~.VectorData`
+        A :class:`~.VectorData` object, from which the projected phase map should be calculated.
+    mode: {'z', 'y', 'x', 'x-tilt', 'y-tilt', 'rot-tilt'}, optional
+        Projection mode which determines the :class:`~.pyramid.projector.Projector` subclass, which
+        is used for the projection. Default is a simple projection along the `z`-direction.
+    b_0 : float, optional
+        Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
+    **kwargs : additional arguments
+        Additional arguments like `dim_uv`, 'tilt' or 'rotation', which are passed to the
+        projector-constructor, defined by the `mode`.
+
+    Returns
+    -------
+    phase_map : :class:`~pyramid.phasemap.PhaseMap`
+        The calculated phase map as a :class:`~.PhaseMap` object.
+
+    """
+    _log.debug('Calling pm')
+    # Determine projection mode:
+    if mode == 'rot-tilt':
+        projector = RotTiltProjector(mag_data.dim, **kwargs)
+    elif mode == 'x-tilt':
+        projector = XTiltProjector(mag_data.dim, **kwargs)
+    elif mode == 'y-tilt':
+        projector = YTiltProjector(mag_data.dim, **kwargs)
+    elif mode in ['x', 'y', 'z']:
+        projector = SimpleProjector(mag_data.dim, axis=mode, **kwargs)
+    else:
+        raise ValueError("Invalid mode (use 'x', 'y', 'z', 'x-tilt', 'y-tilt' or 'rot-tilt')")
+    # Project:
+    mag_proj = projector(mag_data)
+    # Set up phasemapper and map phase:
+    phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv, b_0=b_0))
+    phase_map = phasemapper(mag_proj)
+    # Get mask from magdata:
+    phase_map.mask = mag_proj.get_mask()[0, ...]
+    # Return phase:
+    return phase_map

pyramid/utils/reconstruction_2d_from_phasemap.py

+# -*- coding: utf-8 -*-
+"""Reconstruct a magnetization distributions from a single phase map."""
+
+import numpy as np
+
+from jutil.taketime import TakeTime
+
+from .. import reconstruction
+from ..phasemap import PhaseMap
+from ..dataset import DataSet
+from ..projector import SimpleProjector
+from ..regularisator import FirstOrderRegularisator
+from ..forwardmodel import ForwardModel
+from ..costfunction import Costfunction
+from .pm import pm
+
+
+def reconstruction_2d_from_phasemap(filename, b_0=1, lam=1E-3, max_iter=100, ramp_order=1,
+                                    plot_results=False, ar_dens=None):
+
+    # Load phase map:
+    phase_map = PhaseMap.load_from_hdf5(filename)
+    # Construct regularisator, forward model and costfunction:
+    dim = (1,) + phase_map.dim_uv
+    data = DataSet(phase_map.a, dim, b_0)
+    data.append(phase_map, SimpleProjector(dim))
+    data.set_3d_mask()
+    fwd_model = ForwardModel(data, ramp_order)
+    reg = FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
+    cost = Costfunction(fwd_model, reg)
+    # Reconstruct:
+    with TakeTime('reconstruction time'):
+        mag_data_rec = reconstruction.optimize_linear(cost, max_iter=max_iter)
+        param_cache = cost.fwd_model.ramp.param_cache
+    if ramp_order >= 1:
+        offset, ramp = param_cache[0][0], (param_cache[1][0], param_cache[2][0])
+    elif ramp_order == 0:
+        offset, ramp = param_cache[0][0], (0, 0)
+    else:
+        offset, ramp = 0, (0, 0)
+    # Plot stuff:
+    if plot_results:
+        if ar_dens is None:
+            ar_dens = np.max(dim) // 128
+        mag_data_rec.quiver_plot('Reconstructed Distribution', ar_dens=ar_dens)
+        phase_map.display_combined('Input Phase')
+        phase_map -= fwd_model.ramp(index=0)
+        phase_map.display_combined('Input Phase (ramp corrected)')
+        phase_map_rec = pm(mag_data_rec)
+        title = 'Reconstructed Phase'
+        if ramp_order >= 0:
+            print('offset:', offset)
+            title += ', fitted Offset: {:.2g} [rad]'.format(offset)
+        if ramp_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 ramp_order is not None:
+            fwd_model.ramp(0).display_combined('Fitted Ramp')
+    # Return reconstructed magnetisation distribution and cost function:
+    return mag_data_rec, cost

scripts/reconstruction/reconstruction_3d_from_magdata.py → pyramid/utils/reconstruction_3d_from_magdata.py

 # -*- coding: utf-8 -*-
 """Reconstruct a magnetization distributions from phase maps created from it."""
 
-import multiprocessing as mp
-
-import matplotlib.pyplot as plt
-
 import numpy as np
 
-import pyramid as pr
+import multiprocessing as mp
+
 from jutil.taketime import TakeTime
 
+from .. import reconstruction
+from ..fielddata import VectorData
+from ..dataset import DataSet
+from ..projector import XTiltProjector, YTiltProjector
+from ..ramp import Ramp
+from ..regularisator import FirstOrderRegularisator
+from ..forwardmodel import ForwardModel, DistributedForwardModel
+from ..costfunction import Costfunction
 
-###################################################################################################
-mag_name = 'magdata_mc_alternating_vortices'
-dim_uv = None
-angles = np.linspace(-90, 90, num=19)
-axes = {'x': True, 'y': True}
-b_0 = 1
-noise = 0
-lam = 1E-5
-max_iter = 100
-use_internal_mask = True
-offset_max = 0
-ramp_max = 0
-order = 1
-show_input = True
-ar_dens = 1 # 'auto'
-nprocs = mp.cpu_count()  # or 1 for non-multiprocessing
-###################################################################################################
 
+def reconstruction_3d_from_magdata(filename, b_0=1, lam=1E-3, max_iter=100, ramp_order=1,
+                                   angles=np.linspace(-90, 90, num=19), dim_uv=None,
+                                   axes=(True, True), noise=0, offset_max=0, ramp_max=0,
+                                   use_internal_mask=True, plot_results=False, plot_input=False,
+                                   ar_dens=None, multicore=True):
 
-if __name__ == '__main__':
-    mp.freeze_support()
-    # Load magnetization distribution:
-    mag_data = pr.VectorData.load_from_hdf5(mag_name + '.hdf5')
+    mag_data = VectorData.load_from_hdf5(filename)
     dim = mag_data.dim
-    if ar_dens == 'auto':
-        ar_dens = np.ceil(np.max(dim) / 64)
-    # Construct data set and regularisator:
-    data = pr.DataSet(mag_data.a, mag_data.dim, b_0)
+    # Load magnetization distribution:
+    if ar_dens is None:
+        ar_dens = np.max(dim) // 64
+    data = DataSet(mag_data.a, mag_data.dim, b_0)
     # Construct projectors:
     projectors = []
+    # Construct data set and regularisator:
     for angle in angles:
         angle_rad = angle * np.pi / 180
-        if axes['x']:
-            projectors.append(pr.XTiltProjector(mag_data.dim, angle_rad, dim_uv))
-        if axes['y']:
-            projectors.append(pr.YTiltProjector(mag_data.dim, angle_rad, dim_uv))
-    # Add projectors and construct according phase maps:
+        if axes[0]:
+            projectors.append(XTiltProjector(mag_data.dim, angle_rad, dim_uv))
+        if axes[1]:
+            projectors.append(YTiltProjector(mag_data.dim, angle_rad, dim_uv))
     data.projectors = projectors
     data.phase_maps = data.create_phase_maps(mag_data)
+    # Add projectors and construct according phase maps:
     for i, phase_map in enumerate(data.phase_maps):
         offset = np.random.uniform(-offset_max, offset_max)
         ramp_u = np.random.uniform(-ramp_max, ramp_max)
         ramp_v = np.random.uniform(-ramp_max, ramp_max)
-        phase_map += pr.Ramp.create_ramp(phase_map.a, phase_map.dim_uv, (offset, ramp_u, ramp_v))
+        phase_map += Ramp.create_ramp(phase_map.a, phase_map.dim_uv, (offset, ramp_u, ramp_v))
     # Add noise if necessary:
     if noise != 0:
         for i, phase_map in enumerate(data.phase_maps):
@@ -66,25 +57,27 @@ if __name__ == '__main__':
     else:
         data.set_3d_mask()  # Construct mask from 2D phase masks!
     # Construct regularisator, forward model and costfunction:
-    if nprocs > 1:
-        fwd_model = pr.DistributedForwardModel(data, ramp_order=order, nprocs=nprocs)
+    if multicore:
+        mp.freeze_support()
+        fwd_model = DistributedForwardModel(data, ramp_order=ramp_order, nprocs=mp.cpu_count())
     else:
-        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)
+        fwd_model = ForwardModel(data, ramp_order=ramp_order)
+    reg = FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
+    cost = Costfunction(fwd_model, reg)
     # Reconstruct and save:
     with TakeTime('reconstruction time'):
-        mag_data_rec = pr.reconstruction.optimize_linear(cost, max_iter=max_iter)
-        param_cache = cost.fwd_model.ramp.param_cache
+        mag_data_rec = reconstruction.optimize_linear(cost, max_iter=max_iter)
+    # Finalize ForwardModel (returns workers if multicore):
     fwd_model.finalize()
-    mag_name_rec = mag_name.replace('magdata', 'magdata_rec')
-    mag_data_rec.save_to_hdf5(mag_name_rec + '.hdf5', overwrite=True)
-    # Plot stuff:
-    data.display_mask(ar_dens=ar_dens)
-    mag_data.quiver_plot3d('Original Distribution', ar_dens=ar_dens)
-    mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=ar_dens)
-    mag_data_rec.quiver_plot3d('Reconstructed Distribution', ar_dens=ar_dens, coloring='amplitude')
     # Plot input:
-    if show_input:
+    if plot_input:
         data.display_phase()
-    plt.show()
+    # Plot results:
+    if plot_results:
+        data.display_mask(ar_dens=ar_dens)
+        mag_data.quiver_plot3d('Original Distribution', ar_dens=ar_dens)
+        mag_data_rec.quiver_plot3d('Reconstructed Distribution (angle)', ar_dens=ar_dens)
+        mag_data_rec.quiver_plot3d('Reconstructed Distribution (amplitude)',
+                                   ar_dens=ar_dens, coloring='amplitude')
+    # Return reconstructed magnetisation distribution and cost function:
+    return mag_data_rec, cost

scripts/reconstruction/reconstruction_2d_from_phasemap.py

-# -*- coding: utf-8 -*-
-"""Reconstruct a magnetization distributions from a single phase map."""
-
-import matplotlib.pyplot as plt
-import numpy as np
-
-import pyramid as pr
-from jutil.taketime import TakeTime
-
-
-###################################################################################################
-phase_name = 'phasemap_dm3_14p5kx_m0150mT_q3_pha_sb400_sc512_magn'
-b_0 = 1  # in T
-lam = 1E-1
-max_iter = 100
-buffer_pixel = 0
-order = 1
-###################################################################################################
-
-
-# Load phase map:
-phase_map = pr.PhaseMap.load_from_hdf5(phase_name + '.hdf5')
-phase_map.pad((buffer_pixel, buffer_pixel))
-dim = (1,) + phase_map.dim_uv
-
-ar_dens = np.max(dim) // 128
-
-# Construct regularisator, forward model and costfunction:
-data = pr.DataSet(phase_map.a, dim, b_0)
-data.append(phase_map, pr.SimpleProjector(dim))
-data.append(phase_map, pr.SimpleProjector(dim))
-data.set_3d_mask()
-
-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 = pr.reconstruction.optimize_linear(cost, max_iter=max_iter)
-    param_cache = cost.fwd_model.ramp.param_cache
-if order >= 1:
-    offset, ramp = param_cache[0][0], (param_cache[1][0], param_cache[2][0])
-elif order == 0:
-    offset, ramp = param_cache[0][0], (0, 0)
-else:
-    offset, ramp = 0, (0, 0)
-mag_data_buffer = mag_data_rec.copy()
-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_hdf5(mag_name + '.hdf5', overwrite=True)
-
-# Plot stuff:
-mag_data_rec.quiver_plot('Reconstructed Distribution', ar_dens=ar_dens)
-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)
-title = 'Reconstructed Phase'
-if order >= 0:
-    print('offset:', offset)
-    title += ', fitted Offset: {:.2g} [rad]'.format(offset)
-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')
-plt.show()