update code

docs/examples/droplet_fusion/droplet_fusion.rst

@@ -93,18 +93,32 @@ Fit the data and plot the result::
     force = force_selection.data
 
     popt, pcov = curve_fit(relaxation_model, time, force, [0.1, force[0], 0, 0])
-    print(popt)
     plt.figure()
-    plt.plot(time,force)
-    plt.plot(time,relaxation_model(time,*popt),label=f'$\\tau$ = {popt[0]:0.2f}s')
+    plt.plot(time, force)
+    plt.plot(time, relaxation_model(time,*popt), label=fr"$\tau$"  = {popt[0]:0.2f}s')
     plt.ylabel(r"x-coordinate ($\mu$m)")
     plt.xlabel("Time (s)")
     plt.legend()
     plt.show()
 
 .. image:: fit.png
 
-The relaxation time obtained from the fit is 0.05 seconds. 
+The array :math:`popt` contains all the fitted parameters::
+
+    >>>> print(popt)
+
+    [ 0.04754597  5.39347958 -4.5188036  -6.2716303 ]
+
+The first parameter in :math:`popt` is :math:`𝜏` and the other 3 parameters are :math:`a`, :math:`b` and :math:`c` respectively, as defined in the model above.
+The matrix :math:`pcov` is the covariance matrix and the standard deviation errors in the fitted parameters can be obtained as 
+
+    >>>> np.sqrt(np.diag(pcov))
+
+    [0.00021354, 0.01104559, 0.0218276 , 0.00779505]
+
+See the :hyperlink:`documentation on curve_fit <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html>` for more details on the fitting procedure.
+
+The relaxation time obtained from the fit is 0.0475 +- 0.0002 seconds. 
 
 Now, we will proceed to determine the size of the droplets before the fusion event.
 
@@ -172,40 +186,57 @@ The threshold may need to be optimized for your data::
 
 .. image:: image_segmentation.png
 
-Estimate the size of the droplet by measuring the length of the longest axis and the shortest axis and taking the average::
+For this scan, the fast axis is along the horizontal coordinate (you can check the direction of the fast axis by typing :attr:`scan.fast_axis <lumicks.pylake.scan.Scan.fast_axis>`). 
+Therefore, we estimate the size of the droplets by looking at the width of the identified objects:::
+
+    def get_center_and_width(scan, mask, axis):
+        """Grabs the center and width along the fast scanning axis"""
+        widths = np.sum(mask, axis=axis)
+        max_width = np.max(widths)
+
+        # Grab the position
+        coordinate_weighted_mask = np.indices(mask.shape)[axis] * mask
+        centers = np.sum(coordinate_weighted_mask, axis=axis) / np.clip(np.sum(mask, axis=axis), 1, np.inf)
+        
+        # Since some scanlines can have the same width, we'd want the vertical position to be the average of these
+        max_scanline = int(np.mean(np.nonzero(max_width == widths)[0]))
+
+        if axis:
+            center = (centers[max_scanline], max_scanline)
+        else:
+            center = (max_scanline, centers[max_scanline])
+        
+        return center, max_width
+
+
+    def plot_width(scan, center, width, axis):
+        plt.plot(center[0], center[1], "ko")
+        if axis:
+            plt.plot([center[0] - 0.5 * width, center[0] + 0.5 * width], [center[1], center[1]])
+        else:
+            plt.plot([center[0], center[0]], [center[1] - 0.5 * width, center[1] + 0.5 * width])
+
+    average_droplet_radii = np.array([]) 
 
-    regions = regionprops(label_img)
-    
     fig, ax = plt.subplots()
     ax.imshow(image, cmap=plt.cm.gray)
-    average_droplet_radii = np.array([])  # List of arrays to collect all droplet radii
-
-    for n, props in enumerate(regions):
-        y0, x0 = props.centroid
-        orientation = props.orientation
-        x1 = x0 + math.cos(orientation) * 0.5 * props.axis_minor_length
-        y1 = y0 - math.sin(orientation) * 0.5 * props.axis_minor_length
-        x2 = x0 - math.sin(orientation) * 0.5 * props.axis_major_length
-        y2 = y0 - math.cos(orientation) * 0.5 * props.axis_major_length
-
-        ax.plot((x0, x1), (y0, y1), "-r", linewidth=2.5)
-        ax.plot((x0, x2), (y0, y2), "-r", linewidth=2.5)
-        ax.plot(x0, y0, ".g", markersize=15)
-
-        minr, minc, maxr, maxc = props.bbox
-        bx = (minc, maxc, maxc, minc, minc)
-        by = (minr, minr, maxr, maxr, minr)
-        ax.plot(bx, by, "-b", linewidth=2.5)
-       
-        average_droplet_radii = np.append(
-            average_droplet_radii, 0.25 * (props.minor_axis_length + props.major_axis_length)
-        )
-    plt.title("The short and long axis of the identified objects")
+    ax.imshow(image, cmap=plt.cm.gray)
+    axis = 1 if scan.fast_axis == "X" else 0
+    center, width = get_center_and_width(scan, label_img == 1, axis)
+    average_droplet_radii = np.append(
+            average_droplet_radii, 0.5* width)
+    plot_width(scan, center, width, axis)
+    center, width = get_center_and_width(scan, label_img == 2, axis)
+    average_droplet_radii = np.append(
+            average_droplet_radii, 0.5* width)
+    plot_width(scan, center, width, axis)
+
+    plt.title("Width along the fast axis")
     plt.show()
 
 .. image:: droplet_size_pixels.png
 
-The array `average_droplet_radii` contains the average of the large and small radius for all the droplets in the image, in pixels. 
+The array `average_droplet_radii` contains the radii of both droplets in the image, in pixels. 
 Below we are multiplying this array by the pixel size in micron to obtain the radii in micron::
 
     average_droplet_radii_um = average_droplet_radii * scan.pixelsize_um[0]
@@ -214,7 +245,7 @@ The average radii for the droplets in micrometers are::
 
     >>> print(average_droplet_radii)
 
-    [11.05117733 10.73977514]
+    [11.5 11]
 
 We now determined the relaxation time as well as the droplet radii. The next step is to measure these two quantities for many different fusion events, plot 𝜏 vs average radius and determine the slope.