Updated prediction module

sl_main.py

@@ -441,7 +441,18 @@ def livestack(self):
             if self.refresh is not None:
                 self.refresh.cancel()
             self.newimage.stop()
-
+
+    def show_suggestions(self):
+        reg = predict_regions(self.data.rawimg, self.data.imgdims)
+        self.imgdisplay.img.autoscale(False)
+        for y,x,a in reg:
+            w = np.sqrt(a)*2.5
+            xlist = [x-w,x-w,x+w,x+w,x-w]
+            ylist = [y-w,y+w,y+w,y-w,y-w]
+            ylist = [self.data.imgdims[1]-q for q in ylist]
+            self.imgdisplay.img.plot(xlist, ylist, color='green')
+        self.imgdisplay.canvas.draw()
+
     def setcommands(self):
         # For setting .hdr file
         self.browsebtn.config(command = self.hdrdialog)
@@ -458,6 +469,7 @@ def setcommands(self):
         self.showroischk.config(command = lambda: self.imgdisplay.redraw(self))
         self.autoalignchk.config(command = self.genI0IT)
         self.autosegment.config(command = self.autoselect_rois)
+        self.suggestnext.config(command = self.show_suggestions)
 
         # Image nav buttons
         self.imgfirst.config(command = lambda: self.setxim(ind = 0))

sl_proc.py

@@ -1,44 +1,17 @@
 import numpy as np
 
-#from scipy.ndimage import fourier_shift
 from scipy.ndimage.filters import gaussian_filter
-from scipy.ndimage.interpolation import shift #, zoom
+from scipy.ndimage.interpolation import shift
 from scipy.ndimage.morphology import binary_erosion
 
 from skimage.feature import register_translation
-from skimage.filters import sobel
+from skimage.filters import sobel, threshold_otsu
+from skimage.measure import regionprops
+from skimage.morphology import label
+from skimage.io import imsave
 
-def otsu(img):
-    """Calculates a threshold for a integer-
-    valued image using Otsu's method."""
-
-    tmp = np.ndarray.flatten(img.astype(np.int))
-
-    mymin = np.min(tmp)
-    mymax = np.max(tmp)
-    myavg = np.mean(tmp)
-
-    bins = np.array([n for n in range(mymin, mymax+1)])
-    vals = np.zeros_like(bins, dtype = np.float)
-    for x in tmp:
-        vals[x-mymin] += 1.0
-    vals /= np.sum(vals)
-
-    variance = np.zeros_like(bins, dtype = np.float)
-    wlow = vals[0]
-    mulow = mymin
-    whigh = 1.0 - wlow
-    muhigh = (myavg - wlow*mulow)/whigh
-    variance[0] = wlow*whigh*(mulow - muhigh)*(mulow - muhigh)
-
-    for n in range(1, mymax - mymin):
-        wlow += vals[n]
-        mulow = np.sum( bins[:n]*vals[:n] )/np.sum(vals[:n])
-        whigh = 1.0 - wlow
-        muhigh = (myavg - wlow*mulow)/whigh
-        variance[n] = wlow*whigh*(mulow - muhigh)*(mulow - muhigh)
-
-    return np.float(mymin + np.argmax(variance))
+import keras
+upsampler = keras.models.load_model('upsample_model.h5')
 
 def calculate_shift(imgA, imgB, pxwidth):
     """Uses the registration algorithm from scikit-image
@@ -49,7 +22,7 @@ def calculate_shift(imgA, imgB, pxwidth):
     filtA = gaussian_filter(imgA, min(sigma, 3.0))
     filtB = gaussian_filter(imgB, min(sigma, 3.0))
 
-    thr = otsu(filtB)
+    thr = threshold_otsu(filtB)
     edgecheck = 1*(filtB < thr)
     # Check whether the object lies along a boundary of the image,
     # and if so apply a Sobel filter to use edges to align
@@ -62,70 +35,6 @@ def calculate_shift(imgA, imgB, pxwidth):
     shift = register_translation(filtA, filtB, upsample_factor = min(1000.0, 100.0/pxwidth))
     return [-shift[0][1], shift[0][0]]
 
-# def calculate_shift_scipy(imgA, imgB, pxwidth, zfmax = None):
-#     """Given two images and the width of each pixel (in microns),
-#     applies a Gaussian filter to both images, using a Gaussian
-#     with a width of 0.15 microns (i.e., approximately 3 times
-#     the spatial resolution of a typical STXM image). Calculates
-#     the 2D cross-correlation and returns the shift that maximizes
-#     the cross-correlation, upsampling to calculate the shift to
-#     within 10 nanometers or 1/10 pixel, whichever is finer.
-#     """
-#     
-#     # Apply Gaussian filter
-#     sigma = 0.15/pxwidth
-#     filtA = gaussian_filter(imgA, sigma, mode = 'nearest')
-#     filtB = gaussian_filter(imgB, sigma, mode = 'nearest')
-#     
-#     # Calculate shift
-#     if zfmax is None:
-#         lim = max(1.0, 10.0, pxwidth/0.01)
-#     else:
-#         lim = zfmax
-# 
-#     if lim == 1.0:
-#         tmpA = filtA[::]
-#         tmpB = filtB[::]
-#     else:
-#         zf = lim/2.0
-#         tmpA = zoom(filtA, zf)
-#         tmpB = zoom(filtB, zf)
-#     
-#     image_product = np.fft.fft2(tmpA) * np.fft.fft2(tmpB).conj()
-#     corr = np.fft.fftshift(np.fft.ifft2(image_product))            
-#         
-#     s = corr.shape
-#     origin = [(q - (q%2))/2 for q in s]
-#     
-#     x = np.argmax(corr.real)
-#     center = [int(x/s[1]), x%s[1]]
-# 
-#     finalsh = [np.float(origin[0] - center[0])/zf, np.float(origin[1] - center[1])/zf]
-#     print(lim, zf, x, finalsh)
-#     
-#     if lim != 1.0:
-#         tmpA = zoom(tmpA, 2.0)
-#         tmpB = zoom(shift(tmpB, [-finalsh[0], -finalsh[1]], mode = 'nearest'), 2.0)
-#         me = np.mean(tmpA)
-#         st = np.std(tmpA)
-#         nrmA = (tmpA - me)/st
-#         
-#         corr = np.zeros((3,3), dtype = np.float)
-#         for ii in range(3):
-#             rollB = np.roll(tmpB, ii - 1, axis = 0)
-#             for jj in range(3):
-#                 nrmB = np.roll(rollB, jj - 1, axis = 1)
-#                 corr[ii][jj] += np.sum(nrmA*nrmB)
-#         
-#         x = np.argmax(corr)
-#         center = [int(x/3), x%3]
-#         
-#         finalsh[0] += np.float(1 - center[0])/lim
-#         finalsh[1] += np.float(1 - center[1])/lim
-#     
-#     
-#     return finalsh
-
 def alignoneimage(img, sh):
     """Given an image and x,y shifts in pixels, returns a
     shifted image in which pixels beyond the boundaries of
@@ -183,8 +92,8 @@ def genmap(raw, shift):
     """Given two images and the shift between then, converts
     both to OD and then takes the difference of the two images
     to generate an elemental map."""
-    thr_0 = otsu(raw[0])
-    thr_1 = otsu(raw[1])
+    thr_0 = threshold_otsu(raw[0])
+    thr_1 = threshold_otsu(raw[1])
 
     fl_0 = np.ndarray.flatten(raw[0])
     I0_0 = np.mean([x for x in fl_0 if x > thr_0])
@@ -205,7 +114,7 @@ def autoseg(dataset, bdy):
     regions with a boundary with width of the specified
     number of pixels."""
 
-    thr = otsu(dataset.rawimg)
+    thr = threshold_otsu(dataset.rawimg)
 
     i0_mask = (dataset.displayimg > thr)*(dataset.keeppx)
     if int(bdy/2) + bdy%2 > 0:
@@ -245,5 +154,37 @@ def regridlinescan(rawimg, energies, dims):
 
     return(newgrid)
 
+def predict_regions(rawimg, dims):
+    """Given an image and dimensions, predicts an upsampled version
+    and returns centroids of regions with (a) OD between 0.5 and 1.5
+    and (b) area at least 1 um^2."""
+
+    # Apply neural net upsampling
+    imgin = np.array([[rawimg]]).astype(np.float32)
+    usamp = upsampler.predict(imgin)[0,0]
+
+    # Convert to OD
+    thr = threshold_otsu(usamp)
+    myi0 = np.mean(usamp[usamp>thr])
+    usampOD = np.log(myi0/usamp)
+
+    mask = (usampOD>0.5)*(usampOD<1.5)
+    a = (4*rawimg.shape[0]-usampOD.shape[0])//2
+    b = (4*rawimg.shape[1]-usampOD.shape[1])//2
+    mask = np.pad(mask, ((a,a),(b,b)), mode='constant')
+    reg = regionprops(label(mask))
+    okcentroids = []
+
+    for r in reg:
+        scaledarea = r.area*dims[2]*dims[3]/16
+        if scaledarea>1:
+            okcentroids.append([r.centroid[0]*dims[2]/4,
+                                r.centroid[1]*dims[3]/4,
+                                scaledarea])
+
+    okcentroids.sort(key=lambda ca: -ca[2])
+
+    return okcentroids[:10]
+
 if __name__ == '__main__':
     main()