Skip to content
Snippets Groups Projects
Commit 76697dc7 authored by Jan Caron's avatar Jan Caron
Browse files

Modified reconstructor and deleted now unused files

parent f99371b8
No related branches found
No related tags found
No related merge requests found
# -*- 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
# -*- coding: utf-8 -*-
"""
Unittests for pyramid.
"""
"""Unittests for pyramid."""
import unittest
from test_compliance import TestCaseCompliance
......
# -*- 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)
......@@ -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))
......
# -*- 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)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment