From 15b9227c7cc09cb9c88c52920b95a804b33a9138 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Sun, 24 Aug 2014 23:16:39 +0200
Subject: [PATCH] Added index converter class to convert matrix indices into 3D
coordinates
---
pyramid/converter.py | 68 +++++++++++++++++++++++++++++++++
pyramid/regularisator.py | 16 ++++++++
scripts/collaborations/zi_an.py | 22 ++++++-----
3 files changed, 96 insertions(+), 10 deletions(-)
create mode 100644 pyramid/converter.py
diff --git a/pyramid/converter.py b/pyramid/converter.py
new file mode 100644
index 0000000..3303227
--- /dev/null
+++ b/pyramid/converter.py
@@ -0,0 +1,68 @@
+# -*- 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))
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index b3621b0..5ade7c6 100644
--- a/pyramid/regularisator.py
+++ b/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))
diff --git a/scripts/collaborations/zi_an.py b/scripts/collaborations/zi_an.py
index 85df772..3d9b03e 100644
--- a/scripts/collaborations/zi_an.py
+++ b/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)
--
GitLab