Implemented projection with tilt along one axis and unittests for all modules!

pyramid/holoimage.py

@@ -69,7 +69,8 @@ def holo_image(phase_map, density=1):
     phase_grad_y, phase_grad_x = np.gradient(phase, phase_map.res, phase_map.res)
     phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
     phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)
-    phase_magnitude = np.sin(phase_magnitude/phase_magnitude.max() * pi / 2)
+    if phase_magnitude.max() != 0:
+        phase_magnitude = np.sin(phase_magnitude/phase_magnitude.max() * pi / 2)
     # Color code the angle and create the holography image:
     rgba = HOLO_CMAP(phase_angle)
     rgb = (255.999 * img_holo.T * phase_magnitude.T * rgba[:, :, :3].T).T.astype(np.uint8)

pyramid/magdata.py

@@ -104,6 +104,22 @@ class MagData:
             self.magnitude = magnitude
             self.dim = dim
 
+    def get_mask(self, threshold=0):
+        '''Mask all pixels where the amplitude of the magnetization lies above `threshold`.
+
+        Parameters
+        ----------
+        threshold : float, optional
+            A pixel only gets masked, if it lies above this threshold . The default is 0.
+
+        Returns
+        -------
+        mask : :class:`~numpy.ndarray` (N=3, boolean)
+            Mask of the pixels where the amplitude of the magnetization lies above `threshold`.
+
+        '''
+        return np.sqrt(np.sum(np.array(self.magnitude)**2, axis=0)) > threshold
+
     def get_vector(self, mask):
         '''Returns the magnetic components arranged in a vector, specified by a mask.
 
@@ -145,22 +161,6 @@ class MagData:
         self.magnitude[1][mask] = vector[count:2*count]  # y-component
         self.magnitude[0][mask] = vector[2*count:]  # z-component
 
-    def get_mask(self, threshold=0):
-        '''Mask all pixels where the amplitude of the magnetization lies above `threshold`.
-
-        Parameters
-        ----------
-        threshold : float, optional
-            A pixel only gets masked, if it lies above this threshold . The default is 0.
-
-        Returns
-        -------
-        mask : :class:`~numpy.ndarray` (N=3, boolean)
-            Mask of the pixels where the amplitude of the magnetization lies above `threshold`.
-
-        '''
-        return np.sqrt(np.sum(np.array(self.magnitude)**2, axis=0)) > threshold
-
     def scale_down(self, n=1):
         '''Scale down the magnetic distribution by averaging over two pixels along each axis.
 
@@ -355,8 +355,8 @@ class MagData:
             axis = fig.add_subplot(1, 1, 1, aspect='equal')
         axis.quiver(mag_slice_u, mag_slice_v, pivot='middle', angles='xy', scale_units='xy',
                    scale=1, headwidth=6, headlength=7)
-        axis.set_xlim(0, np.shape(mag_slice_u)[1])
-        axis.set_ylim(0, np.shape(mag_slice_u)[0])
+        axis.set_xlim(-1, np.shape(mag_slice_u)[1])
+        axis.set_ylim(-1, np.shape(mag_slice_u)[0])
         axis.set_title(title, fontsize=18)
         axis.set_xlabel(u_label, fontsize=15)
         axis.set_ylabel(v_label, fontsize=15)

pyramid/phasemap.py

@@ -39,7 +39,7 @@ class PhaseMap:
         always in `rad`.
 
     '''
-    
+
     UNITDICT = {'rad': 1E0,
                 'mrad': 1E3,
                 'µrad': 1E6}
@@ -76,7 +76,7 @@ class PhaseMap:
 
         Parameters
         ----------
-        unit : {'rad', 'mrad', 'µrad'}, optional
+        unit : {'rad', 'mrad'}, optional
             Set the unit of the phase map. This is important for the :func:`~.display` function,
             because the phase is scaled accordingly. Does not change the phase itself, which is
             always in `rad`.
@@ -86,6 +86,7 @@ class PhaseMap:
         None
 
         '''
+        assert unit in ['rad', 'mrad']
         self.unit = unit
 
     @classmethod
@@ -105,7 +106,7 @@ class PhaseMap:
         '''
         with open(filename, 'r') as phase_file:
             phase_file.readline()  # Headerline is not used
-            res = int(phase_file.readline()[13:-4])
+            res = float(phase_file.readline()[13:-4])
             phase = np.loadtxt(filename, delimiter='\t', skiprows=2)
         return PhaseMap(res, phase)
 

pyramid/phasemapper.py

@@ -16,6 +16,8 @@ from numpy import pi
 
 import pyramid.numcore as nc
 
+from scipy.signal import fftconvolve
+
 
 PHI_0 = -2067.83    # magnetic flux in T*nm²
 H_BAR = 6.626E-34  # Planck constant in J*s
@@ -145,6 +147,144 @@ def phase_mag_real(res, projection, method='discs', b_0=1, jacobi=None):
     return phase
 
 
+def phase_mag_real_conv(res, projection, method='disc', b_0=1, jacobi=None):
+    '''Calculate the magnetic phase from magnetization data (real space approach).
+
+    Parameters
+    ----------
+    res : float
+        The resolution of the grid (grid spacing) in nm.
+    projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
+        The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
+        of the magnetization and the thickness projection for the resulting 2D-grid.
+    method : {'disc', 'slab'}, optional
+        Specifies the elemental geometry to use for the pixel field.
+        The default is 'disc', because of the smaller computational overhead.
+    b_0 : float, optional
+        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+        The default is 1.
+    jacobi : :class:`~numpy.ndarray` (N=2), optional
+        Specifies the matrix in which to save the jacobi matrix. The jacobi matrix will not be
+        calculated, if no matrix is specified (default), resulting in a faster computation.
+
+    Returns
+    -------
+    phase : :class:`~numpy.ndarray` (N=2)
+        The phase as a 2-dimensional array.
+
+    '''  # TODO: Docstring!!!
+
+    # Function for creating the lookup-tables:
+    def phi_lookup(method, n, m, res, b_0):
+        if method == 'slab':
+            def F_h(n, m):
+                a = np.log(res**2 * (n**2 + m**2))
+                b = np.arctan(n / m)
+                return n*a - 2*n + 2*m*b
+            return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
+                                - F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
+        elif method == 'disc':
+            in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
+            return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
+
+    # Process input parameters:
+    v_dim, u_dim = np.shape(projection[0])
+    v_mag, u_mag = projection[:-1]
+    coeff = -b_0 * res**2 / (2*PHI_0)
+
+    # Create lookup-tables for the phase of one pixel:
+    u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
+    v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
+    uu, vv = np.meshgrid(u, v)
+    u_phi = phi_lookup(method, uu, vv, res, b_0)
+    v_phi = phi_lookup(method, vv, uu, res, b_0)
+
+    # Return the phase:
+    return fftconvolve(u_mag, u_phi, 'same') - fftconvolve(v_mag, v_phi, 'same')
+
+#    def fftconvolve(in1, in2, mode="full"):
+#    """Convolve two N-dimensional arrays using FFT.
+#
+#    Convolve `in1` and `in2` using the fast Fourier transform method, with
+#    the output size determined by the `mode` argument.
+#
+#    This is generally much faster than `convolve` for large arrays (n > ~500),
+#    but can be slower when only a few output values are needed, and can only
+#    output float arrays (int or object array inputs will be cast to float).
+#
+#    Parameters
+#    ----------
+#    in1 : array_like
+#        First input.
+#    in2 : array_like
+#        Second input. Should have the same number of dimensions as `in1`;
+#        if sizes of `in1` and `in2` are not equal then `in1` has to be the
+#        larger array.
+#    mode : str {'full', 'valid', 'same'}, optional
+#        A string indicating the size of the output:
+#
+#        ``full``
+#           The output is the full discrete linear convolution
+#           of the inputs. (Default)
+#        ``valid``
+#           The output consists only of those elements that do not
+#           rely on the zero-padding.
+#        ``same``
+#           The output is the same size as `in1`, centered
+#           with respect to the 'full' output.
+#
+#    Returns
+#    -------
+#    out : array
+#        An N-dimensional array containing a subset of the discrete linear
+#        convolution of `in1` with `in2`.
+#
+#    """
+#    in1 = asarray(in1)
+#    in2 = asarray(in2)
+#
+#    if rank(in1) == rank(in2) == 0:  # scalar inputs
+#        return in1 * in2
+#    elif not in1.ndim == in2.ndim:
+#        raise ValueError("in1 and in2 should have the same rank")
+#    elif in1.size == 0 or in2.size == 0:  # empty arrays
+#        return array([])
+#
+#    s1 = array(in1.shape)
+#    s2 = array(in2.shape)
+#    complex_result = (np.issubdtype(in1.dtype, np.complex) or
+#                      np.issubdtype(in2.dtype, np.complex))
+#    size = s1 + s2 - 1
+#
+#    if mode == "valid":
+#        for d1, d2 in zip(s1, s2):
+#            if not d1 >= d2:
+#                warnings.warn("in1 should have at least as many items as in2 in "
+#                              "every dimension for 'valid' mode.  In scipy "
+#                              "0.13.0 this will raise an error",
+#                              DeprecationWarning)
+#
+#    # Always use 2**n-sized FFT
+#    fsize = 2 ** np.ceil(np.log2(size)).astype(int)
+#    print('fsize =', fsize)
+#    print('s1 =', s1)
+#    print('s2 =', s2)
+#    fslice = tuple([slice(0, int(sz)) for sz in size])
+#    if not complex_result:
+#        ret = irfftn(rfftn(in1, fsize) *
+#                     rfftn(in2, fsize), fsize)[fslice].copy()
+#        ret = ret.real
+#    else:
+#        ret = ifftn(fftn(in1, fsize) * fftn(in2, fsize))[fslice].copy()
+#
+#    if mode == "full":
+#        return ret
+#    elif mode == "same":
+#        return _centered(ret, s1)
+#    elif mode == "valid":
+#        return _centered(ret, abs(s1 - s2) + 1)
+
+
 def phase_elec(res, projection, v_0=1, v_acc=30000):
     '''Calculate the electric phase from magnetization distributions.
 

pyramid/projector.py

@@ -8,6 +8,12 @@ directions with the use of transfer functions (work in progress).
 
 """
 
+import time
+
+import numpy as np
+from numpy import pi
+
+import itertools
 
 from pyramid.magdata import MagData
 
@@ -49,3 +55,84 @@ def simple_axis_projection(mag_data, axis='z', threshold=0):
                       mag_data.magnitude[1].sum(2),  # y_mag -> u_mag
                       mag_data.get_mask(threshold).sum(2))  # thickness profile
     return projection
+
+
+def single_tilt_projection(mag_data, tilt=0, threshold=0):
+    # TODO: Docstring!!!
+    assert isinstance(mag_data, MagData), 'Parameter mag_data has to be a MagData object!'
+#    assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as string)!'
+
+    # Set starting variables:
+    dim = (mag_data.dim[0], mag_data.dim[2])
+    z_mag, y_mag, x_mag = mag_data.magnitude
+    mask = mag_data.get_mask()
+    projection = (np.zeros((mag_data.dim[1], mag_data.dim[2])),
+                  np.zeros((mag_data.dim[1], mag_data.dim[2])),
+                  np.zeros((mag_data.dim[1], mag_data.dim[2])))
+
+    def get_position(p, m, b, size):
+        x, y = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
+        return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
+
+    def get_impact(pos, r, size):
+        return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
+                if 0 <= x < size]
+
+    def get_weight(delta, rho):
+        a, b = delta-rho, delta+rho
+        if a >= 1 or b <= -1:  # TODO: Should not be necessary!
+            print 'Outside of bounds:', delta
+            return 0
+        # Upper boundary:
+        if b >= 1:
+            w_b = 0.5
+        else:
+            w_b = (b*np.sqrt(1-b**2) + np.arctan(b/np.sqrt(1-b**2))) / pi
+        # Lower boundary:
+        if a <= -1:
+            w_a = -0.5
+        else:
+            w_a = (a*np.sqrt(1-a**2) + np.arctan(a/np.sqrt(1-a**2))) / pi
+        return w_b - w_a
+
+    # Creating coordinate list of all voxels:
+    xi = range(dim[0])
+    yj = range(dim[1])
+    ii, jj = np.meshgrid(xi, yj)
+    voxels = list(itertools.product(yj, xi))
+
+    # Calculate positions along the projected pixel coordinate system:
+    direction = (-np.cos(tilt), np.sin(tilt))
+    center = (dim[0]/2., dim[1]/2.)
+    m = direction[0]/(direction[1]+1E-30)
+    b = center[0] - m * center[1]
+    positions = get_position(voxels, m, b, dim[0])
+
+    # Calculate weights:
+    r = 1/np.sqrt(np.pi)  # radius of the voxel circle
+    rho = 0.5 / r  # TODO: ratio of radii
+    weights = []
+    for i, voxel in enumerate(voxels):
+        voxel_weights = []
+        impacts = get_impact(positions[i], r, dim[0])
+        for impact in impacts:
+            distance = np.abs(impact+0.5 - positions[i])
+            delta = distance / r
+            voxel_weights.append((impact, get_weight(delta, rho)))
+        weights.append((voxel, voxel_weights))
+
+    # Calculate projection with the calculated weights for the voxels:
+    for i, weight in enumerate(weights):
+        voxel = weights[i][0]
+        voxel_weights = weights[i][1]
+        for voxel_weight in voxel_weights:
+            pixel, weight = voxel_weight
+            # Component parallel to tilt axis (':' goes over all slices):
+            projection[0][:, pixel] += weight * y_mag[voxel[0], :, voxel[1]]
+            # Component perpendicular to tilt axis:
+            projection[1][:, pixel] += weight * (x_mag[voxel[0], :, voxel[1]]*np.cos(tilt)
+                                               + z_mag[voxel[0], :, voxel[1]]*np.sin(tilt))
+            # Thickness profile:
+            projection[2][:, pixel] += weight * mask[voxel[0], :, voxel[1]]
+
+    return projection

scripts/compare methods/compare_methods.py

@@ -32,18 +32,18 @@ def compare_methods():
     b_0 = 1.1    # in T
     res = 10.0  # in nm
     phi = pi/4
-    padding = 12
+    padding = 3
     density = 1
-    dim = (16, 128, 128)  # in px (z, y, x)
+    dim = (8, 512, 512)  # in px (z, y, x)
     # Create magnetic shape:
     geometry = 'slab'
     if geometry == 'slab':
-        center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
-        width = (dim[0]/2, dim[1]/2., dim[2]/2.)  # in px (z, y, x)
+        center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z,y,x) index starts at 0!
+        width = (dim[0]/2, dim[1]/2, dim[2]/2)  # in px (z, y, x)
         mag_shape = mc.Shapes.slab(dim, center, width)
         phase_ana = an.phase_mag_slab(dim, res, phi, center, width, b_0)
     elif geometry == 'disc':
-        center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
+        center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z,y,x) index starts at 0!
         radius = dim[1]/4  # in px
         height = dim[0]/2  # in px
         mag_shape = mc.Shapes.disc(dim, center, radius, height)
@@ -68,21 +68,35 @@ def compare_methods():
     start_time = time.time()
     phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc', b_0))
     print 'Time for real space approach (Disc):', time.time() - start_time
+    start_time = time.time()
+    phase_map_slab_conv = PhaseMap(res, pm.phase_mag_real_conv(res, projection, 'slab', b_0))
+    print 'Time for real space approach (Slab Convolve):', time.time() - start_time
+    start_time = time.time()
+    phase_map_disc_conv = PhaseMap(res, pm.phase_mag_real_conv(res, projection, 'disc', b_0))
+    print 'Time for real space approach (Disc Convolve):', time.time() - start_time
     # Display the combinated plots with phasemap and holography image:
     hi.display_combined(phase_map_ana, density, 'Analytic Solution')
     hi.display_combined(phase_map_fft, density, 'Fourier Space')
     hi.display_combined(phase_map_slab, density, 'Real Space (Slab)')
     hi.display_combined(phase_map_disc, density, 'Real Space (Disc)')
+    hi.display_combined(phase_map_slab_conv, density, 'Real Space (Slab Convolve)')
+    hi.display_combined(phase_map_disc_conv, density, 'Real Space (Disc Convolve)')
     # Plot differences to the analytic solution:
     phase_map_diff_fft = PhaseMap(res, phase_map_ana.phase-phase_map_fft.phase)
     phase_map_diff_slab = PhaseMap(res, phase_map_ana.phase-phase_map_slab.phase)
     phase_map_diff_disc = PhaseMap(res, phase_map_ana.phase-phase_map_disc.phase)
+    phase_map_diff_slab_conv = PhaseMap(res, phase_map_ana.phase-phase_map_slab_conv.phase)
+    phase_map_diff_disc_conv = PhaseMap(res, phase_map_ana.phase-phase_map_disc_conv.phase)
     RMS_fft = np.std(phase_map_diff_fft.phase)
     RMS_slab = np.std(phase_map_diff_slab.phase)
     RMS_disc = np.std(phase_map_diff_disc.phase)
+    RMS_slab_conv = np.std(phase_map_diff_slab_conv.phase)
+    RMS_disc_conv = np.std(phase_map_diff_disc_conv.phase)
     phase_map_diff_fft.display('Fourier Space (RMS = {:3.2e})'.format(RMS_fft))
     phase_map_diff_slab.display('Real Space (Slab) (RMS = {:3.2e})'.format(RMS_slab))
     phase_map_diff_disc.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_disc))
+    phase_map_diff_slab_conv.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_slab_conv))
+    phase_map_diff_disc_conv.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_disc_conv))
 
 
 if __name__ == "__main__":

scripts/compare methods/compare_pixel_fields.py

@@ -8,7 +8,7 @@ import traceback
 import sys
 import os
 
-from numpy import pi
+import numpy as np
 
 import pyramid.magcreator as mc
 import pyramid.projector as pj
@@ -30,20 +30,63 @@ def compare_pixel_fields():
         os.makedirs(directory)
     # Input parameters:
     res = 1.0  # in nm
-    phi = pi/2  # in rad
-    dim = (1, 101, 101)
+    phi = 0  # in rad
+    dim = (1, 32, 32)
     pixel = (0,  int(dim[1]/2),  int(dim[2]/2))
-    limit = 0.25
+    limit = 0.35
+
+    def get_fourier_kernel():
+        PHI_0 = 2067.83    # magnetic flux in T*nm²
+        b_0 = 1
+        coeff = - 1j * b_0 * res**2 / (2*PHI_0)
+        nyq = 1 / res  # nyquist frequency
+        f_u = np.linspace(0, nyq/2, dim[2]/2.+1)
+        f_v = np.linspace(-nyq/2, nyq/2, dim[1], endpoint=False)
+        f_uu, f_vv = np.meshgrid(f_u, f_v)
+        phase_fft = coeff * f_vv / (f_uu**2 + f_vv**2 + 1e-30)
+        # Transform to real space and revert padding:
+        phase_fft = np.fft.ifftshift(phase_fft, axes=0)
+        phase_fft_kernel = np.fft.fftshift(np.fft.irfft2(phase_fft), axes=(0, 1))
+        return phase_fft_kernel
+
     # Create magnetic data, project it, get the phase map and display the holography image:
     mag_data = MagData(res, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, pixel), phi))
     mag_data.save_to_llg(directory + '/mag_dist_single_pixel.txt')
     projection = pj.simple_axis_projection(mag_data)
-    phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'), 'mrad')
-    phase_map_slab.display('Phase of one Pixel (Slab)', limit=limit)
+    # Kernel of a disc in real space:
     phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc'), 'mrad')
     phase_map_disc.display('Phase of one Pixel (Disc)', limit=limit)
-    phase_map_diff = PhaseMap(res, phase_map_disc.phase - phase_map_slab.phase, 'mrad')
-    phase_map_diff.display('Phase difference of one Pixel (Disc - Slab)')
+    # Kernel of a slab in real space:
+    phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'), 'mrad')
+    phase_map_slab.display('Phase of one Pixel (Slab)', limit=limit)
+    # Kernel of the Fourier method:
+    phase_map_fft = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=0), 'mrad')
+    phase_map_fft.display('Phase of one Pixel (Fourier)', limit=limit)
+    # Kernel of the Fourier method, calculated directly:
+    phase_map_fft_kernel = PhaseMap(res, get_fourier_kernel(), 'mrad')
+    phase_map_fft_kernel.display('Phase of one Pixel (Fourier Kernel)', limit=limit)
+    # Kernel differences:
+    print 'Fourier Kernel, direct and indirect method are identical:', \
+          np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
+    phase_map_diff = PhaseMap(res, phase_map_fft.phase - phase_map_disc.phase, 'mrad')
+    phase_map_diff.display('Phase difference of one Pixel (Disc - Fourier)')
+
+    phase_inv_fft_r = np.real(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))
+    phase_map_inv_fft = PhaseMap(res, phase_inv_fft_r)
+    phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, real)')
+
+    phase_inv_fft_i = np.imag(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))
+    phase_map_inv_fft = PhaseMap(res, phase_inv_fft_i)
+    phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, imag)')
+
+    phase_inv_disc_r = np.real(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))
+    phase_map_inv_disc = PhaseMap(res, phase_inv_disc_r)
+    phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, real)')
+
+    phase_inv_disc_i = np.imag(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))
+    phase_map_inv_disc = PhaseMap(res, phase_inv_disc_i)
+    phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, imag)')
+
 
 if __name__ == "__main__":
     try:

scripts/create distributions/create_sample.py

@@ -21,7 +21,7 @@ def create_sample():
         None
 
     '''
-    
+
     directory = '../../output/magnetic distributions'
     if not os.path.exists(directory):
         os.makedirs(directory)

scripts/gui/create_logo.py

@@ -7,13 +7,17 @@
 #
 # WARNING! All changes made in this file will be lost!
 
+
 from PyQt4 import QtCore, QtGui
+from matplotlibwidget import MatplotlibWidget
+
 
 try:
     _fromUtf8 = QtCore.QString.fromUtf8
 except AttributeError:
     _fromUtf8 = lambda s: s
 
+
 class Ui_CreateLogoWidget(object):
     def setupUi(self, CreateLogoWidget):
         CreateLogoWidget.setObjectName(_fromUtf8("CreateLogoWidget"))
@@ -64,5 +68,3 @@ class Ui_CreateLogoWidget(object):
         self.xLabel.setText(QtGui.QApplication.translate("CreateLogoWidget", "X [px] :", None, QtGui.QApplication.UnicodeUTF8))
         self.yLabel.setText(QtGui.QApplication.translate("CreateLogoWidget", "Y [px] :", None, QtGui.QApplication.UnicodeUTF8))
         self.logoPushButton.setText(QtGui.QApplication.translate("CreateLogoWidget", "Create Logo", None, QtGui.QApplication.UnicodeUTF8))
-
-from matplotlibwidget import MatplotlibWidget

scripts/interactive_setup.py

@@ -11,4 +11,4 @@ from pyramid.magdata import MagData
 from pyramid.phasemap import PhaseMap
 
 import os
-os.chdir('C:\Users\Jan\Home\PhD Thesis\Pyramid\output')
+os.chdir('C:\\Users\\Jan\\Home\\PhD Thesis\\Pyramid\\tests')

scripts/kernel_comparison.py

+# -*- coding: utf-8 -*-
+"""
+Created on Thu Sep 12 13:22:43 2013
+
+@author: Jan
+"""
+
+from pylab import *
+
+
+f = 1
+a, b = f * 32, f * 64
+f_u = np.linspace(0, 1./3, a)
+f_v = np.linspace(-1./6., 1./6., b, endpoint=False)
+f_uu, f_vv = np.meshgrid(f_u, f_v)
+phase_fft = 1j * f_vv / (f_uu**2 + f_vv**2 + 1e-30)
+# Transform to real space and revert padding:
+phase_fft = np.fft.ifftshift(phase_fft, axes=0)
+ss = np.fft.irfft2(phase_fft)
+print ss
+pcolormesh(np.fft.fftshift(np.fft.fftshift(ss, axes=1), axes=0), cmap='RdBu')
+colorbar()
+show()

scripts/test_arb_projection.py

+# -*- coding: utf-8 -*-
+"""
+Created on Tue Sep 03 12:55:40 2013
+
+@author: Jan
+"""
+
+'''2D Case Bresenham's line algorithm'''
+
+import numpy as np
+from numpy import pi
+import itertools
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MaxNLocator, NullLocator
+
+
+dim = (4, 4)
+offset = (dim[0]/2., dim[1]/2.)
+phi = 1.000000001*pi/4
+field = np.zeros(dim)
+
+
+def get_position(p, m, b, size):
+    x, y = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
+    return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
+
+
+def get_impact(pos, r, size):
+    return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int) if 0 <= x < size]
+
+
+def get_weight(d, r):
+    return (1 - d/r) * (d < r)
+
+
+direction = (-np.cos(phi), np.sin(phi))
+xi = range(dim[0])
+yj = range(dim[1])
+ii, jj = np.meshgrid(xi, yj)
+
+r = 1/np.sqrt(np.pi)
+
+m = direction[0]/direction[1]
+b = offset[0] - m * offset[1]
+
+voxels = list(itertools.product(yj, xi))
+
+positions = get_position(voxels, m, b, dim[0])
+
+weights = []
+for i, voxel in enumerate(voxels):
+    voxel_weights = []
+    impacts = get_impact(positions[i], r, dim[0])
+    for impact in impacts:
+        distance = np.abs(impact+0.5 - positions[i])
+        voxel_weights.append((impact, get_weight(distance, r)))
+    weights.append((voxel, voxel_weights))
+
+pixels = np.floor(positions).astype(int)
+
+pixel_hits = zip(set(pixels), [list(pixels).count(i) for i in set(pixels)])
+
+
+def Y(x):
+    return direction[0]/direction[1] * (x - offset[1]) + offset[0]
+
+
+def Y_perp(x):
+    return - direction[1]/direction[0] * (x - offset[1]) + offset[0]
+
+
+def X(y):
+    return direction[1]/direction[0] * (y - offset[0]) + offset[1]
+
+
+x = np.linspace(-2, 2, 10)
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1, aspect='equal')
+axis.pcolormesh(field, cmap='PuRd')
+axis.grid(which='both', color='k', linestyle='-')
+x = np.linspace(0, dim[1])
+y = Y(x)
+y_perp = Y_perp(x)
+axis.plot(x, y, '-r', linewidth=2)
+axis.plot(x, y_perp, '-g', linewidth=2)
+axis.set_xlim(0, dim[1])
+axis.set_ylim(0, dim[0])
+axis.xaxis.set_major_locator(MaxNLocator(nbins=dim[1], integer=True))
+axis.yaxis.set_major_locator(MaxNLocator(nbins=dim[0], integer=True))
+
+for i, p in enumerate(voxels):
+    if 0 <= positions[i] < dim[0]:
+        color = 'k'
+    else:
+        color = 'r'
+    plt.annotate('{:1.1f}'.format(positions[i]), p, size=8, color=color)
+
+plt.show()
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1, aspect='equal')
+axis.scatter(positions, 0.5*np.ones_like(positions))
+axis.grid(which='both', color='k', linestyle='-')
+axis.vlines((0, dim[0]), 0, 1, colors='r', linewidth=3)
+axis.set_xlim(-int(0.3*dim[0]), int(1.3*dim[0]))
+axis.set_ylim(0, 1)
+axis.xaxis.set_major_locator(MaxNLocator(nbins=int(1.6*dim[0]), integer=True))
+axis.yaxis.set_major_locator(NullLocator())
+
+for i, px in enumerate(pixel_hits):
+    plt.annotate('{:d}'.format(px[1]), (px[0], 0.6), size=8)
+
+plt.show()

scripts/test_methods.py

@@ -24,43 +24,46 @@ def compare_methods():
         None
 
     '''
-#    # Input parameters:
-#    res = 10.0  # in nm
-#    phi = pi/4
-#    density = 1
-#    dim = (64, 64, 64)  # in px (z, y, x)
-#    # Create magnetic shape:
-#    geometry = 'sphere'
-#    if geometry == 'slab':
-#        center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
-#        width = (dim[0]/2, dim[1]/2., dim[2]/2.)  # in px (z, y, x)
-#        mag_shape = mc.Shapes.slab(dim, center, width)
-#    elif geometry == 'disc':
-#        center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
-#        radius = dim[1]/4  # in px
-#        height = dim[0]/2  # in px
-#        mag_shape = mc.Shapes.disc(dim, center, radius, height)
-#    elif geometry == 'sphere':
-#        center = (32, 32, 32)  # in px (z, y, x) index starts with 0!
-#        radius = 16  # in px
-#        mag_shape = mc.Shapes.sphere(dim, center, radius)
-#    # Project the magnetization data:
-#    mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
+    # Input parameters:
+    res = 10.0  # in nm
+    phi = 0
+    density = 1
+    dim = (128, 128, 128)  # in px (z, y, x)
+    # Create magnetic shape:
+    geometry = 'slab'
+    if geometry == 'slab':
+        center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
+        width = (dim[0]/2, dim[1]/2., dim[2]/2.)  # in px (z, y, x)
+        mag_shape = mc.Shapes.slab(dim, center, width)
+    elif geometry == 'disc':
+        center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
+        radius = dim[1]/4  # in px
+        height = dim[0]/2  # in px
+        mag_shape = mc.Shapes.disc(dim, center, radius, height)
+    elif geometry == 'sphere':
+        center = (dim[0]/2, dim[1]/2, dim[2]/2)  # in px (z, y, x) index starts with 0!
+        radius = dim[0]/4  # in px
+        mag_shape = mc.Shapes.sphere(dim, center, radius)
+    # Project the magnetization data:
+    mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi, theta=0))
 
-    density = 0.3
+#    density = 0.3
+#    mag_data = MagData.load_from_llg('../output/magnetic distributions/mag_dist_sphere.txt')
 
-    mag_data = MagData.load_from_llg('../output/magnetic distributions/mag_dist_sphere.txt')
     res = mag_data.res
     mag_data.quiver_plot()
 #    mag_data.quiver_plot3d()
-    projection = pj.simple_axis_projection(mag_data)
+    import time
+    start = time.time()
+    projection = pj.single_tilt_projection(mag_data, tilt=pi/4)
+    print 'Total projection time:', time.time() - start
     # Construct phase maps:
     phase_map_mag = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=1))
-    phase_map_elec = PhaseMap(res, pm.phase_elec(res, projection, v_0 = 3))
+    phase_map_elec = PhaseMap(res, pm.phase_elec(res, projection, v_0=3))
     # Display the combinated plots with phasemap and holography image:
     hi.display_combined(phase_map_mag, density, title='Magnetic Phase')
     hi.display_combined(phase_map_elec, density, title='Electric Phase')
-    
+
     phase_map = PhaseMap(res, phase_map_mag.phase+phase_map_elec.phase)
     hi.display_combined(phase_map, density)
 

tests/test_analytic.py

@@ -2,32 +2,57 @@
 """Testcase for the analytic module."""
 
 
+import os
 import unittest
+import numpy as np
+from numpy import pi
+
 import pyramid.analytic as an
 
 
 class TestCaseAnalytic(unittest.TestCase):
 
     def setUp(self):
-        pass
+        self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_analytic/')
+        self.dim = (4, 4, 4)
+        self.res = 10.0
+        self.phi = pi/4
+        self.center = (self.dim[0]/2-0.5, self.dim[1]/2-0.5, self.dim[2]/2-0.5)
+        self.radius = self.dim[2]/4
 
     def tearDown(self):
-        pass
-
-    def test_template(self):
-        pass
+        self.path = None
+        self.dim = None
+        self.res = None
+        self.phi = None
+        self.center = None
+        self.radius = None
 
     def test_phase_mag_slab(self):
-        pass
+        width = (self.dim[0]/2, self.dim[1]/2, self.dim[2]/2)
+        phase = an.phase_mag_slab(self.dim, self.res, self.phi, self.center, width)
+        reference = np.load(os.path.join(self.path, 'ref_phase_slab.npy'))
+        np.testing.assert_equal(phase, reference, 'Unexpected behavior in phase_mag_slab()')
 
     def test_phase_mag_disc(self):
-        pass
+        radius = self.dim[2]/4
+        height = self.dim[2]/2
+        phase = an.phase_mag_disc(self.dim, self.res, self.phi, self.center, radius, height)
+        reference = np.load(os.path.join(self.path, 'ref_phase_disc.npy'))
+        np.testing.assert_equal(phase, reference, 'Unexpected behavior in phase_mag_disc()')
 
     def test_phase_mag_sphere(self):
-        pass
+        radius = self.dim[2]/4
+        phase = an.phase_mag_sphere(self.dim, self.res, self.phi, self.center, radius)
+        reference = np.load(os.path.join(self.path, 'ref_phase_sphere.npy'))
+        np.testing.assert_equal(phase, reference, 'Unexpected behavior in phase_mag_sphere()')
 
     def test_phase_mag_vortex(self):
-        pass
+        radius = self.dim[2]/4
+        height = self.dim[2]/2
+        phase = an.phase_mag_vortex(self.dim, self.res, self.center, radius, height)
+        reference = np.load(os.path.join(self.path, 'ref_phase_vort.npy'))
+        np.testing.assert_equal(phase, reference, 'Unexpected behavior in phase_mag_vortex()')
 
 if __name__ == '__main__':
     suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseAnalytic)

tests/test_compliance.py

@@ -2,10 +2,11 @@
 """Testcase for the magcreator module."""
 
 
-import datetime
-import sys
 import os
+import sys
+import datetime
 import unittest
+
 import pep8
 
 
@@ -13,32 +14,32 @@ class TestCaseCompliance(unittest.TestCase):
     """Class for checking compliance of pyramid."""  # TODO: Docstring
 
     def setUp(self):
-        pass
+        self.path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0]  # Pyramid dir
 
     def tearDown(self):
-        pass
+        self.path = None
 
     def get_files_to_check(self, rootdir):
         filepaths = []
         for root, dirs, files in os.walk(rootdir):
             for filename in files:
                 if ((filename.endswith('.py') or filename.endswith('.pyx'))
-                and root != os.path.join('scripts', 'gui')):
+                and root != os.path.join(self.path, 'scripts', 'gui')):
                     filepaths.append(os.path.join(root, filename))
         return filepaths
 
     def test_pep8(self):
         # TODO: Docstring
-        files = self.get_files_to_check('pyramid') \
-            + self.get_files_to_check('scripts') \
-            + self.get_files_to_check('tests')
+        files = self.get_files_to_check(os.path.join(self.path, 'pyramid')) \
+            + self.get_files_to_check(os.path.join(self.path, 'scripts')) \
+            + self.get_files_to_check(os.path.join(self.path, 'tests'))
         ignores = ('E226', 'E128')
         pep8.MAX_LINE_LENGTH = 99
         pep8style = pep8.StyleGuide(quiet=False)
         pep8style.options.ignore = ignores
 
         stdout_buffer = sys.stdout
-        with open(os.path.join('output', 'pep8_log.txt'), 'w') as sys.stdout:
+        with open(os.path.join(self.path, 'output', 'pep8_log.txt'), 'w') as sys.stdout:
             print '<<< PEP8 LOGFILE >>>'
             print 'RUN:', datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
             print 'IGNORED RULES:', ', '.join(ignores)

tests/test_holoimage.py

@@ -2,29 +2,42 @@
 """Testcase for the holoimage module."""
 
 
+import os
 import unittest
+import numpy as np
+from numpy import pi
+
+from pyramid.phasemap import PhaseMap
 import pyramid.holoimage as hi
 
 
 class TestCaseHoloImage(unittest.TestCase):
 
     def setUp(self):
-        pass
+        self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_holoimage')
+        phase = np.zeros((4, 4))
+        phase[1:-1, 1:-1] = pi/4
+        self.phase_map = PhaseMap(10.0, phase)
 
     def tearDown(self):
-        pass
+        self.path = None
+        self.phase_map = None
 
     def test_holo_image(self):
-        pass
-
-    def test_make_color_wheel(self):
-        pass
-
-    def test_display(self):
-        pass
+        img = hi.holo_image(self.phase_map)
+        arr = np.array(img.getdata(), np.uint8).reshape(img.size[1], img.size[0], 3)
+        holo_img_r, holo_img_g, holo_img_b = arr[..., 0], arr[..., 1], arr[..., 2]
+        ref_holo_img_r = np.loadtxt(os.path.join(self.path, 'ref_holo_img_r.txt'))
+        ref_holo_img_g = np.loadtxt(os.path.join(self.path, 'ref_holo_img_g.txt'))
+        ref_holo_img_b = np.loadtxt(os.path.join(self.path, 'ref_holo_img_b.txt'))
+        hi.display(img)
+        np.testing.assert_equal(holo_img_r, ref_holo_img_r,
+                                'Unexpected behavior in holo_image() (r-component)!')
+        np.testing.assert_equal(holo_img_g, ref_holo_img_g,
+                                'Unexpected behavior in holo_image() (g-component)!')
+        np.testing.assert_equal(holo_img_b, ref_holo_img_b,
+                                'Unexpected behavior in holo_image() (b-component)!')
 
-    def test_display_combined(self):
-        pass
 
 if __name__ == '__main__':
     suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseHoloImage)