New background correction option and improved junge calculation

docs/notebooks/background_correction.ipynb

@@ -192,6 +192,72 @@
     "np.abs(imbg_before - processing_pipeline.data['imbg']).sum()\n"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Correct images by subtracting vs dividing average background\n",
+    "The CorrectBackgroundAccurate class have two modes for a subtracting background correction (divide_bg=False), \n",
+    "or dividing background correction (divide_bg=True) that provides the corrected image ('im_corrected') for further analysis.\n",
+    "For dividing background mode, the zero-value pixels of the average background image are initially rescaled to 1/255 to prevent division by zero.\n",
+    "For more information, refer to: https://doi.org/10.1016/j.marpolbul.2016.11.063).\n",
+    "You can select the subtracting/dividing correction modes used in the pipeline to process raw images, as illustrated below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Corrected image uisng subtracting method\n",
+    "im_corrected_subtract = processing_pipeline.data['im_corrected'].copy()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Update pipeline config and background step with dividing background correction (divide_bg=True)\n",
+    "\n",
+    "pipeline_config['steps']['correctbackground']['divide_bg'] = True\n",
+    "\n",
+    "# Run the first N images to creat the background\n",
+    "for filename in image_files[:NUM_IMAGES_FOR_BACKGROUND]:\n",
+    "    processing_pipeline.run(filename)\n",
+    "\n",
+    "# Now process one of the raw images\n",
+    "processing_pipeline.run(image_files[NUM_IMAGES_FOR_BACKGROUND])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Corrected image uisng dividing mode\n",
+    "im_corrected_division = processing_pipeline.data['im_corrected'].copy()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot raw, im_corrected by subtracing and dividing averaged background image\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(3*6, 5))\n",
+    "axes[0].imshow(processing_pipeline.data['imraw'])\n",
+    "axes[0].set_title(f'Raw image #{NUM_IMAGES_FOR_BACKGROUND}')\n",
+    "axes[1].imshow(im_corrected_subtract)\n",
+    "axes[1].set_title('Corrected image by background subtraction')\n",
+    "axes[2].imshow(im_corrected_division)\n",
+    "axes[2].set_title('Corrected image by background division')"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

pyopia/__init__.py

@@ -1 +1 @@
-__version__ = '2.6.2'
+__version__ = '2.7.0'

pyopia/background.py

@@ -95,9 +95,11 @@ def shift_bgstack_fast(bgstack, imbg, imnew):
     return bgstack, imbg
 
 
-def correct_im_accurate(imbg, imraw):
+def correct_im_accurate(imbg, imraw, divide_bg=False):
     '''
-    Corrects raw image by subtracting the background and scaling the output
+    Corrects raw image by subtracting or dividing the background and scaling the output
+
+    For dividing method see: https://doi.org/10.1016/j.marpolbul.2016.11.063)
 
     There is a small chance of clipping of imc in both crushed blacks and blown
     highlights if the background or raw images are very poorly obtained
@@ -108,17 +110,26 @@ def correct_im_accurate(imbg, imraw):
         background averaged image
     imraw : float64
         raw image
+    divide_bg : (bool, optional)
+        If True, the correction will be performed by dividing the raw image by the background
+        Default to False
 
     Returns
     -------
     im_corrected : float64
         corrected image, same type as input
     '''
 
-    im_corrected = imraw - imbg
-    im_corrected += (1 / 2 - np.percentile(im_corrected, 50))
-
-    im_corrected += 1 - im_corrected.max()
+    if divide_bg:
+        imbg = np.clip(imbg, a_min=1/255, a_max=None)   # Clipping the zero_value pixels
+        im_corrected = imraw / imbg
+        im_corrected += (1 / 2 - np.percentile(im_corrected, 50))
+        im_corrected -= im_corrected.min()   # Shift the negative values to zero
+        im_corrected = np.clip(im_corrected, a_min=0, a_max=1)
+    else:
+        im_corrected = imraw - imbg
+        im_corrected += (1 / 2 - np.percentile(im_corrected, 50))
+        im_corrected += 1 - im_corrected.max()   # Shift the positive values exceeding unity to one
 
     return im_corrected
 
@@ -207,7 +218,7 @@ class CorrectBackgroundAccurate():
     ----------
     bgshift_function : (string, optional)
         Function used to shift the background. Defaults to passing (i.e. static background)
-        Available options are 'accurate', 'fast', or 'pass' to apply a statick background correction:
+        Available options are 'accurate', 'fast', or 'pass' to apply a static background correction:
 
         :func:`pyopia.background.shift_bgstack_accurate`
 
@@ -220,6 +231,10 @@ class CorrectBackgroundAccurate():
         The key in Pipeline.data of the image to be background corrected.
         Defaults to 'imraw'
 
+    divide_bg : (bool)
+        If True, it performs background correction by dividing the raw image by the background.
+        Default to False.
+
     Returns
     -------
     data : :class:`pyopia.pipeline.Data`
@@ -256,10 +271,11 @@ class CorrectBackgroundAccurate():
     Then you could use :class:`pyopia.pipeline.CorrectBackgroundNone` if you need to instead.
     '''
 
-    def __init__(self, bgshift_function='pass', average_window=1, image_source='imraw'):
+    def __init__(self, bgshift_function='pass', average_window=1, image_source='imraw', divide_bg=False):
         self.bgshift_function = bgshift_function
         self.average_window = average_window
         self.image_source = image_source
+        self.divide_bg = divide_bg
 
     def _build_background_step(self, data):
         '''Add one layer to the background stack from the raw image in data pipeline, and update the background image.'''
@@ -284,7 +300,7 @@ def __call__(self, data):
             data['skip_next_steps'] = True
             return data
 
-        data['im_corrected'] = correct_im_accurate(data['imbg'], data[self.image_source])
+        data['im_corrected'] = correct_im_accurate(data['imbg'], data[self.image_source], divide_bg=self.divide_bg)
 
         match self.bgshift_function:
             case 'pass':

pyopia/statistics.py

@@ -534,7 +534,7 @@ def get_j(dias, number_distribution):
 
     # use polyfit to obtain the slope of the ditriubtion in log-space (which is
     # assumed near-linear in most parts of the ocean)
-    p = np.polyfit(np.log(dias[ind]), np.log(number_distribution[ind]), 1)
+    p = np.polyfit(np.log(dias[ind]), np.log(number_distribution[ind], where=number_distribution[ind] > 0), 1)
     junge_slope = p[0]
     return junge_slope