From bad56c3a1805c5403d41ab69892d6cdb0b1870e6 Mon Sep 17 00:00:00 2001
From: pc494 <pc494@cam.ac.uk>
Date: Wed, 28 Feb 2018 11:42:15 +0000
Subject: [PATCH 1/2] First commit of sub-pixel finder

---
 pyxem/utils/expt_utils.py | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/pyxem/utils/expt_utils.py b/pyxem/utils/expt_utils.py
index 31f3a890f..1e3575f06 100644
--- a/pyxem/utils/expt_utils.py
+++ b/pyxem/utils/expt_utils.py
@@ -102,6 +102,16 @@ def _polar2cart(r, theta):
     x = r * np.sin(theta)
     return x, y
 
+def _capped_gaussian(x,prefactor,sigma,center):
+    """
+    Worker functions for : subpixel_beam_finder()
+    Creates a 1D Gaussian model, and then truncates all values
+    greate than 1 to unity. There is no discretisation here.
+    """
+    model_value = prefactor*np.exp(-(x-center)**2/(2*sigma**2))
+    model_value[np.greater(model_value,np.ones_like(model_value))] = 1 #finite dynamic range
+    return np.ravel(model_value) #convert this to 1D
+
 def radial_average(z, center,cython=True):
     """Calculate the radial profile by azimuthal averaging about a specified
     center.
@@ -459,6 +469,34 @@ def refine_beam_position(z, start, radius):
 
     return c
 
+def subpixel_beam_finder(single_pattern,half_shape):
+    """
+    This routine is designed to find, to sub-pixel accuracy, the 
+    center of a pattern that has saturated a detector. The noise
+    model is Gaussian. The input dp should be approximately (within about 3 
+    pixels) centered. Use one of the other centering methods to achieve this.
+    
+    single_pattern : numpy array : A dp (for .map purposes)
+    half_shape     : int         : An int for the approx center of the pattern
+    """
+    
+    # Set up 
+    hs = half_shape #readability
+    size = np.int(half_shape/7) #this prevents fitting to the far out noise
+    pattern = single_pattern[hs-size:hs+size,hs-size:hs+size]
+    
+    if np.max(pattern) > 1:
+        raise ValueError('Patterns should be normalised to max intensity')
+    
+    #fitting
+    from scipy.optimize import curve_fit as cf
+    i_array = np.array([np.arange(hs-size,hs+size),np.arange(hs-size,hs+size),np.arange(hs-size,hs+size)])
+    x_center = cf(_capped_gaussian,i_array,np.ravel(pattern[size:size+3,:].T),p0=[2,5,hs])[0][2]
+    y_center = cf(_capped_gaussian,i_array,np.ravel(pattern[:,size:size+3].T),p0=[2,5,hs])[0][2]
+    
+    return (x_center,y_center)
+
+
 def enhance_gauss_sauvola(z, sigma_blur, sigma_enhance, k, window_size, threshold, morph_opening=True):
     z = z.astype(np.float64)
     im1 = ndi.gaussian_filter(z, sigma=sigma_blur, mode='mirror')

From 8854090b3c359f08148e717507f74736e331e6b6 Mon Sep 17 00:00:00 2001
From: pc494 <pc494@cam.ac.uk>
Date: Wed, 28 Feb 2018 12:06:19 +0000
Subject: [PATCH 2/2] Changing from tuple return to np.array return

---
 pyxem/utils/expt_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyxem/utils/expt_utils.py b/pyxem/utils/expt_utils.py
index 1e3575f06..0a5e32c52 100644
--- a/pyxem/utils/expt_utils.py
+++ b/pyxem/utils/expt_utils.py
@@ -494,7 +494,7 @@ def subpixel_beam_finder(single_pattern,half_shape):
     x_center = cf(_capped_gaussian,i_array,np.ravel(pattern[size:size+3,:].T),p0=[2,5,hs])[0][2]
     y_center = cf(_capped_gaussian,i_array,np.ravel(pattern[:,size:size+3].T),p0=[2,5,hs])[0][2]
     
-    return (x_center,y_center)
+    return np.asarray([x_center,y_center])
 
 
 def enhance_gauss_sauvola(z, sigma_blur, sigma_enhance, k, window_size, threshold, morph_opening=True):