From 76697dc7b75f0eda68edac068e8c257114edb718 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Thu, 6 Jun 2013 10:06:56 +0200
Subject: [PATCH] Modified reconstructor and deleted now unused files

---
 pyramid/reconstructor.py             |  26 +++---
 pyramid/test/run_tests.py            |   4 +-
 scripts/invert2.py                   | 115 ---------------------------
 scripts/reconstruct_random_pixels.py |  37 ++-------
 scripts/reconstruct_test.py          |  73 -----------------
 5 files changed, 22 insertions(+), 233 deletions(-)
 delete mode 100644 scripts/invert2.py
 delete mode 100644 scripts/reconstruct_test.py

diff --git a/pyramid/reconstructor.py b/pyramid/reconstructor.py
index c1120f3..105804f 100644
--- a/pyramid/reconstructor.py
+++ b/pyramid/reconstructor.py
@@ -1,11 +1,6 @@
 # -*- coding: utf-8 -*-
-"""
-Created on Tue May 14 15:19:33 2013
+"""Reconstruct magnetic distributions with given phasemaps"""
 
-@author: Jan
-"""  # TODO: Docstring
-
-# TODO: Implement
 
 import numpy as np
 import pyramid.projector as pj
@@ -15,7 +10,16 @@ from scipy.optimize import leastsq
 
 
 def reconstruct_simple_lsqu(phase_map, mask, b_0):
-    # TODO: Docstring!
+    '''Reconstruct a magnetic distribution where the positions of the magnetized voxels are known
+    from a single phase_map using the least square method (only works for slice thickness = 1)
+    Arguments:
+        phase_map - a PhaseMap object, from which to reconstruct the magnetic distribution
+        mask      - a boolean matrix representing the positions of the magnetized voxels (3D)
+        b_0       - magnetic induction corresponding to a magnetization Mo in T (default: 1)
+    Returns:
+        the reconstructed magnetic distribution (as a MagData object)
+
+    '''
     # Read in parameters:
     y_m = phase_map.phase.reshape(-1)  # Measured phase map as a vector
     res = phase_map.res  # Resolution
@@ -24,14 +28,13 @@ def reconstruct_simple_lsqu(phase_map, mask, b_0):
     lam = 1e-6  # Regularisation parameter
     # Create empty MagData object for the reconstruction:
     mag_data_rec = MagData(res, (np.zeros(dim), np.zeros(dim), np.zeros(dim)))
-    ############################# FORWARD MODEL ###################################################
+
     # Function that returns the phase map for a magnetic configuration x:
     def F(x):
         mag_data_rec.set_vector(mask, x)
         phase = pm.phase_mag_real(res, pj.simple_axis_projection(mag_data_rec), 'slab', b_0)
         return phase.reshape(-1)
-    ############################# FORWARD MODEL ###################################################
-    ############################# RECONSTRUCTION ##################################################
+
     # Cost function which should be minimized:
     def J(x_i):
         y_i = F(x_i)
@@ -40,6 +43,5 @@ def reconstruct_simple_lsqu(phase_map, mask, b_0):
         return np.concatenate([term1, term2])
     # Reconstruct the magnetization components:
     x_rec, _ = leastsq(J, np.zeros(3*count))
-    ############################# RECONSTRUCTION ##################################################
     mag_data_rec.set_vector(mask, x_rec)
-    return mag_data_rec
\ No newline at end of file
+    return mag_data_rec
diff --git a/pyramid/test/run_tests.py b/pyramid/test/run_tests.py
index 16554d9..a3c6eb0 100644
--- a/pyramid/test/run_tests.py
+++ b/pyramid/test/run_tests.py
@@ -1,7 +1,5 @@
 # -*- coding: utf-8 -*-
-"""
-Unittests for pyramid.
-"""
+"""Unittests for pyramid."""
 
 import unittest
 from test_compliance import TestCaseCompliance
diff --git a/scripts/invert2.py b/scripts/invert2.py
deleted file mode 100644
index 5837248..0000000
--- a/scripts/invert2.py
+++ /dev/null
@@ -1,115 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Create random magnetic distributions."""
-
-import random as rnd
-import pdb, traceback, sys
-import numpy as np
-from numpy import pi
-import pylab
-from copy import deepcopy
-
-import pyramid.magcreator as mc
-from pyramid.magdata import MagData
-import pyramid.phasemapper as pm
-import pyramid.projector as pj
-import pyramid.holoimage as hi
-from pyramid.phasemap import PhaseMap
-from scipy.optimize import leastsq
-
-
-def create_random_distribution():
-    '''Calculate, display and save a random magnetic distribution to file.
-    Arguments:
-        None
-    Returns:
-        None
-
-    '''
-    # Input parameters:
-    count = 200
-    dim = (1, 32, 32)
-    b_0 = 1    # in T
-    res = 10 # in nm
-    rnd.seed(12)
-    # Create lists for magnetic objects:
-    mag_shape_list = np.zeros((count,) + dim)
-    beta_list      = np.zeros(count)
-    magnitude_list = np.zeros(count)
-    for i in range(count):
-        pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
-        mag_shape_list[i,...] = mc.Shapes.pixel(dim, pixel)
-        beta_list[i] = 2*pi*rnd.random()
-        magnitude_list[i] = rnd.random()
-    # Create magnetic distribution
-    magnitude = mc.create_mag_dist_comb(mag_shape_list, beta_list, magnitude_list)
-    mag_data = MagData(res, magnitude)
-#    mag_data.quiver_plot()
-#    mag_data.save_to_llg('output/mag_dist_random_pixel.txt')
-
-
-    magx, magy, magz = mag_data.magnitude
-    maskx, masky, maskz = magx != 0, magy != 0, magz !=0
-    x_t = np.concatenate([magx[maskx], magy[masky], magz[maskz]])
-    print "x_t", x_t
-
-    def F(x):
-        mag_data_temp = MagData(deepcopy(mag_data.res), deepcopy(mag_data.magnitude))
-        magx, magy, magz = mag_data_temp.magnitude
-        maskx, masky, maskz = magx != 0, magy != 0, magz !=0
-#        print maskx.sum() + masky.sum() + maskz.sum()
-        assert len(x) == maskx.sum() + masky.sum() + maskz.sum()
-        magx[maskx] = x[:maskx.sum()]
-        magy[masky] = x[maskx.sum():maskx.sum() + masky.sum()]
-        magz[maskz] = x[maskx.sum() + masky.sum():]
-        projection = pj.simple_axis_projection(mag_data_temp)
-        phase_map_slab = PhaseMap(res, pm.phase_mag_real_ANGLE(res, projection, 'slab', b_0))
-        return phase_map_slab.phase.reshape(-1)
-
-    y_m = F(x_t)
-    print "y_m", y_m
-
-    lam = 1e-6
-
-    def J(x_i):
-#        print "x_i", x_i
-        y_i = F(x_i)
-#        dd1 = np.zeros(mx.shape)
-#        dd2 = np.zeros(mx.shape)
-        term1 = (y_i - y_m)
-        term2 = lam * x_i
-#        dd1[:, :-1, :] += np.diff(mx, axis=1)
-#        dd1[:, -1, :] += np.diff(mx, axis=1)[:, -1, :]
-#        dd1[:, :, :-1] += np.diff(mx, axis=2)
-#        dd1[:, :, -1] += np.diff(mx, axis=2)[:, :, -1]
-
-#        dd2[:, :-1, :] += np.diff(my, axis=1)
-#        dd2[:, -1, :] += np.diff(my, axis=1)[:, -1, :]
-#        dd2[:, :, :-1] += np.diff(my, axis=2)
-#        dd2[:, :, -1] += np.diff(my, axis=2)[:, :, -1]
-
-#        result = np.concatenate([term1, np.sqrt(abs(dd1.reshape(-1))), np.sqrt(abs(dd2.reshape(-1)))])
-#        result = np.concatenate([term1, np.sqrt(abs(dd1.reshape(-1)))])
-#        print result
-        return np.concatenate([term1, term2])
-
-    x_f, _ = leastsq(J, np.zeros(x_t.shape))
-    y_f = F(x_f)
-#    print "y_m", y_m
-#    print "y_f", y_f
-#    print "dy", y_f - y_m
-#    print "x_t", x_t
-#    print "x_f", x_f
-#    print "dx", x_f - x_t
-#    pylab.show()
-    projection = pj.simple_axis_projection(mag_data)
-    phase_map  = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
-    hi.display(hi.holo_image(phase_map, 10))
-
-
-if __name__ == "__main__":
-    try:
-        create_random_distribution()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)
diff --git a/scripts/reconstruct_random_pixels.py b/scripts/reconstruct_random_pixels.py
index 4c42eb9..15a20d4 100644
--- a/scripts/reconstruct_random_pixels.py
+++ b/scripts/reconstruct_random_pixels.py
@@ -12,7 +12,7 @@ import pyramid.phasemapper as pm
 import pyramid.projector as pj
 import pyramid.holoimage as hi
 from pyramid.phasemap import PhaseMap
-from scipy.optimize import leastsq
+import pyramid.reconstructor as rc
 
 
 def reconstruct_random_distribution():
@@ -26,10 +26,11 @@ def reconstruct_random_distribution():
     # Input parameters:
     n_pixel = 10
     dim = (32, 32, 32)
-    b_0 = 1    # in T
+    b_0 = 1 # in T
     res = 10.0 # in nm
     rnd.seed(12)
     threshold = 0
+
     # Create lists for magnetic objects:
     mag_shape_list = np.zeros((n_pixel,) + dim)
     beta_list      = np.zeros(n_pixel)
@@ -52,34 +53,10 @@ def reconstruct_random_distribution():
     x_mask = abs(x_mag) > threshold
     y_mask = abs(y_mag) > threshold
     mask = np.logical_or(np.logical_or(x_mask, y_mask), z_mask)
-    # True values for the magnetisation informations, condensed into one vector:
-    x_t = mag_data.get_vector(mask)
-    # Create empty MagData object for the reconstruction:
-    mag_data_rec = MagData(res, (np.zeros(dim), np.zeros(dim), np.zeros(dim)))
-    ############################# FORWARD MODEL ###################################################
-    # Function that returns the phase map for a magnetic configuration x:
-    def F(x):
-        mag_data_rec.set_vector(mask, x)
-        phase = pm.phase_mag_real(res, pj.simple_axis_projection(mag_data_rec), 'slab', b_0)
-        return phase.reshape(-1)
-    ############################# FORWARD MODEL ###################################################
-    # Get a vector containing the measured phase at the specified places:
-    y_m = F(x_t)
-    print "y_m", y_m
-    ############################# RECONSTRUCTION ##################################################
-    lam = 1e-6  # Regularisation parameter
-    # Cost function which should be minimized:
-    def J(x_i):
-        y_i = F(x_i)
-        term1 = (y_i - y_m)
-        term2 = lam * x_i
-        return np.concatenate([term1, term2])
-    # Reconstruct the magnetization components:
-    x_f, _ = leastsq(J, np.zeros(x_t.shape))
-    ############################# RECONSTRUCTION ##################################################
-    # Save the reconstructed values in the MagData object:
-    y_f = F(x_f)
-    print "y_f", y_f
+    
+    # Reconstruct the magnetic distribution:
+    mag_data_rec = rc.reconstruct_simple_lsqu(phase_map, mask, b_0)
+
     # Display the reconstructed phase map and holography image:
     projection_rec = pj.simple_axis_projection(mag_data_rec)
     phase_map_rec = PhaseMap(res, pm.phase_mag_real(res, projection_rec, 'slab', b_0))
diff --git a/scripts/reconstruct_test.py b/scripts/reconstruct_test.py
deleted file mode 100644
index 4d1b67b..0000000
--- a/scripts/reconstruct_test.py
+++ /dev/null
@@ -1,73 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Create random magnetic distributions."""
-
-import random as rnd
-import pdb, traceback, sys
-import numpy as np
-from numpy import pi
-
-import pyramid.magcreator as mc
-from pyramid.magdata import MagData
-import pyramid.phasemapper as pm
-import pyramid.projector as pj
-import pyramid.holoimage as hi
-from pyramid.phasemap import PhaseMap
-import pyramid.reconstructor as rc
-
-
-def reconstruct_random_distribution():
-    '''Calculate and reconstruct a random magnetic distribution.
-    Arguments:
-        None
-    Returns:
-        None
-
-    '''
-    # Input parameters:
-    n_pixel = 10
-    dim = (32, 32, 32)
-    b_0 = 1    # in T
-    res = 10.0 # in nm
-    rnd.seed(12)
-    threshold = 0
-
-    # Create lists for magnetic objects:
-    mag_shape_list = np.zeros((n_pixel,) + dim)
-    beta_list      = np.zeros(n_pixel)
-    magnitude_list = np.zeros(n_pixel)
-    for i in range(n_pixel):
-        pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
-        mag_shape_list[i,...] = mc.Shapes.pixel(dim, pixel)
-        beta_list[i] = 2*pi*rnd.random()
-        magnitude_list[i] = rnd.random()
-    # Create magnetic distribution:
-    magnitude = mc.create_mag_dist_comb(mag_shape_list, beta_list, magnitude_list)
-    mag_data = MagData(res, magnitude)
-    # Display phase map and holography image:
-    projection = pj.simple_axis_projection(mag_data)
-    phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0))
-    hi.display_combined(phase_map, 10, 'Generated Distribution')
-
-    # Get the locations of the magnetized pixels (mask):
-    z_mag, y_mag, x_mag = mag_data.magnitude
-    z_mask = abs(z_mag) > threshold
-    x_mask = abs(x_mag) > threshold
-    y_mask = abs(y_mag) > threshold
-    mask = np.logical_or(np.logical_or(x_mask, y_mask), z_mask)
-    
-    # Reconstruct the magnetic distribution:
-    mag_data_rec = rc.reconstruct_simple_lsqu(phase_map, mask, b_0)
-
-    # Display the reconstructed phase map and holography image:
-    projection_rec = pj.simple_axis_projection(mag_data_rec)
-    phase_map_rec = PhaseMap(res, pm.phase_mag_real(res, projection_rec, 'slab', b_0))
-    hi.display_combined(phase_map_rec, 10, 'Reconstructed Distribution')
-
-
-if __name__ == "__main__":
-    try:
-        reconstruct_random_distribution()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)
-- 
GitLab