Changed Plotting behavior, fixed analyticel solution (no more offset)

Holoimaging and phasemap now have separate display functions
Display_combine implemented in main.py
Changed sqrt(x^2 + y^2) to np.hypot(x,y)
magcreator now uses <= instead of <
pdb starts on error (implemented in main.py)

Python/magneticimaging/analytic.py

@@ -28,27 +28,32 @@ def plot_phase(phase, res, name):
 def phasemap_slab(dim, res, beta, center, width, b_0):
     '''INPUT VARIABLES'''
     y_dim, x_dim = dim
-    y0, x0 = res * center[0], res * center[1]
-    Ly, Lx = res *  width[0], res *  width[1]
+    # y0, x0 have to be in the center of a pixel, hence: cellindex + 0.5
+    y0 = res * (center[0] + 0.5)
+    x0 = res * (center[1] + 0.5)
+    # Ly, Lx have to be odd, because the slab borders should not lie in the
+    # center of a pixel (so L/2 can't be an integer)
+    Ly = res * ( width[0] + (width[0]+1)%2) 
+    Lx = res * ( width[1] + (width[1]+1)%2)   
     
     '''COMPUTATION MAGNETIC PHASE SHIFT (REAL SPACE) SLAB'''
-    
+
     coeff = b_0 * res / ( 4 * PHI_0 )
      
     def F0(x,y):
         a = np.log(x**2 + y**2)
         b = np.arctan(x / y)
         return x*a - 2*x + 2*y*b   
-      
+    
     def phiMag(x,y):
-        return coeff * ( - np.cos(beta) * ( F0(x-x0-Lx/2,y-y0-Ly/2)
-                                           -F0(x-x0+Lx/2,y-y0-Ly/2)
-                                           -F0(x-x0-Lx/2,y-y0+Ly/2)
-                                           +F0(x-x0+Lx/2,y-y0+Ly/2) )
-                         + np.sin(beta) * ( F0(y-y0-Ly/2,x-x0-Lx/2)
-                                           -F0(y-y0+Ly/2,x-x0-Lx/2)
-                                           -F0(y-y0-Ly/2,x-x0+Lx/2)
-                                           +F0(y-y0+Ly/2,x-x0+Lx/2) ) )
+        return coeff * ( - np.cos(beta) * ( F0(x-x0-Lx/2, y-y0-Ly/2)
+                                           -F0(x-x0+Lx/2, y-y0-Ly/2)
+                                           -F0(x-x0-Lx/2, y-y0+Ly/2)
+                                           +F0(x-x0+Lx/2, y-y0+Ly/2) )
+                         + np.sin(beta) * ( F0(y-y0-Ly/2, x-x0-Lx/2)
+                                           -F0(y-y0+Ly/2, x-x0-Lx/2)
+                                           -F0(y-y0-Ly/2, x-x0+Lx/2)
+                                           +F0(y-y0+Ly/2, x-x0+Lx/2) ) )
              
     '''CREATE COORDINATE GRIDS'''
     x = np.linspace(res/2,x_dim*res-res/2,num=x_dim)
@@ -69,9 +74,9 @@ def phasemap_disc(dim, res, beta, center, radius, b_0):
     coeff = - pi * res * b_0 / ( 2 * PHI_0 )
       
     def phiMag(x,y):
-        r = np.sqrt((x-x0) ** 2 + (y-y0) ** 2)      
+        r = np.hypot(x-x0, y-y0)      
         result = coeff * ((y-y0) * np.cos(beta) - (x-x0) * np.sin(beta))
-        in_or_out = 1 * (r < R) + (R / r) ** 2 * (r > R)    
+        in_or_out = 1 * (r <= R) + (R / r) ** 2 * (r > R)    
         result *= in_or_out
         return result
     

Python/magneticimaging/holoimage.py

@@ -30,15 +30,14 @@ CDICT = {'red':   [(0.00, 1.0, 0.0),
 HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
 
 
-def holo_image(phase, res, density, title):
-    '''Display a holography image with color-encoded gradient direction.
+def holo_image(phase, res, density):
+    '''Returns a holography image with color-encoded gradient direction.
     Arguments:
         phase   - the phasemap that should be displayed
         res     - the resolution of the phasemap
         density - the factor for determining the number of contour lines
-        title   - the title of the plot
     Returns:
-        None
+        Image
         
     '''
     img_holo = (1 + np.cos(density * phase * pi/2)) /2
@@ -47,8 +46,11 @@ def holo_image(phase, res, density, title):
     
     phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
     
-    phase_magnitude = np.sqrt(phase_grad_x ** 2 + phase_grad_y ** 2)    
-    phase_magnitude /= np.amax(phase_magnitude)
+    # TODO: Delete
+#    import pdb; pdb.set_trace()      
+    
+    phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)    
+    phase_magnitude /= phase_magnitude.max()
     phase_magnitude = np.sin(phase_magnitude * pi / 2)
     
     cmap = HOLO_CMAP
@@ -59,15 +61,13 @@ def holo_image(phase, res, density, title):
     green *= 255.999 * img_holo * phase_magnitude
     blue  *= 255.999 * img_holo * phase_magnitude
     rgb = np.dstack((red, green, blue)).astype(np.uint8)
+    # TODO: Which one?
+    rgb = (255.999 * img_holo.T * phase_magnitude.T 
+           * rgba_img[:, :, :3].T).T.astype(np.uint8)
     img = Image.fromarray(rgb)
     
-    fig = plt.figure()
-    ax = fig.add_subplot(111, aspect='equal')
-    ax.imshow(img)
-    ax.set_title(title + ' - Holography Image')
-    ax.set_xlabel('x-axis')
-    ax.set_ylabel('y-axis')
-
+    return img
+    
 
 def make_color_wheel():
     '''Display a color wheel for the gradient direction.
@@ -93,7 +93,9 @@ def make_color_wheel():
     red   *= 255.999 * color_wheel_magnitude
     green *= 255.999 * color_wheel_magnitude
     blue  *= 255.999 * color_wheel_magnitude
-    rgb = np.dstack((red, green, blue)).astype(np.uint8)
+    #rgb = np.dstack((red, green, blue)).astype(np.uint8)
+    # TODO Evtl. einfacher:
+    rgb = (255.999 * color_wheel_magnitude.T * rgba_img[:, :, :3].T).T.astype(np.uint8)
     color_wheel = Image.fromarray(rgb)
     
     fig = plt.figure()
@@ -101,4 +103,17 @@ def make_color_wheel():
     ax.imshow(color_wheel)
     ax.set_title('Color Wheel')
     ax.set_xlabel('x-axis')
-    ax.set_ylabel('y-axis')
\ No newline at end of file
+    ax.set_ylabel('y-axis')
+    
+
+def display_holo(holo, title, axis=None):
+    # TODO: Docstring
+    if axis == None:    
+        fig = plt.figure()
+        axis = fig.add_subplot(1,1,1, aspect='equal')
+        
+    axis.imshow(holo)
+    
+    axis.set_title(title)
+    axis.set_xlabel('x-axis')
+    axis.set_ylabel('y-axis')
\ No newline at end of file

Python/magneticimaging/magcreator.py

@@ -17,8 +17,8 @@ def slab(dim, params):
         
     '''
     center, width = params
-    shape_mag = np.array([[abs(x - center[1]) < width[1] / 2
-                       and abs(y - center[0]) < width[0] / 2
+    shape_mag = np.array([[abs(x - center[1]) <= width[1] / 2
+                       and abs(y - center[0]) <= width[0] / 2
                        for x in range(dim[1])] for y in range(dim[0])])
     return shape_mag
 
@@ -35,7 +35,7 @@ def disc(dim, params):
     '''
     center, radius = params
     shape_mag = np.array([[(np.sqrt((x - center[1]) ** 2 +
-                                    (y - center[0]) ** 2) < radius)
+                                    (y - center[0]) ** 2) <= radius)
                                     for x in range(dim[1])] 
                                     for y in range(dim[0])])    
     return shape_mag
@@ -76,13 +76,11 @@ def alternating_filaments(dim, params):
     shape_mag = np.zeros(dim)
     if x_or_y == 'y':
         # TODO: List comprehension:
-        for i in range(dim[1]):
-            if i % spacing == 0:
-                shape_mag[:, i] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1
+        for i in range(0, dim[1], spacing):
+            shape_mag[:, i] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1
     else:
-        for i in range(dim[0]):
-            if i % spacing == 0:
-                shape_mag[i, :] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1
+        for i in range(0, dim[0], spacing):
+            shape_mag[i, :] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1            
     return shape_mag
 
     

Python/magneticimaging/phasemap.py

@@ -68,7 +68,8 @@ def real_space(mag_data, b_0=1, v_0=0, v_acc=30000):
     '''
     # TODO: Expand docstring!
 
-#    import pdb; pdb.set_trace() # TODO   
+    # TODO: Delete
+#    import pdb; pdb.set_trace()    
     
     res  = mag_data.res
     x_dim, y_dim, z_dim = mag_data.dim
@@ -102,21 +103,53 @@ def real_space(mag_data, b_0=1, v_0=0, v_acc=30000):
     yF = np.linspace(-(y_dim-1), y_dim-1, num=2*y_dim-1)
     xxF, yyF = np.meshgrid(xF,yF)
     
-    F_cos_part = F_part(xxF, yyF)
-    F_sin_part = F_part(yyF, xxF)
-    display(F_cos_part, res, 'F_cos_part')
-    display(F_sin_part, res, 'F_sin_part')
+    F_part_cos = F_part(xxF, yyF)
+    F_part_sin = F_part(yyF, xxF)
+    
+    display_phase(F_part_cos, res, 'F_part_cos')
+    display_phase(F_part_sin, res, 'F_part_sin')      
     
     phase = np.zeros((y_dim,x_dim))
     
+    for j in y:
+        for i in x:
+            phase += phiMag(xx, yy, i, j, coeff[j,i], beta[j,i])
+            
+    return phase
+    
+    
+    xF = np.linspace(-(x_dim-1), x_dim-1, num=2*x_dim-1)
+    yF = np.linspace(-(y_dim-1), y_dim-1, num=2*y_dim-1)
+    xxF, yyF = np.meshgrid(xF,yF)
+    
+    F_part_cos = F_part(xxF, yyF)
+    F_part_sin = F_part(yyF, xxF)
+    
+    
+    display_phase(F_part_cos, res, 'F_part_cos')
+    display_phase(F_part_sin, res, 'F_part_sin')    
+    
+    def phiMag2(xx, yy, i, j):
+        #x_ind = xxF[yy]
+        
+        return coeff[j,i] * ( - np.cos(beta[j,i]) * F_part_cos[yy.min()-j:(2*yy.max()-1)-j, xx.min()-i:(2*xx.max()-1)-i]
+                              + np.sin(beta[j,i]) * F_part_sin[xx.min()-i:(2*xx.max()-1)-i, yy.min()-j:(2*yy.max()-1)-j] )
+    
+    
+    phase2 = np.zeros((y_dim*x_dim, y_dim, x_dim))
+    
     for j in range(y_dim):
         for i in range(x_dim):
-            phase += phiMag(xx, yy, x[i], y[j], coeff[j][i], beta[j][i])
-  
-    return phase
-	
+            phase2 += phiMag2(xx, yy, 0, 0)
+                
+   
+    
+    phase2 = np.sum(phase2, axis=0)
+    
+    return phase2
+    
 	
-def display(phase, res, title):
+def display_phase(phase, res, title, axis=None):
     '''Display the phasemap as a colormesh.
     Arguments:
         phase - the phasemap that should be displayed
@@ -125,20 +158,25 @@ def display(phase, res, title):
     Returns:
         None
         
-    '''    
-    fig = plt.figure()
-    ax = fig.add_subplot(111, aspect='equal')
+    '''
+    if axis == None:
+        fig = plt.figure()
+        axis = fig.add_subplot(1,1,1, aspect='equal')
     
-    plt.pcolormesh(phase, cmap='gray')
+    im = plt.pcolormesh(phase, cmap='gray')
 
-    ticks = ax.get_xticks()*res
-    ax.set_xticklabels(ticks.astype(int))
-    ticks = ax.get_yticks()*res
-    ax.set_yticklabels(ticks.astype(int))
+    ticks = axis.get_xticks()*res
+    axis.set_xticklabels(ticks.astype(int))
+    ticks = axis.get_yticks()*res
+    axis.set_yticklabels(ticks.astype(int))
 
-    ax.set_title(title)
-    ax.set_xlabel('x-axis [nm]')
-    ax.set_ylabel('y-axis [nm]')
+    axis.set_title(title)
+    axis.set_xlabel('x-axis [nm]')
+    axis.set_ylabel('y-axis [nm]')
+    
+    fig = plt.gcf()
+    fig.subplots_adjust(right=0.85)
+    cbar_ax = fig.add_axes([0.9, 0.15, 0.02, 0.7])
+    fig.colorbar(im, cax=cbar_ax)
     
-    plt.colorbar()
     plt.show()
\ No newline at end of file

Python/main.py

@@ -5,13 +5,14 @@ Created on Wed Apr 03 11:15:38 2013
 @author: Jan
 """
 
-
+import matplotlib.pyplot as plt
 import magneticimaging.magcreator as mc
 import magneticimaging.dataloader as dl
 import magneticimaging.phasemap as pm
 import magneticimaging.holoimage as hi
 import magneticimaging.analytic as an
 import time
+import pdb, traceback, sys
 from numpy import pi
     
 
@@ -33,21 +34,21 @@ def phase_from_mag():
     padding = 20
     density = 10
     
-    dim = (20, 20)  # in px (y,x)
+    dim = (50, 50)  # in px (y,x)
     res = 10.0  # in nm
     beta = pi/4
     
-    plot_mag_distr = False
+    plot_mag_distr = True
     
     # Slab:
     shape_fun = mc.slab
-    center = (10, 10)  # in px (y,x)
-    width = (10, 10)  # in px (y,x)
+    center = (24, 24)  # in px (y,x) index starts with 0!
+    width  = (25, 25)  # in px (y,x)
     params = (center, width)
 #    # Disc:
 #    shape_fun = mc.disc
-#    center = (30, 30)  # in px (y,x)
-#    radius = 15
+#    center = (4, 4)  # in px (y,x)
+#    radius = 2
 #    params = (center, radius)
 #    # Filament:
 #    shape_fun = mc.filament
@@ -80,39 +81,58 @@ def phase_from_mag():
     phase_fft = pm.fourier_space(mag_data, b_0, padding)   
     tocf = time.clock()
     print 'Time for Fourier Space Approach: ' + str(tocf - ticf)
-    pm.display(phase_fft, mag_data.res, 'Fourier Space Approach')
-    hi.holo_image(phase_fft, mag_data.res, density, 'Fourier Space Approach')
+    holo_fft = hi.holo_image(phase_fft, mag_data.res, density)
+    display_combined(phase_fft, mag_data.res, holo_fft, 
+                     'Fourier Space Approach')
     # numerical solution Real Space:
     ticr = time.clock()
     phase_real = pm.real_space(mag_data, b_0)
     tocr = time.clock()
     print 'Time for Real Space Approach:    ' + str(tocr - ticr)
     print 'Fourier Approach is ' + str((tocr-ticr) / (tocf-ticf)) + ' faster!'
-    pm.display(phase_real, mag_data.res, 'Real Space Approach')
-    hi.holo_image(phase_real, mag_data.res, density, 'Real Space Approach')
+    holo_real = hi.holo_image(phase_real, mag_data.res, density)
+    display_combined(phase_real, mag_data.res, holo_real, 
+                     'Real Space Approach')
     
     '''ANALYTIC SOLUTION'''
     # analytic solution slab:
     phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
-    pm.display(phase, res, 'Analytic Solution - Slab')
-    hi.holo_image(phase, res, density, 'Analytic Solution - Slab')
+    holo  = hi.holo_image(phase, res, density)
+    display_combined(phase, res, holo, 'Analytic Solution - Slab')
 #    # analytic solution disc:
 #    phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
-#    pm.display(phase, res, 'Analytic Solution - Disc')
-#    hi.holo_image(phase, res, density, 'Analytic Solution - Disc')
+#    holo  = hi.holo_image(phase, res, density)
+#    display_combined(phase, res, holo, 'Analytic Solution - Disc')
 #    # analytic solution sphere:
 #    phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
-#    pm.display(phase, res, 'Analytic Solution - Sphere')
-#    hi.holo_image(phase, res, density, 'Analytic Solution - Sphere')
+#    holo  = hi.holo_image(phase, res, density)
+#    display_combined(phase, res, holo, 'Analytic Solution - Sphere')
+
     
     '''DIFFERENCES'''
     diff_real_to_ana = phase_real - phase
     diff_fft_to_ana  = phase_fft  - phase 
-    diff_real_to_fft = phase_real - phase_fft
-    pm.display(diff_real_to_ana, res, 'Difference: Analytic - Real')
-    pm.display(diff_fft_to_ana,  res, 'Difference: Analytic - Fourier')
-    pm.display(diff_real_to_fft, res, 'Difference: Real - Fourier') 
+    diff_real_to_fft = phase_fft - phase_real
+    pm.display_phase(diff_real_to_ana, res, 'Difference: Analytic - Real')
+    pm.display_phase(diff_fft_to_ana,  res, 'Difference: Analytic - Fourier')
+    pm.display_phase(diff_real_to_fft, res, 'Difference: Fourier - Real') 
+ 
+    
+def display_combined(phase, res, holo, title):
+    fig = plt.figure(figsize=(14, 7))
+    fig.suptitle(title, fontsize=20)
+    
+    holo_axis = fig.add_subplot(1,2,1, aspect='equal')
+    hi.display_holo(holo, 'Holography Image', holo_axis)
+    
+    phase_axis = fig.add_subplot(1,2,2, aspect='equal')
+    pm.display_phase(phase, res, 'Phasemap', phase_axis)
     
     
 if __name__ == "__main__":
-    phase_from_mag()
\ No newline at end of file
+    try:
+        phase_from_mag()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)
\ No newline at end of file