Merge pull request #24 from KitwareMedical/new-butterfly-clarius

ENH: more robust preprocessing for Clarius and Butterfly

itkpocus/itkpocus/butterfly.py

@@ -127,10 +127,17 @@ def _calc_crop(npvid):
     # center of image
     cr, cc = (np.array(npimg.shape) / 2.0).astype('int')
 
-    sr = int(0.05 * npimg.shape[0])
+    #sr = int(0.05 * npimg.shape[0])
+    sr = cr
     er = cr
-
-    while er < npimg.shape[0]-1 and len(np.argwhere(npimg[er,:])) > npimg.shape[1] * 0.25:
+    row_cutoff = npimg.shape[1] * 0.25
+
+    while sr > 0 and len(np.argwhere(npimg[sr,:])) > row_cutoff:
+        sr -= 1
+
+    sr += 5 # get rid of gray bar
+
+    while er < npimg.shape[0]-1 and len(np.argwhere(npimg[er,:])) > row_cutoff:
         er += 1
 
     er -= 15 # get rid of angling shadowbox on butterfly

itkpocus/itkpocus/clarius.py

@@ -2,6 +2,8 @@
 import skvideo
 import skvideo.io
 from skimage.transform import SimilarityTransform
+from skimage.restoration import inpaint
+from skimage.morphology import dilation, disk
 from scipy.signal import find_peaks
 from itkpocus.util import get_framerate
 import itk
@@ -10,47 +12,117 @@
 Preprocessing and device-specific IO for the Clarius HD C7.
 '''
 
-def _find_spacing_and_crop(npimg):
+def _get_rightmost_ruler_from_mask(mask):
     '''
-    Find the ruler in the image to calculate the physical dimensions and then crop to the B-mode data.
+    Finds the ruler as long as it is the rightmost thing on the mask.
     
     Parameters
     ----------
-    npimg : ndarray
-        MxNx3
-    
+    mask : ndarray
+        MxN where the ruler and other overlay elements are non-zero values
+        
     Returns
-    =======
-    spacing
-    crop : ndarray
-        [[ymin,ymax],[xmin,xmax]]
+    -------
+    spacing : float
+    ruler_col : int
+    ruler_peaks : ndarray of int
     '''
-    colsum = np.sum(npimg[:,:,0], axis=0)
-    col_peaks, col_props = find_peaks(colsum, prominence=20)
+    col_prominence=20 # strength of ruler signal column-wise
+    row_prominence=1 # strength of each ruler tick row-wise
+    row_width=5
+    consistency_threshold = 10 # max number of pixels a ruler step can differ in
+    outlier_threshold = 0.1
+    tick_dist = 10.0 # physical distance between ticks in mm
+
+    colsum = np.sum(mask, axis=0)
+    col_peaks, col_props = find_peaks(colsum, prominence=col_prominence)
     ruler_col = col_peaks[-1]
-    ruler_peaks, ruler_props = find_peaks(npimg[:,ruler_col,0])
+    ruler_peaks, ruler_props = find_peaks(mask[:,ruler_col], prominence=row_prominence, width=row_width)
+
+    diffs = np.abs(ruler_peaks[1:] - ruler_peaks[:-1])
+    median_diff = np.median(diffs)
+    outliers_perc =  len(np.argwhere(diffs - median_diff > consistency_threshold)) / len(diffs)
+    if outliers_perc > outlier_threshold:
+        raise ValueError('Max percentage of outlier ruler ticks greater than threshold (ruler likely not detected or obscured): {} > {}'.format(outliers_perc, outlier_threshold))
+
     # TODO: stress test this with max zoomed video, not sure if the ruler would switch to 0.5 cm spacing or < 4 circles would be on the video
-    spacing = 10.0 / np.max(ruler_peaks[1:] - ruler_peaks[0:len(ruler_peaks)-1]) # 10 mm per ruler spacing, we get the max peak distance because half-circles on the top and bottom shift the peaks
-
-    # find left/right bounds
-    # just going assume the ultrasound is in the middle of the image and work outwards to find the letterboxing
-    midcol = npimg.shape[1]/2.0 # not worried if non-integer
-    zerocol = np.argwhere(colsum == 0)
-    leftbound = np.max(zerocol[zerocol < midcol])+1
-    rightbound = np.min(zerocol[midcol < zerocol])-1
-
-    midrow = npimg.shape[0]/2.0
-    rowsum = np.sum(npimg[:,leftbound:rightbound+1,0], axis=1)
-    zerorow = np.argwhere(rowsum == 0)
-
-    if len(zerorow) == 0: # fills the whole image vertically
-        topbound=0
-        bottombound=npimg.shape[0]-1
-    else:
-        topbound = np.max(zerorow[zerorow < midrow])
-        bottombound = np.min(zerorow[midrow < zerorow])
-
-    return spacing, np.array([[topbound, bottombound], [leftbound, rightbound]])
+    spacing = tick_dist / median_diff
+
+    return spacing, ruler_col, ruler_peaks
+
+def _inpaint(npimg, mask):
+    '''
+    Inpaint the locations in npimg corresponding to true values in mask.
+    
+    Parameters
+    ----------
+    npimg : ndarray
+        MxN
+    mask : ndarray
+        MxN
+        
+    Returns
+    -------
+    ndarray
+    '''
+    return inpaint.inpaint_biharmonic(npimg, mask)
+
+def _find_crop(npimg_clean):
+    '''
+    Parameters
+    ----------
+    npimg_clean : ndarray
+        MxN with overlay elements removed/inpainted, intensities between 0 and 1.0
+    
+    Returns
+    -------
+    ndarray
+        [[toprow, bottomrow], [leftcol, rightcol]]
+    '''
+
+    bg_cutoff = 10.0 / 255.0
+    col_cutoff_perc = 0.1
+    row_cutoff_perc = 0.1
+
+    non_bg = npimg_clean > bg_cutoff
+
+    colsum = np.sum(non_bg, axis=0)
+    col_cutoff = col_cutoff_perc * np.max(colsum)
+
+    fgcol = np.argwhere(colsum > col_cutoff)
+    leftbound = np.min(fgcol)
+    rightbound = np.max(fgcol)
+
+    rowsum = np.sum(non_bg[:,leftbound:rightbound+1], axis=1)
+    row_cutoff = row_cutoff_perc * np.max(rowsum)
+    fgrow = np.argwhere(rowsum > row_cutoff)
+
+    topbound = np.min(fgrow)
+    bottombound = np.max(fgrow)
+
+    return np.array([[topbound, bottombound], [leftbound, rightbound]])
+
+def _get_overlay_mask(npimg, version=None, threshold=250):
+    '''
+    Identify the overlay and annotation elements by brightness (threshold) and color.
+    
+    Parameters
+    ----------
+    npimg : ndarray
+        MxNx3
+    version : None
+        reserved for future use
+    threshold : int
+        brightness threshold for overlay elements, if too low, B-mode data may mistakenly masked
+        
+    Returns
+    -------
+    ndarray
+        False is background and True is overlay
+    '''
+    mask = np.logical_not((npimg[:,:,0] == npimg[:,:,1]) & (npimg[:,:,1] == npimg[:,:,2])) | (npimg[:,:,0] >= threshold)
+    mask = dilation(mask, disk(5))
+    return mask
 
 def preprocess_image(npimg, version=None):
     '''
@@ -68,9 +140,13 @@ def preprocess_image(npimg, version=None):
     itkImage[itk.F,2]
         cropped image scaled to 0.0 to 1.0 (MxN) with physical spacing
     '''
-    spacing, crop = _find_spacing_and_crop(npimg)
-    img = itk.image_from_array((npimg[crop[0,0]:crop[0,1]+1, crop[1,0]:crop[1,1]+1, 0] / 255.0).astype('float32'))
+    mask = _get_overlay_mask(npimg)
+    npimg_clean = _inpaint(npimg[:,:,0], mask)
+    crop = _find_crop(npimg_clean)
+    spacing, ruler_col, ruler_peaks = _get_rightmost_ruler_from_mask(mask)
+    img = itk.image_from_array(npimg_clean[crop[0,0]:crop[0,1]+1, crop[1,0]:crop[1,1]+1])
     img.SetSpacing([spacing, spacing])
+
     return img, { 'spacing' : [spacing, spacing], 'crop' : crop }
 
 def preprocess_video(npvid, framerate=1, version=None):
@@ -83,8 +159,19 @@ def preprocess_video(npvid, framerate=1, version=None):
         reserved for future use
     '''
     npmean = np.mean(npvid, axis=0)
-    spacing, crop = _find_spacing_and_crop(npmean)
-    vid = itk.image_from_array((npvid[:,crop[0,0]:crop[0,1]+1, crop[1,0]:crop[1,1]+1,0] / 255.0).astype('float32').squeeze())
+
+    mask = _get_overlay_mask(npmean)
+    spacing, ruler_col, ruler_peaks = _get_rightmost_ruler_from_mask(mask)
+
+    npmean_clean = _inpaint(npmean[:,:,0], mask)
+    crop = _find_crop(npmean_clean)
+    mask_crop = mask[crop[0,0]:crop[0,1]+1, crop[1,0]:crop[1,1]+1]
+
+    npvid2 = npvid[:,crop[0,0]:crop[0,1]+1, crop[1,0]:crop[1,1]+1,0] / 255.0
+    for j in range(npvid2.shape[0]):
+        npvid2[j] = _inpaint(npvid2[j], mask_crop)
+
+    vid = itk.image_from_array(npvid2)
     vid.SetSpacing([spacing, spacing, framerate])
 
     return vid, { 'spacing' : [spacing, spacing, framerate], 'crop' : crop }

itkpocus/tests/test_butterfly.py

@@ -22,7 +22,7 @@ def test_load_and_preprocess_image(self):
         self.assertAlmostEqual(spacing[0], 0.047619, delta=SPACING_DELTA)
         self.assertAlmostEqual(spacing[1], 0.047619, delta=SPACING_DELTA)
 
-        self.assertAlmostEqual(crop[0,0], 54)
+        self.assertAlmostEqual(crop[0,0], 31)
         self.assertAlmostEqual(crop[0,1], 1048)
         print(crop)
 

itkpocus/tests/test_clarius.py

@@ -15,14 +15,30 @@ def test_load_and_preprocess_image(self):
         #  with self.subTest(i=i):
         #      self.assertEqual()...
 
-        img, meta_dict = itkpocus.clarius.load_and_preprocess_image('./tests/data/clarius_axial-lateral-resolution-2020.png')
+        img, meta_dict = itkpocus.clarius.load_and_preprocess_image('./tests/data/clarius-ocular-image.png')
         spacing = meta_dict['spacing']
         crop = meta_dict['crop']
-        self.assertAlmostEqual(spacing[0], 0.7142857142857143, delta=SPACING_DELTA)
-        self.assertAlmostEqual(spacing[1], 0.7142857142857143, delta=SPACING_DELTA)
+        self.assertAlmostEqual(spacing[0], 0.2564, delta=SPACING_DELTA)
+        self.assertAlmostEqual(spacing[1], 0.2564, delta=SPACING_DELTA)
 
-        self.assertAlmostEqual(crop[0,0], 159)
-        self.assertAlmostEqual(crop[0,1], 800)
+        self.assertAlmostEqual(crop[0,0], 87)
+        self.assertAlmostEqual(crop[0,1], 872)
+        print(crop)
+
+    def test_load_and_preprocess_video(self):
+    # consider using subTests in the future for different versions/images
+    #  with self.subTest(i=i):
+    #      self.assertEqual()...
+
+        img, meta_dict = itkpocus.clarius.load_and_preprocess_video('tests/data/clarius-ocular-short.mp4')
+        spacing = meta_dict['spacing']
+        crop = meta_dict['crop']
+        self.assertAlmostEqual(spacing[0], 0.2564, delta=SPACING_DELTA)
+        self.assertAlmostEqual(spacing[1], 0.2564, delta=SPACING_DELTA)
+        self.assertAlmostEqual(spacing[2], 0.04, delta=SPACING_DELTA)
+
+        self.assertAlmostEqual(crop[0,0], 79)
+        self.assertAlmostEqual(crop[0,1], 881)
         print(crop)
 
 if __name__ == '__main__':