Added index converter class to convert matrix indices into 3D coordinates

pyramid/converter.py

+# -*- coding: utf-8 -*-
+"""
+Created on Tue Aug 19 08:48:45 2014
+
+@author: Jan
+"""  # TODO: Docstring
+
+# TODO: put into other class
+# TODO: use 3 components (more complex)
+# TODO: take masks into account
+
+
+import numpy as np
+
+
+class IndexConverter3components(object):
+
+    def __init__(self, dim):
+        self.dim = dim
+        self.size_3d = np.prod(dim)
+        self.size_2d = dim[2]*dim[1]
+
+    def ind_to_coord(self, ind):
+        m, remain = int(ind/self.size_3d), ind % self.size_3d
+        z, remain = int(remain/self.size_2d), remain%self.size_2d
+        y, x = int(remain/self.dim[2]), remain%self.dim[2]
+        coord = m, z, y, x
+        return coord
+
+    def coord_to_ind(self, coord):
+        z, y, x = coord
+        ind = [i*self.size_3d + z*self.size_2d + y*self.dim[2] + x for i in range(3)]
+        return ind
+
+    def find_neighbour_ind(self, coord):
+        z, y, x = coord
+        t, d = (z-1, y, x), (z+1, y, x)  # t: top, d: down
+        f, b = (z, y-1, x), (z, y+1, x)  # f: front, b: back
+        l, r = (z, y, x-1), (z, y, x+1)  # l: left, r: right
+        neighbours = [self.coord_to_ind(i) for i in [t, d, f, b, l, r]]
+        return np.reshape(np.swapaxes(neighbours, 0, 1), (3, 3, 2))
+
+
+class IndexConverter(object):
+
+    def __init__(self, dim):
+        self.dim = dim
+        self.size_2d = dim[2]*dim[1]
+
+    def ind_to_coord(self, ind):
+        z, remain = int(ind/self.size_2d), ind%self.size_2d
+        y, x = int(remain/self.dim[2]), remain%self.dim[2]
+        coord = z, y, x
+        return coord
+
+    def coord_to_ind(self, coord):
+        z, y, x = coord
+        ind = z*self.size_2d + y*self.dim[2] + x
+        return ind
+
+    def find_neighbour_ind(self, coord):
+        z, y, x = coord
+        t, d = (z-1, y, x), (z+1, y, x)  # t: top, d: down
+        f, b = (z, y-1, x), (z, y+1, x)  # f: front, b: back
+        l, r = (z, y, x-1), (z, y, x+1)  # l: left, r: right
+        neighbours = [self.coord_to_ind(i) for i in [t, d, f, b, l, r]]
+        return neighbours
+        return np.reshape(neighbours, (3, 2))

pyramid/regularisator.py

@@ -66,6 +66,22 @@ class FirstOrderRegularisator(Regularisator):
     # TODO: Docstring!
 
     def __init__(self, lam, size, x_a=None):
+
         Sa_inv = lam * eye(size)
+
+
+
+
+
+        csr_matrix(coo_matrix((np.tile(data, dim_rot), (rows, columns)), shape=(size_2d, size_3d)))
+
+        term2 = []
+        for i in range(3):
+            component = mag_data[i, ...]
+            for j in range(3):
+                if component.shape[j] > 1:
+                    term2.append(np.diff(component, axis=j).reshape(-1))
+
+
         super(FirstOrderRegularisator, self).__init__(Sa_inv, x_a)
         self.LOG.debug('Created '+str(self))

scripts/collaborations/zi_an.py

@@ -46,21 +46,23 @@ order = 1
 if smoothed_pictures:
     dm3_2_mag = dm3.DM3(PATH+'Output333_hw512.dm3').image
     dm3_4_mag = dm3.DM3(PATH+'Output334_hw512.dm3').image
+    dm3_2_ele = dm3.DM3(PATH+'Output335.dm3').image
+    dm3_4_ele = dm3.DM3(PATH+'Output336.dm3').image
 else:
     dm3_2_mag = dm3.DM3(PATH+'07_0102mag60_q3_pha_01_sb280_sc512_vf3_med5.dm3').image
     dm3_4_mag = dm3.DM3(PATH+'18a_0102mag_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
-dm3_2_ele = dm3.DM3(PATH+'07_0102ele60_q3_pha_01_sb280_sc512_vf3_med5.dm3').image
-dm3_4_ele = dm3.DM3(PATH+'18a_0102ele_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
+    dm3_2_ele = dm3.DM3(PATH+'07_0102ele60_q3_pha_01_sb280_sc512_vf3_med5.dm3').image
+    dm3_4_ele = dm3.DM3(PATH+'18a_0102ele_ub140_62k_q3_pha_01_sb180_sc512_vf3_med5.dm3').image
 # Construct phase maps and masks
 phase_map_2 = PhaseMap(a, np.array(dm3_2_mag.resize(dim_small))+0.101)
 phase_map_2.display_combined(density=density, interpolation=inter)
 plt.savefig(PATH+'phase_map_2part.png')
-mask_2 = np.expand_dims(np.where(np.array(dm3_2_ele.resize(dim_small)) > threshold,
+mask_2 = np.expand_dims(np.where(np.array(dm3_2_ele.resize(dim_small)) >= threshold,
                                  True, False), axis=0)
 phase_map_4 = PhaseMap(a, np.array(dm3_4_mag.resize(dim_small))-2.546)
 phase_map_4.display_combined(density=density, interpolation=inter)
 plt.savefig(PATH+'phase_map_4part.png')
-mask_4 = np.expand_dims(np.where(np.array(dm3_4_ele.resize(dim_small)) > threshold,
+mask_4 = np.expand_dims(np.where(np.array(dm3_4_ele.resize(dim_small)) >= threshold,
                                  True, False), axis=0)
 # Reconstruct the magnetic distribution:
 tic = clock()
@@ -78,12 +80,12 @@ phase_map_rec_4.display_combined('Reconstr. Distribution', density=density, inte
 plt.savefig(PATH+'phase_map_4part_rec.png')
 # Plot the magnetization:
 axis = (mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
 plt.savefig(PATH+'mag_data_2part.png')
 axis = (mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot()
-axis.set_xlim(20, 45)
-axis.set_ylim(20, 45)
+#axis.set_xlim(20, 45)
+#axis.set_ylim(20, 45)
 plt.savefig(PATH+'mag_data_4part.png')
 # Display the Phase:
 phase_diff_2 = phase_map_rec_2-phase_map_2
@@ -96,13 +98,13 @@ plt.savefig(PATH+'phase_map_4part_diff.png')
 print 'Average difference (2 cubes):', np.average(phase_diff_2.phase)
 print 'Average difference (4 cubes):', np.average(phase_diff_4.phase)
 # Plot holographic contour maps with overlayed magnetic distributions:
-axis = phase_map_rec_2.display_holo('Magnetization Overlay', density=density, interpolation=inter)
+axis = phase_map_rec_2.display_holo('Magnetization Overlay', density=0.1, interpolation=inter)
 (mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot(axis=axis)
 axis = plt.gca()
 axis.set_xlim(20, 45)
 axis.set_ylim(20, 45)
 plt.savefig(PATH+'phase_map_2part_holo.png')
-axis = phase_map_rec_4.display_holo('Magnetization Overlay', density=density, interpolation=inter)
+axis = phase_map_rec_4.display_holo('Magnetization Overlay', density=0.1, interpolation=inter)
 (mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot(axis=axis)
 axis = plt.gca()
 axis.set_xlim(20, 45)