OSCR

Clinical and pathologic correlations of machine learning quantification of Aβ deposits across 3 brain regions of decedents with Alzheimer disease.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match
  1. [1] § METHODS › Statistical analysis ↔ notebook/3) Analysis - Tissue Separation.ipynb, lines 479–551 · score 0.51 · standard deviations, computation, square, mixed, binary, zero

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,223 lines · 46 KB · GPL-3.0 · 1 match

  1. # %% [markdown]
  2. # ### 3 Analysis - Tissue Separation
  3. #
  4. # A sliding window approach was applied on the whole slide images (WSIs) with the trained CNN model to generate confidence heatmaps.
  5. #
  6. # The WSIs were tiled into 1536 x 1536 images patches in the prepocessing steps. A sliding window approached was applied on evey image patch of a WSI. At each time, the CNN model took a 256x256 pixels region as input, forward propagated and generated a prediction score for cored plaque, diffuse plaque and CAA respectively. By systematically sliding the input region across the entire 1536 x 1536 image patch, the prediction scores were saved and ploted as prediction confidence heatmap for this patch. The heatmap for the WSI was obtained by doing this on all image patches of it.
  7. #
  8. # Based on the heatmap prediction scores (saved as npy), we analyzed the plaque density distribution
  10. # %%
  11. import os, glob, datetime
  12. from time import time
  13. from tqdm import tqdm
  15. import numpy as np
  16. import matplotlib
  17. import matplotlib.pyplot as plt
  18. from PIL import Image
  19. Image.MAX_IMAGE_PIXELS = None
  21. import cppyy
  22. cppyy.load_library('gSLICr/build/libgSLICr_lib.so')
  23. cppyy.load_library('/usr/local/cuda-10.0/targets/x86_64-linux/lib/libcudart.so')
  24. cppyy.add_include_path('/usr/local/cuda-10.0/targets/x86_64-linux/include') # Add cuda library
  25. cppyy.include('gSLICr/gSLICr_Lib/gSLICr.h')
  26. from cppyy.gbl import gSLICr
  27. import cppyy.ll
  28. from skimage.segmentation import find_boundaries, mark_boundaries
  29. from skimage.measure import regionprops, label
  30. from skimage.morphology import binary_dilation, convex_hull_image
  32. # %%
  33. IMG_DIR = 'data/outputs/thres_tissue_png/'
  34. SAVE_DIR = 'data/outputs/tissue_mask_png/'
  36. filenames = glob.glob(IMG_DIR + '*.png')
  37. filenames = [filename.split('/')[-1] for filename in filenames]
  38. print(filenames)
  40. # %% [markdown]
  41. # ### Final Pipeline
  43. # %%
  44. # Separate stains
  45. #from skimage.color import separate_stains # Too much memory burden, separate computation tasks
  46. from skimage.util import img_as_float
  47. from tempfile import mkdtemp
  48. import shutil
  50. def color_deconv(rgb, conv_matrix):
  51. """RGB to stain color space conversion using color deconvolution.
  52. Original implementation in scikit-image
  53. https://github.com/scikit-image/scikit-image/blob/81363d32b94763150125dcf2c8909533c37567d9/skimage/color/colorconv.py#L1364
  55. Parameters
  56. ----------
  57. rgb : array_like
  58. The image in RGB format, in a 3-D array of shape ``(.., .., 3)``.
  59. conv_matrix: ndarray
  60. The stain separation matrix as described by G. Landini [1]_.
  61. Returns
  62. -------
  63. out : ndarray
  64. The image in stain color channel, in a 2-D array of shape
  65. ``(.., ..)``.
  66. """
  67. #rgb = img_as_float(rgb, force_copy=True)
  68. rgb = rgb.astype(dtype='float32')
  69. rgb += 2
  70. rgb = -np.log10(rgb, dtype='float32')
  71. out_shape = (rgb.shape[0], rgb.shape[1], 2)
  72. #rgb = np.reshape(rgb, (-1, 3)) @ (conv_matrix[:,0:2]) # 0:2 because only H&E channels are needed
  73. rgb = np.matmul(np.reshape(rgb, (-1, 3)), (conv_matrix[:,0:2]), dtype='float32')
  74. return np.reshape(rgb, out_shape)
  76. # Scale [0,1] to [0,255]
  77. def uint8scale(arr, _min=0, _max=1):
  78. return ((arr - _min) * 1/(_max - _min) * 255).astype('uint8')
  80. # Scale arr to [new_min, new_max]
  81. def scale_range(arr, new_min, new_max):
  82. return (arr - arr.min()) * (new_max - new_min) / (arr.max() - arr.min()) + new_min
  84. # Scale arr to have new_mean and new_std
  85. def scale_meanstd(arr, new_mean, new_std):
  86. return (arr - np.mean(arr)) / np.std(arr) * new_std + new_mean
  88. # Apply contrast threshold using QuPath data
  89. def apply_contrast_threshold(arr, thres_min, thres_max):
  90. arr[arr < thres_min] = thres_min
  91. arr[arr > thres_max] = thres_max
  92. return scale_range(arr, 0, 1)
  94. # Create a memmap numpy array on disk for float64 WSI
  95. # dirpath = mkdtemp()
  96. # fp = os.path.join(dirpath, 'arr.dat')
  97. # print(dirpath)
  99. startT = time()
  101. # Divide the WSI into regular slices to avoid insufficient memory
  102. # Cannot divide WSI into slices due to contrast threshold and scale_meanstd
  103. # slice_size = 10000
  104. rgb_img = np.array(Image.open('data/outputs/norm_png/NA4077-02_AB.png'))
  105. height, width, _ = rgb_img.shape
  106. print(rgb_img.shape)
  108. # Corresponding stain matrix
  109. #stain_mat = np.array([[0.572, 0.589, 0.572], [0.281, 0.514, 0.811], [0.488, -0.804, 0.341]]) # NA3077
  110. stain_mat = np.array([[0.624, 0.662, 0.415], [0.469, 0.612, 0.637], [0.616, -0.743, 0.261]]) # ref: NA5002_2
  111. stain_mat = np.linalg.inv(stain_mat)
  113. # Corresponding scale and contrast threshold constants
  114. consts = []
  115. consts.append({'mean':0.037, 'std':0.12, 'thres_min':-0.045568, 'thres_max':1.051931})
  116. consts.append({'mean':0.091, 'std':0.106, 'thres_min':0.01, 'thres_max':0.3})
  118. # t = tqdm(total=2)
  119. # for i in range(2):
  121. # Separate stains
  122. # i=0: Hematoxylin, i=1: Eosin, i=2: Residual
  123. print('Load Image: %s' % (time()-startT))
  124. mask_img = color_deconv(rgb_img, stain_mat)
  125. del rgb_img
  126. print('Color Deconv: %s' % (time()-startT))
  128. norm_img = []
  129. for i in range(2):
  130. # Scale using mean and standard deviation
  131. mask_img[:,:,i] = scale_meanstd(mask_img[:,:,i], consts[i]['mean'], consts[i]['std'])
  133. # Apply contrast threshold
  134. mask_img[:,:,i] = apply_contrast_threshold(mask_img[:,:,i], consts[i]['thres_min'], consts[i]['thres_max'])
  136. # Convert into uint8
  137. norm_img.append(np.asarray(uint8scale(mask_img[:,:,i]), dtype='uint8'))
  138. del mask_img
  140. # Convert numpy array into PIL image and save to local
  141. # Divide the WSI into regular slices to avoid PIL buffer overflow
  142. slice_size = 30000
  143. iters = np.uint8(np.ceil([height / slice_size, width / slice_size]))
  145. for i in range(2):
  146. mask_img = Image.new('L', (width, height))
  147. for row in range(iters[0]):
  148. for col in range(iters[1]):
  149. # Get start and end pixel location
  150. start_r, start_c = row * slice_size, col * slice_size
  151. end_r, end_c = min(height, start_r + slice_size), min(width, start_c + slice_size)
  153. # Paste the slice mask to image mask
  154. mask_slice = Image.fromarray(norm_img[i][start_r:end_r, start_c:end_c], 'L')
  155. mask_img.paste(mask_slice, (start_c, start_r))
  156. # Save masked image
  157. if i == 0:
  158. mask_img.save('data/outputs/norm_png/NA4077-02_AB_Hemato_4.png')
  159. elif i == 1:
  160. mask_img.save('data/outputs/norm_png/NA4077-02_AB_Eosin_4.png')
  161. del mask_slice, mask_img
  163. print('Total: %s' % (time()-startT))
  166. # %%
  167. # Function definition for each step
  169. # Helper func: Scale [0,1] to [0,255]
  170. def uint8scale(arr, _min=0, _max=1):
  171. """Scale an array to [0,255].
  173. Parameters
  174. -----------
  175. arr : ndarray
  176. Input array.
  177. _min : int, optional
  178. Minimum value within arr.
  179. _max : int, optional
  180. Maximum value within arr.
  182. Returns
  183. -------
  184. output : ndarray
  185. The scaled array.
  186. """
  187. return ((arr - _min) * 1/(_max - _min) * 255).astype('uint8')
  189. def run_slic_seg(img, n_segments=80000, n_iter=5, compactness=0.01, enforce_connectivity=True, slic_zero=True):
  190. """Run gpu SLICO in C++.
  192. Parameters
  193. -----------
  194. img : 2D ndarray
  195. Input image, must be grayscale.
  196. n_segments : int, optional
  197. The (approximate) number of labels in the segmented output image.
  198. n_iter : int, optional
  199. The number of iterations of k-means.
  200. compactness : float, optional
  201. Balances color proximity and space proximity. Higher values give
  202. more weight to space proximity, making superpixel shapes more
  203. square/cubic. In SLICO mode, this is the initial compactness.
  204. This parameter depends strongly on image contrast and on the
  205. shapes of objects in the image. We recommend exploring possible
  206. values on a log scale, e.g., 0.01, 0.1, 1, 10, 100, before
  207. refining around a chosen value.
  208. enforce_connectivity : bool, optional
  209. Whether the generated segments are connected or not.
  210. slic_zero : bool, optional
  211. Run SLIC-zero, the zero-parameter mode of SLIC. [2]_
  213. Returns
  214. -------
  215. segments_slic : 2D ndarray
  216. Integer mask indicating SLIC segment labels.
  218. References
  219. ----------
  220. .. [1] Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi,
  221. Pascal Fua, and Sabine Süsstrunk, SLIC Superpixels Compared to
  222. State-of-the-art Superpixel Methods, TPAMI, May 2012.
  223. .. [2] http://ivrg.epfl.ch/research/superpixels#SLICO
  224. """
  225. height, width = img.shape
  226. # Create a Run_SLIC C++ class
  227. run_slic = gSLICr.Run_SLIC()
  229. # Specify settings
  230. gslic_settings = run_slic.get_settings()
  231. gslic_settings = cppyy.bind_object(gslic_settings, gSLICr.objects.settings)
  232. gslic_settings.img_size.x = width
  233. gslic_settings.img_size.y = height
  234. gslic_settings.no_segs = n_segments
  235. gslic_settings.no_iters = n_iter
  236. gslic_settings.coh_weight = compactness
  237. gslic_settings.do_enforce_connectivity = enforce_connectivity
  238. gslic_settings.slic_zero = slic_zero
  239. gslic_settings.color_space = gSLICr.GRAY
  240. gslic_settings.seg_method = gSLICr.GIVEN_NUM
  241. run_slic.print_settings()
  243. # Start running gSLIC
  244. run_slic.run(img)
  246. # Get the segmentation mask
  247. segments_slic = run_slic.get_mask()
  248. segments_slic.reshape((height*width,))
  249. segments_slic = np.frombuffer(segments_slic, dtype=np.intc, count=height*width).reshape((height, width)).copy()
  251. # Destruct the C++ class
  252. run_slic.__destruct__()
  253. return segments_slic
  255. def slic_mean_intensity_thres(img, segments_slic, thres=60):
  256. """Thresholding based on SLIC segments' mean intensity.
  258. Parameters
  259. -----------
  260. img : 2D ndarray
  261. Input image, must be grayscale.
  262. segments_slic : 2D ndarray
  263. Integer mask indicating SLIC segment labels.
  264. thres : int, optional
  265. Threshold value of mean intensity.
  267. Returns
  268. -------
  269. mask_img : 2D ndarray
  270. Binary integer mask indicating segment labels.
  271. """
  272. mask_img = np.zeros_like(img, dtype='uint8')
  273. for region in regionprops(segments_slic, img):
  274. if region.mean_intensity >= thres:
  275. mask_img.flat[np.ravel_multi_index(region.coords.transpose(), mask_img.shape)] = 255
  276. return mask_img
  278. def get_connected_comp(mask_img, structure_size=9):
  279. """Get connected components inside an image.
  281. Parameters
  282. -----------
  283. mask_img : 2D ndarray
  284. Input image, must be binary.
  285. structure_size : int, optional
  286. Structure size for binary dilation before connected component labeling.
  288. Returns
  289. -------
  290. output : 2D ndarray
  291. Integer mask indicating connected component labels.
  292. """
  293. return label(binary_dilation(mask_img, selem=np.ones((structure_size, structure_size), dtype='int')), connectivity=2)
  295. def find_tissue(img, connected_comp, prop=0.2):
  296. """Find the tissue component within a connected component map.
  298. Parameters
  299. -----------
  300. img : 2D ndarray
  301. Input image, must be grayscale.
  302. connected_comp : 2D ndarray
  303. Integer mask indicating connected component labels.
  304. prop : float, optional
  305. Proportional occupied area of the entire image to be considered as tissue.
  307. Returns
  308. -------
  309. tissue_img : 2D ndarray
  310. Binary integer mask indicating the tissue.
  311. """
  312. height, width = img.shape
  313. area_thres = prop * height * width
  314. tissue_img = np.zeros_like(img, dtype='uint8')
  315. for region in regionprops(connected_comp, img):
  316. if region.area >= area_thres:
  317. tissue_img.flat[np.ravel_multi_index(region.coords.transpose(), tissue_img.shape)] = 255
  318. return tissue_img
  320. def find_refine_boundaries(img, segments_slic, tissue_img):
  321. """Find and refine boundaries of the tissue image using convex hull.
  323. Parameters
  324. -----------
  325. img : 2D ndarray
  326. Input image, must be grayscale.
  327. segments_slic : 2D ndarray
  328. Integer mask indicating SLIC segment labels.
  329. tissue_img : 2D ndarray
  330. Binary integer mask indicating the tissue.
  332. Returns
  333. -------
  334. tissue_img : 2D ndarray
  335. Binary integer refined mask indicating the tissue.
  336. """
  337. # Find boundaries
  338. boundaries = find_boundaries(tissue_img, connectivity=2, mode='thick')
  339. labels = np.unique(segments_slic[boundaries])
  340. # Convex hull on boundaries
  341. for region in regionprops(segments_slic, img):
  342. if region.label in labels:
  343. min_r, min_c, max_r, max_c = region.bbox
  344. img_tile = tissue_img[min_r:max_r, min_c:max_c]
  345. outlier_hull_img = convex_hull_image(img_tile)
  346. tissue_img[min_r:max_r, min_c:max_c] = (outlier_hull_img * 255).astype('uint8')
  347. return tissue_img
  349. def save_grayscale_img(img, filename):
  350. """Save a grayscale image locally as filename.
  352. Parameters
  353. -----------
  354. img : 2D ndarray
  355. Input image, must be grayscale.
  356. filename : string
  357. File path and file name to save.
  358. """
  359. height, width = img.shape
  360. # Convert numpy array into PIL image and save to local
  361. # Divide the WSI into regular slices to avoid PIL buffer overflow
  362. slice_size = 30000
  363. iters = np.uint8(np.ceil([height / slice_size, width / slice_size]))
  365. save_img = Image.new('L', (width, height))
  366. for row in range(iters[0]):
  367. for col in range(iters[1]):
  368. # Get start and end pixel location
  369. start_r, start_c = row * slice_size, col * slice_size
  370. end_r, end_c = min(height, start_r + slice_size), min(width, start_c + slice_size)
  372. # Paste the slice mask to image mask
  373. save_slice = Image.fromarray(img[start_r:end_r, start_c:end_c], 'L')
  374. save_img.paste(save_slice, (start_c, start_r))
  375. # Save masked image
  376. save_img.save(filename)
  378. # %%
  379. # Final Pipeline
  380. currentDT = datetime.datetime.now(); print(str(currentDT)); totalT = 0.0; startT = time()
  382. img = np.array(Image.open('data/outputs/norm_png/NA3777-02_AB_Eosin.png'))
  383. # Binary Thresholding to improve contrast
  384. img = np.where(img < 15, 0, 255).astype('uint8')
  385. print("Loaded Image: %s" % (time()-startT)); totalT += (time()-startT); startT = time()
  387. segments_slic = run_slic_seg(img)
  388. print("Total SLIC: %s" % (time()-startT)); totalT += (time()-startT); startT = time()
  390. # Mean Intensity Thresholding
  391. mask_img = slic_mean_intensity_thres(img, segments_slic)
  392. print("Mean Intensity Thresholding: %s" % (time()-startT)); totalT += (time()-startT); startT = time()
  394. # Connected components
  395. connected_comp = get_connected_comp(mask_img)
  396. del mask_img
  397. print("Connected component labeling: %s" % (time()-startT)); totalT += (time()-startT); startT = time()
  399. # Find tissue component
  400. tissue_img = find_tissue(img, connected_comp)
  401. del connected_comp
  402. print("Find tissue component: %s" % (time()-startT)); totalT += (time()-startT); startT = time()
  404. # Find and refine boundaries
  405. tissue_img = find_refine_boundaries(img, segments_slic, tissue_img)
  406. del img, segments_slic
  407. print("Find and refine boundaries: %s" % (time()-startT)); totalT += (time()-startT); startT = time()
  409. # Save image
  410. save_grayscale_img(tissue_img, 'NA3777-02_AB_Eosin_SLICO_0.01_iter=5_80000_finalmask.png')
  411. print("Saved: %s" % (time()-startT)); totalT += (time()-startT)
  412. print("Total: %s" % totalT)
  414. # %% [markdown]
  415. # ### Timing Results
  416. # | SLIC #Segments | Load Image | SLIC | Mean Intensity Thresholding | Connected Component | Find Tissue | Find Boundaries | Convex Hull | Save Output | Total |
  417. # |:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|
  418. # | 40000 | 27.6471 | 232.6231 | 57.3725 | 106.8111 | 58.1894 | 137.0627 | 33.8033 | 48.8813 | 702.3917 |
  419. # | 60000 | 30.0168 | 206.1811 | 60.3673 | 110.6402 | 59.3389 | 140.0935 | 31.8562 | 49.2528 | 687.7481 |
  420. # | 80000 | 27.9608 | 249.4539 | 58.0493 | 107.5019 | 57.5752 | 136.3321 | 30.8568 | 49.1053 | 716.8363 |
  422. # %% [markdown]
  423. # ## Pipeline Development
  424. # ### Latest Scikit-image Canny Algorithm
  426. # %%
  427. """
  428. canny.py - Canny Edge detector
  430. Reference: Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
  431. Pattern Analysis and Machine Intelligence, 8:679-714, 1986
  433. Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
  434. Website: http://www.cellprofiler.org
  435. Copyright (c) 2003-2009 Massachusetts Institute of Technology
  436. Copyright (c) 2009-2011 Broad Institute
  437. All rights reserved.
  438. Original author: Lee Kamentsky
  439. """
  441. import numpy as np
  442. import scipy.ndimage as ndi
  443. from scipy.ndimage import generate_binary_structure
  444. from skimage.filters import gaussian
  445. from skimage import dtype_limits, img_as_float
  446. from skimage._shared.utils import assert_nD
  449. def smooth_with_function_and_mask(image, function, mask):
  450. """Smooth an image with a linear function, ignoring masked pixels
  452. Parameters
  453. ----------
  454. image : array
  455. Image you want to smooth.
  456. function : callable
  457. A function that does image smoothing.
  458. mask : array
  459. Mask with 1's for significant pixels, 0's for masked pixels.
  461. Notes
  462. ------
  463. This function calculates the fractional contribution of masked pixels
  464. by applying the function to the mask (which gets you the fraction of
  465. the pixel data that's due to significant points). We then mask the image
  466. and apply the function. The resulting values will be lower by the
  467. bleed-over fraction, so you can recalibrate by dividing by the function
  468. on the mask to recover the effect of smoothing from just the significant
  469. pixels.
  470. """
  471. bleed_over = function(mask.astype(float))
  472. masked_image = np.zeros(image.shape, image.dtype)
  473. masked_image[mask] = image[mask]
  474. smoothed_image = function(masked_image)
  475. output_image = smoothed_image / (bleed_over + np.finfo(float).eps)
  476. return output_image
  479. def canny(image, sigma=1., low_threshold=None, high_threshold=None, mask=None,
  480. use_quantiles=False):
  481. """Edge filter an image using the Canny algorithm.
  483. Parameters
  484. -----------
  485. image : 2D array
  486. Grayscale input image to detect edges on; can be of any dtype.
  487. sigma : float
  488. Standard deviation of the Gaussian filter.
  489. low_threshold : float
  490. Lower bound for hysteresis thresholding (linking edges).
  491. If None, low_threshold is set to 10% of dtype's max.
  492. high_threshold : float
  493. Upper bound for hysteresis thresholding (linking edges).
  494. If None, high_threshold is set to 20% of dtype's max.
  495. mask : array, dtype=bool, optional
  496. Mask to limit the application of Canny to a certain area.
  497. use_quantiles : bool, optional
  498. If True then treat low_threshold and high_threshold as quantiles of the
  499. edge magnitude image, rather than absolute edge magnitude values. If True
  500. then the thresholds must be in the range [0, 1].
  502. Returns
  503. -------
  504. output : 2D array (image)
  505. The binary edge map.
  507. See also
  508. --------
  509. skimage.sobel
  511. Notes
  512. -----
  513. The steps of the algorithm are as follows:
  515. * Smooth the image using a Gaussian with ``sigma`` width.
  517. * Apply the horizontal and vertical Sobel operators to get the gradients
  518. within the image. The edge strength is the norm of the gradient.
  520. * Thin potential edges to 1-pixel wide curves. First, find the normal
  521. to the edge at each point. This is done by looking at the
  522. signs and the relative magnitude of the X-Sobel and Y-Sobel
  523. to sort the points into 4 categories: horizontal, vertical,
  524. diagonal and antidiagonal. Then look in the normal and reverse
  525. directions to see if the values in either of those directions are
  526. greater than the point in question. Use interpolation to get a mix of
  527. points instead of picking the one that's the closest to the normal.
  529. * Perform a hysteresis thresholding: first label all points above the
  530. high threshold as edges. Then recursively label any point above the
  531. low threshold that is 8-connected to a labeled point as an edge.
  533. References
  534. -----------
  535. .. [1] Canny, J., A Computational Approach To Edge Detection, IEEE Trans.
  536. Pattern Analysis and Machine Intelligence, 8:679-714, 1986
  537. .. [2] William Green's Canny tutorial
  538. http://dasl.unlv.edu/daslDrexel/alumni/bGreen/www.pages.drexel.edu/_weg22/can_tut.html
  540. Examples
  541. --------
  542. >>> from skimage import feature
  543. >>> # Generate noisy image of a square
  544. >>> im = np.zeros((256, 256))
  545. >>> im[64:-64, 64:-64] = 1
  546. >>> im += 0.2 * np.random.rand(*im.shape)
  547. >>> # First trial with the Canny filter, with the default smoothing
  548. >>> edges1 = feature.canny(im)
  549. >>> # Increase the smoothing for better results
  550. >>> edges2 = feature.canny(im, sigma=3)
  551. """
  553. #
  554. # The steps involved:
  555. #
  556. # * Smooth using the Gaussian with sigma above.
  557. #
  558. # * Apply the horizontal and vertical Sobel operators to get the gradients
  559. # within the image. The edge strength is the sum of the magnitudes
  560. # of the gradients in each direction.
  561. #
  562. # * Find the normal to the edge at each point using the arctangent of the
  563. # ratio of the Y sobel over the X sobel - pragmatically, we can
  564. # look at the signs of X and Y and the relative magnitude of X vs Y
  565. # to sort the points into 4 categories: horizontal, vertical,
  566. # diagonal and antidiagonal.
  567. #
  568. # * Look in the normal and reverse directions to see if the values
  569. # in either of those directions are greater than the point in question.
  570. # Use interpolation to get a mix of points instead of picking the one
  571. # that's the closest to the normal.
  572. #
  573. # * Label all points above the high threshold as edges.
  574. # * Recursively label any point above the low threshold that is 8-connected
  575. # to a labeled point as an edge.
  576. #
  577. # Regarding masks, any point touching a masked point will have a gradient
  578. # that is "infected" by the masked point, so it's enough to erode the
  579. # mask by one and then mask the output. We also mask out the border points
  580. # because who knows what lies beyond the edge of the image?
  581. #
  582. assert_nD(image, 2)
  583. dtype_max = dtype_limits(image, clip_negative=False)[1]
  585. if low_threshold is None:
  586. low_threshold = 0.1
  587. else:
  588. low_threshold = low_threshold / dtype_max
  590. if high_threshold is None:
  591. high_threshold = 0.2
  592. else:
  593. high_threshold = high_threshold / dtype_max
  595. if mask is None:
  596. mask = np.ones(image.shape, dtype=bool)
  598. def fsmooth(x):
  599. return img_as_float(gaussian(x, sigma, mode='constant'))
  601. smoothed = smooth_with_function_and_mask(image, fsmooth, mask)
  602. jsobel = ndi.sobel(smoothed, axis=1)
  603. isobel = ndi.sobel(smoothed, axis=0)
  604. abs_isobel = np.abs(isobel)
  605. abs_jsobel = np.abs(jsobel)
  606. magnitude = np.hypot(isobel, jsobel)
  608. #
  609. # Make the eroded mask. Setting the border value to zero will wipe
  610. # out the image edges for us.
  611. #
  612. s = generate_binary_structure(2, 2)
  613. eroded_mask = ndi.binary_erosion(mask, s, border_value=0)
  614. eroded_mask = eroded_mask & (magnitude > 0)
  615. #
  616. #--------- Find local maxima --------------
  617. #
  618. # Assign each point to have a normal of 0-45 degrees, 45-90 degrees,
  619. # 90-135 degrees and 135-180 degrees.
  620. #
  621. local_maxima = np.zeros(image.shape, bool)
  622. #----- 0 to 45 degrees ------
  623. pts_plus = (isobel >= 0) & (jsobel >= 0) & (abs_isobel >= abs_jsobel)
  624. pts_minus = (isobel <= 0) & (jsobel <= 0) & (abs_isobel >= abs_jsobel)
  625. pts = pts_plus | pts_minus
  626. pts = eroded_mask & pts
  627. # Get the magnitudes shifted left to make a matrix of the points to the
  628. # right of pts. Similarly, shift left and down to get the points to the
  629. # top right of pts.
  630. c1 = magnitude[1:, :][pts[:-1, :]]
  631. c2 = magnitude[1:, 1:][pts[:-1, :-1]]
  632. m = magnitude[pts]
  633. w = abs_jsobel[pts] / abs_isobel[pts]
  634. c_plus = c2 * w + c1 * (1 - w) <= m
  635. c1 = magnitude[:-1, :][pts[1:, :]]
  636. c2 = magnitude[:-1, :-1][pts[1:, 1:]]
  637. c_minus = c2 * w + c1 * (1 - w) <= m
  638. local_maxima[pts] = c_plus & c_minus
  639. #----- 45 to 90 degrees ------
  640. # Mix diagonal and vertical
  641. #
  642. pts_plus = (isobel >= 0) & (jsobel >= 0) & (abs_isobel <= abs_jsobel)
  643. pts_minus = (isobel <= 0) & (jsobel <= 0) & (abs_isobel <= abs_jsobel)
  644. pts = pts_plus | pts_minus
  645. pts = eroded_mask & pts
  646. c1 = magnitude[:, 1:][pts[:, :-1]]
  647. c2 = magnitude[1:, 1:][pts[:-1, :-1]]
  648. m = magnitude[pts]
  649. w = abs_isobel[pts] / abs_jsobel[pts]
  650. c_plus = c2 * w + c1 * (1 - w) <= m
  651. c1 = magnitude[:, :-1][pts[:, 1:]]
  652. c2 = magnitude[:-1, :-1][pts[1:, 1:]]
  653. c_minus = c2 * w + c1 * (1 - w) <= m
  654. local_maxima[pts] = c_plus & c_minus
  655. #----- 90 to 135 degrees ------
  656. # Mix anti-diagonal and vertical
  657. #
  658. pts_plus = (isobel <= 0) & (jsobel >= 0) & (abs_isobel <= abs_jsobel)
  659. pts_minus = (isobel >= 0) & (jsobel <= 0) & (abs_isobel <= abs_jsobel)
  660. pts = pts_plus | pts_minus
  661. pts = eroded_mask & pts
  662. c1a = magnitude[:, 1:][pts[:, :-1]]
  663. c2a = magnitude[:-1, 1:][pts[1:, :-1]]
  664. m = magnitude[pts]
  665. w = abs_isobel[pts] / abs_jsobel[pts]
  666. c_plus = c2a * w + c1a * (1.0 - w) <= m
  667. c1 = magnitude[:, :-1][pts[:, 1:]]
  668. c2 = magnitude[1:, :-1][pts[:-1, 1:]]
  669. c_minus = c2 * w + c1 * (1.0 - w) <= m
  670. local_maxima[pts] = c_plus & c_minus
  671. #----- 135 to 180 degrees ------
  672. # Mix anti-diagonal and anti-horizontal
  673. #
  674. pts_plus = (isobel <= 0) & (jsobel >= 0) & (abs_isobel >= abs_jsobel)
  675. pts_minus = (isobel >= 0) & (jsobel <= 0) & (abs_isobel >= abs_jsobel)
  676. pts = pts_plus | pts_minus
  677. pts = eroded_mask & pts
  678. c1 = magnitude[:-1, :][pts[1:, :]]
  679. c2 = magnitude[:-1, 1:][pts[1:, :-1]]
  680. m = magnitude[pts]
  681. w = abs_jsobel[pts] / abs_isobel[pts]
  682. c_plus = c2 * w + c1 * (1 - w) <= m
  683. c1 = magnitude[1:, :][pts[:-1, :]]
  684. c2 = magnitude[1:, :-1][pts[:-1, 1:]]
  685. c_minus = c2 * w + c1 * (1 - w) <= m
  686. local_maxima[pts] = c_plus & c_minus
  688. #
  689. #---- If use_quantiles is set then calculate the thresholds to use
  690. #
  691. if use_quantiles:
  692. if high_threshold > 1.0 or low_threshold > 1.0:
  693. raise ValueError("Quantile thresholds must not be > 1.0")
  694. if high_threshold < 0.0 or low_threshold < 0.0:
  695. raise ValueError("Quantile thresholds must not be < 0.0")
  697. high_threshold = np.percentile(magnitude, 100.0 * high_threshold)
  698. low_threshold = np.percentile(magnitude, 100.0 * low_threshold)
  700. #
  701. #---- Create two masks at the two thresholds.
  702. #
  703. high_mask = local_maxima & (magnitude >= high_threshold)
  704. low_mask = local_maxima & (magnitude >= low_threshold)
  706. #
  707. # Segment the low-mask, then only keep low-segments that have
  708. # some high_mask component in them
  709. #
  710. strel = np.ones((3, 3), bool)
  711. labels, count = ndi.label(low_mask, strel)
  712. if count == 0:
  713. return low_mask
  715. sums = (np.array(ndi.sum(high_mask, labels,
  716. np.arange(count, dtype=np.int32) + 1),
  717. copy=False, ndmin=1))
  718. good_label = np.zeros((count + 1,), bool)
  719. good_label[1:] = sums > 0
  720. output_mask = good_label[labels]
  721. return output_mask
  723. # %% [markdown]
  724. # ### Superpixel
  726. # %%
  727. # Superpixel
  728. from skimage.segmentation import slic, mark_boundaries
  729. from skimage.util import img_as_float
  731. img = np.array(Image.open('data/outputs/norm_png/NA3777-02_AB_Eosin.png'))
  732. # Binary Thresholding
  733. img = np.where(img < 15, 0, 255).astype('uint8')
  734. height, width = img.shape
  735. print(img.shape)
  737. # %%
  738. import cppyy
  739. cppyy.load_library('gSLICr/build/libgSLICr_lib.so')
  740. cppyy.load_library('/usr/local/cuda-10.0/targets/x86_64-linux/lib/libcudart.so')
  741. cppyy.add_include_path('/usr/local/cuda-10.0/targets/x86_64-linux/include') # Add cuda library
  742. cppyy.include('gSLICr/gSLICr_Lib/gSLICr.h')
  743. from cppyy.gbl import gSLICr
  744. import cppyy.ll
  745. from skimage.segmentation import mark_boundaries
  747. currentDT = datetime.datetime.now(); print(str(currentDT)); startT = time()
  748. # Create a Run_SLIC C++ class
  749. run_slic = gSLICr.Run_SLIC()
  750. print("Run_SLIC: %s" % (time()-startT))
  752. # Specify settings
  753. gslic_settings = run_slic.get_settings()
  754. gslic_settings = cppyy.bind_object(gslic_settings, gSLICr.objects.settings)
  755. gslic_settings.img_size.x = width
  756. gslic_settings.img_size.y = height
  757. gslic_settings.no_segs = 60000
  758. gslic_settings.no_iters = 5
  759. gslic_settings.coh_weight = 0.01
  760. gslic_settings.do_enforce_connectivity = True
  761. gslic_settings.slic_zero = True
  762. gslic_settings.color_space = gSLICr.GRAY
  763. gslic_settings.seg_method = gSLICr.GIVEN_NUM
  764. run_slic.print_settings()
  766. # Start running gSLIC
  767. run_slic.run(img)
  768. print("run: %s" % (time()-startT))
  770. # Get the segmentation mask
  771. segments_slic = run_slic.get_mask()
  772. print("get_mask: %s" % (time()-startT))
  774. segments_slic.reshape((height*width,))
  775. segments_slic = np.frombuffer(segments_slic, dtype=np.intc, count=height*width).reshape((height, width)).copy()
  776. print("reshape: %s" % (time()-startT))
  778. run_slic.__destruct__()
  779. print("__destruct__: %s" % (time()-startT))
  780. # 81 to 73 after block-dim=32 using slic iter=0
  781. # Down to 50 after changing from Float3Image to FloatImage
  783. np.save('segments_slico_0.01_iter=5_60000', segments_slic)
  784. print("save as npy: %s" % (time()-startT))
  786. # Scale [0,1] to [0,255]
  787. def uint8scale(arr, _min=0, _max=1):
  788. return ((arr - _min) * 1/(_max - _min) * 255).astype('uint8')
  790. print('SLIC number of segments: {}'.format(len(np.unique(segments_slic))))
  791. mark_img = mark_boundaries(img, segments_slic, outline_color=(1,1,0), mode='thick')
  792. print('Marked')
  793. del segments_slic#, img
  794. norm_img = uint8scale(mark_img)
  795. print('Normalized')
  796. del mark_img
  798. save_img = Image.fromarray(norm_img, 'RGB')
  799. save_img.save('NA3777-02_AB_Eosin_SLICO_0.01_iter=5_60000.png')
  800. del save_img
  801. print('Saved')
  803. # %%
  804. # Superpixel
  805. from skimage.segmentation import slic, mark_boundaries
  806. from skimage.util import img_as_float
  808. # Scale [0,1] to [0,255]
  809. def uint8scale(arr, _min=0, _max=1):
  810. return ((arr - _min) * 1/(_max - _min) * 255).astype('uint8')
  812. img = np.array(Image.open('data/outputs/norm_png/NA3777-02_AB_Eosin.png'))
  813. # Binary Thresholding
  814. img = np.where(img < 15, 0, 255).astype('uint8')
  815. height, width = img.shape
  816. print(img.shape)
  818. img = img_as_float(img)
  820. # Max memory 123+16 using float64. 255 segments. Run time: 1324 seconds = 22 minutes
  821. currentDT = datetime.datetime.now(); print(str(currentDT)); startT = time()
  822. segments_slic = slic(img, n_segments=5000, compactness=0.01, slic_zero=True)
  823. print("SLIC: %s" % (time()-startT))
  824. print('SLIC number of segments: {}'.format(len(np.unique(segments_slic))))
  826. del img
  827. img = np.array(Image.open('data/outputs/norm_png/NA3777-02_AB_Eosin.png'))
  828. # Binary Thresholding
  829. img = np.where(img < 15, 0, 255).astype('uint8')
  830. print('Loaded')
  831. mark_img = mark_boundaries(img, segments_slic, outline_color=(1,1,0), mode='thick')
  832. print('Marked')
  833. del segments_slic, img
  834. norm_img = uint8scale(mark_img)
  835. print('Normalized')
  836. del mark_img
  838. save_img = Image.fromarray(norm_img, 'RGB')
  839. save_img.save('NA3777-02_AB_Eosin_SLIC.png')
  841. print('Saved')
  842. print("Total: %s" % (time()-startT))
  844. # %% [markdown]
  845. # ### Mean Intensity Threshold
  847. # %%
  848. segments_slic = np.load('segments_slico_0.01_iter=5_60000.npy')
  849. img = np.array(Image.open('data/outputs/norm_png/NA3777-02_AB_Eosin.png'))
  850. # Binary Thresholding
  851. img = np.where(img < 15, 0, 255).astype('uint8')
  852. height, width = img.shape
  854. # %%
  855. from skimage.measure import regionprops
  857. mean_intensity_thres = 60
  858. mask_img = np.zeros_like(img, dtype='uint8')
  859. for region in regionprops(segments_slic, img):
  860. if region.mean_intensity >= mean_intensity_thres:
  861. mask_img.flat[np.ravel_multi_index(region.coords.transpose(), mask_img.shape)] = 255
  863. # %%
  864. # Save superpixel with average intensity
  865. from skimage.measure import regionprops
  866. segments_slic = np.load('segments_slico_0.1_iter=5_40000.npy')
  867. img = np.array(Image.open('data/outputs/norm_png/NA3777-02_AB_Eosin.png'))
  868. # Binary Thresholding
  869. img = np.where(img < 15, 0, 255).astype('uint8')
  870. height, width = img.shape
  872. mask_img = np.zeros_like(img, dtype='uint8')
  873. for region in regionprops(segments_slic, img):
  874. mask_img.flat[np.ravel_multi_index(region.coords.transpose(), mask_img.shape)] = region.mean_intensity
  876. # Convert numpy array into PIL image and save to local
  877. # Divide the WSI into regular slices to avoid PIL buffer overflow
  878. slice_size = 30000
  879. iters = np.uint8(np.ceil([height / slice_size, width / slice_size]))
  881. save_img = Image.new('L', (width, height))
  882. for row in range(iters[0]):
  883. for col in range(iters[1]):
  884. # Get start and end pixel location
  885. start_r, start_c = row * slice_size, col * slice_size
  886. end_r, end_c = min(height, start_r + slice_size), min(width, start_c + slice_size)
  888. # Paste the slice mask to image mask
  889. save_slice = Image.fromarray(mask_img[start_r:end_r, start_c:end_c], 'L')
  890. save_img.paste(save_slice, (start_c, start_r))
  891. # Save masked image
  892. save_img.save('NA3777-02_AB_Eosin_SLIC_0.1_iter=5_40000_meanintensity.png')
  894. # %% [markdown]
  895. # ### Connected Component
  897. # %%
  898. from skimage.measure import label
  899. from skimage.morphology import binary_dilation
  901. currentDT = datetime.datetime.now(); print(str(currentDT)); startT = time()
  902. connected_comp = label(binary_dilation(mask_img, selem=np.ones((9, 9), dtype='int')), connectivity=2)
  903. print("Connected component labeling: %s" % (time()-startT))
  905. print('Number of connected components: {}'.format(len(np.unique(connected_comp))))
  907. # %%
  908. for region in regionprops(connected_comp, img):
  909. print(region.area)
  911. # %%
  912. mask_img2 = np.zeros_like(img, dtype='uint8')
  913. for region in regionprops(connected_comp, img):
  914. if region.area >= 1e8:
  915. mask_img2.flat[np.ravel_multi_index(region.coords.transpose(), mask_img2.shape)] = 255
  917. # %%
  918. # Save connected component image
  919. # Convert numpy array into PIL image and save to local
  920. # Divide the WSI into regular slices to avoid PIL buffer overflow
  921. slice_size = 30000
  922. iters = np.uint8(np.ceil([height / slice_size, width / slice_size]))
  924. save_img = Image.new('L', (width, height))
  925. for row in range(iters[0]):
  926. for col in range(iters[1]):
  927. # Get start and end pixel location
  928. start_r, start_c = row * slice_size, col * slice_size
  929. end_r, end_c = min(height, start_r + slice_size), min(width, start_c + slice_size)
  931. # Paste the slice mask to image mask
  932. save_slice = Image.fromarray(mask_img2[start_r:end_r, start_c:end_c], 'L')
  933. save_img.paste(save_slice, (start_c, start_r))
  934. # Save masked image
  935. save_img.save('NA3777-02_AB_Eosin_SLIC_0.01_iter=5_40000_connectedcomp=25.png')
  937. # %% [markdown]
  938. # ### Find and Extend Boundaries
  940. # %%
  941. from skimage.segmentation import find_boundaries
  943. boundaries = find_boundaries(mask_img2, connectivity=2, mode='thick')
  945. # %% [markdown]
  946. # ### Associate with SLICO labels
  948. # %%
  949. labels = np.unique(segments_slic[boundaries])
  950. print(len(labels))
  952. # %%
  953. # Keep only medium intensity labels
  954. lower_thres = 60
  955. edge_labels = []
  956. for region in regionprops(segments_slic, img):
  957. mean_inten = region.mean_intensity
  958. if region.label in labels and mean_inten >= lower_thres:
  959. edge_labels.append(region.label)
  960. print(len(edge_labels))
  962. # %%
  963. # Mark boundaries for visualization
  964. del mask_img, mask_img2, connected_comp, boundaries
  966. from skimage.segmentation import mark_boundaries
  968. mark_img = mark_boundaries(img, segments_slic, mode='thick')
  969. # Scale to uint8
  970. mark_img *= 255
  971. mark_img = mark_img.astype('uint8')
  973. # %%
  974. # dev
  975. import matplotlib.patches as patches
  976. # from skimage.feature import canny
  977. from skimage.morphology import convex_hull_image
  978. from IPython import display
  980. #from pyod.models.ocsvm import OCSVM
  981. #from pyod.models.cof import COF
  982. #from pyod.models.abod import ABOD
  983. #from pyod.models.lscp import LSCP
  984. #from pyod.models.loci import LOCI
  986. from pyod.models.knn import KNN
  987. from pyod.models.iforest import IForest
  988. from pyod.models.cblof import CBLOF
  989. from pyod.models.feature_bagging import FeatureBagging
  990. from pyod.models.lof import LOF
  991. from pyod.models.hbos import HBOS
  992. from skimage.morphology import binary_erosion
  993. from pyod.models.pca import PCA
  994. from pyod.models.mcd import MCD
  996. nrows = 10
  997. ncols = 6
  998. fig, axes = plt.subplots(nrows, ncols, figsize=(16,20), dpi=300)
  999. axes = axes.reshape((-1,))
  1000. canny_times = np.zeros(nrows)
  1001. outlier_times = np.zeros(nrows)
  1002. count = 0
  1004. for region in regionprops(segments_slic, img):
  1005. #if region.label in edge_labels:
  1006. if region.label in labels:
  1007. min_r, min_c, max_r, max_c = region.bbox
  1008. img_tile = img[min_r:max_r, min_c:max_c]
  1009. #axes[count%nrows*ncols].imshow(img_tile, cmap='gray')
  1010. axes[count%nrows*ncols].imshow(mask_img[min_r:max_r, min_c:max_c])
  1011. axes[count%nrows*ncols].set_title('Region {} Label {}'.format(count+1, region.label))
  1012. axes[count%nrows*ncols].set_xticklabels([])
  1013. axes[count%nrows*ncols].set_yticklabels([])
  1015. # Rough location within WSI
  1016. rect = patches.Rectangle((min_c, min_r),max_c-min_c,max_r-min_r,linewidth=10,edgecolor='r',facecolor='r')
  1017. axes[count%nrows*ncols+1].clear()
  1018. axes[count%nrows*ncols+1].add_patch(rect)
  1019. axes[count%nrows*ncols+1].set_xlim((0, img.shape[1]))
  1020. axes[count%nrows*ncols+1].set_ylim((0, img.shape[0]))
  1021. axes[count%nrows*ncols+1].invert_yaxis()
  1022. axes[count%nrows*ncols+1].set_title('Approx Location')
  1023. axes[count%nrows*ncols+1].set_xticklabels([])
  1024. axes[count%nrows*ncols+1].set_yticklabels([])
  1026. # # Canny edge detection
  1027. # canny_times[count%nrows] = time()
  1028. # canny_img = canny(img_tile, sigma=7)
  1029. # canny_times[count%nrows] = time() - canny_times[count%nrows]
  1030. # axes[count%nrows*ncols+2].imshow(canny_img, cmap='gray')
  1031. # axes[count%nrows*ncols+2].set_title('Canny')
  1032. # axes[count%nrows*ncols+2].set_xticklabels([])
  1033. # axes[count%nrows*ncols+2].set_yticklabels([])
  1035. # canny_hull_img = convex_hull_image(canny_img)
  1036. # axes[count%nrows*ncols+3].imshow(canny_hull_img, cmap='gray')
  1037. # axes[count%nrows*ncols+3].set_title('Canny Hull')
  1038. # axes[count%nrows*ncols+3].set_xticklabels([])
  1039. # axes[count%nrows*ncols+3].set_yticklabels([])
  1041. # Outlier removal and convex hull
  1042. canny_times[count%nrows] = time()
  1043. #X = np.column_stack(np.where(img_tile==255)) # (row, col) indices
  1044. X = np.column_stack(np.where(binary_erosion(img_tile, selem=np.ones((5,5))))) # (row, col) indices
  1046. model = HBOS(contamination=0.03, n_bins=100, tol=0.5) # 0.01s Very Good and doesn't eat edges
  1047. model.fit(X)
  1048. outlier_img = np.zeros_like(img_tile, dtype='uint8')
  1049. outlier_img.flat[np.ravel_multi_index(X[np.invert(model.labels_.astype('bool'))].transpose(), outlier_img.shape)] = 255
  1050. outlier_hull_img = convex_hull_image(outlier_img)
  1051. canny_times[count%nrows] = time() - canny_times[count%nrows]
  1053. axes[count%nrows*ncols+2].imshow(outlier_img, cmap='gray')
  1054. axes[count%nrows*ncols+2].set_title('HBOS Outlier Removal')
  1055. axes[count%nrows*ncols+2].set_xticklabels([])
  1056. axes[count%nrows*ncols+2].set_yticklabels([])
  1058. axes[count%nrows*ncols+3].imshow(outlier_hull_img, cmap='gray')
  1059. axes[count%nrows*ncols+3].set_title('Outlier Hull')
  1060. axes[count%nrows*ncols+3].set_xticklabels([])
  1061. axes[count%nrows*ncols+3].set_yticklabels([])
  1063. # Outlier removal and convex hull
  1064. outlier_times[count%nrows] = time()
  1065. #X = np.column_stack(np.where(img_tile==255)) # (row, col) indices
  1067. #X = np.column_stack(np.where(binary_erosion(img_tile, selem=np.ones((3,3))))) # (row, col) indices
  1069. #model = OCSVM(contamination=0.03) # > 10 min
  1070. #model = COF(contamination=0.03) # insufficient memory
  1071. #model = ABOD(contamination=0.03, n_neighbors = 5) # nonpython TypingError
  1072. #detector_list = [LOF(n_neighbors=5), LOF(n_neighbors=10), LOF(n_neighbors=15),
  1073. # LOF(n_neighbors=20), LOF(n_neighbors=25), LOF(n_neighbors=30),
  1074. # LOF(n_neighbors=35), LOF(n_neighbors=40), LOF(n_neighbors=45),
  1075. # LOF(n_neighbors=50)]
  1076. #model = LSCP(detector_list, contamination=0.03) # Combined average LOFs
  1077. #model = LOCI(contamination=0.03) # insufficient memory
  1079. #model = MCD(contamination=0.03) # 34s Good but eats edges 38s 31s
  1081. #model = PCA(contamination=0.03) # 0.02s Good but eats edges 0.04s 0.02s 0.04s
  1083. #model = KNN(contamination=0.5, n_neighbors = 9, method = 'mean', n_jobs=-1) # 1s
  1084. #model = IForest(contamination=0.03, n_jobs=-1, behaviour='new') # 3s Good
  1085. #model = CBLOF(n_clusters=9, contamination=0.03,check_estimator=False, n_jobs=-1) # 1.5-5s OK
  1086. #model = FeatureBagging(HBOS(contamination=0.03, n_bins=10, tol=0.5), contamination=0.03, n_jobs=-1) # 9-15s Very Good
  1087. #model = HBOS(contamination=0.03, n_bins=100, tol=0.5) # 0.01s Very Good and doesn't eat edges 0.1s 0.01s 0.02s
  1089. #model = LOF(n_neighbors=35, contamination=0.03, n_jobs=-1) # 1s Very Good
  1091. #model.fit(X)
  1092. #outlier_img = np.zeros_like(img_tile, dtype='uint8')
  1093. #outlier_img.flat[np.ravel_multi_index(X[np.invert(model.labels_.astype('bool'))].transpose(), outlier_img.shape)] = 255
  1094. outlier_hull_img = convex_hull_image(img_tile)
  1095. outlier_times[count%nrows] = time() - outlier_times[count%nrows]
  1097. axes[count%nrows*ncols+4].imshow(outlier_img, cmap='gray')
  1098. axes[count%nrows*ncols+4].set_title('IForest Outlier Removal')
  1099. axes[count%nrows*ncols+4].set_xticklabels([])
  1100. axes[count%nrows*ncols+4].set_yticklabels([])
  1102. axes[count%nrows*ncols+5].imshow(outlier_hull_img, cmap='gray')
  1103. axes[count%nrows*ncols+5].set_title('Outlier Hull')
  1104. axes[count%nrows*ncols+5].set_xticklabels([])
  1105. axes[count%nrows*ncols+5].set_yticklabels([])
  1107. if count%nrows == (nrows-1) or count == len(edge_labels)-1:
  1108. print('Average HBOS Runtime: {} sec'.format(np.mean(canny_times)))
  1109. print('Average IForest Runtime: {} sec'.format(np.mean(outlier_times)))
  1110. display.display(fig)
  1111. input("Press Enter to continue...")
  1112. display.clear_output(wait=True)
  1113. count += 1
  1115. print('Done')
  1117. # %%
  1118. # dev-final
  1119. import matplotlib.patches as patches
  1120. from skimage.morphology import convex_hull_image
  1122. from pyod.models.hbos import HBOS
  1123. from skimage.morphology import binary_erosion
  1125. currentDT = datetime.datetime.now(); print(str(currentDT)); startT = time()
  1126. mask_img3 = np.copy(mask_img2)
  1127. for region in regionprops(segments_slic, img):
  1128. if region.label in labels:
  1129. min_r, min_c, max_r, max_c = region.bbox
  1130. img_tile = img[min_r:max_r, min_c:max_c]
  1132. # Outlier removal and convex hull
  1133. X = np.column_stack(np.where(binary_erosion(img_tile, selem=np.ones((5,5))))) # (row, col) indices
  1135. model = HBOS(contamination=0.03, n_bins=100, tol=0.5) # 0.01s Very Good and doesn't eat edges
  1136. model.fit(X)
  1137. outlier_img = np.zeros_like(img_tile, dtype='uint8')
  1138. outlier_img.flat[np.ravel_multi_index(X[np.invert(model.labels_.astype('bool'))].transpose(), outlier_img.shape)] = 255
  1139. outlier_hull_img = convex_hull_image(outlier_img)
  1140. mask_img3[min_r:max_r, min_c:max_c] = (outlier_hull_img * 255).astype('uint8')
  1142. print("Outlier removal and convex hull: %s" % (time()-startT))
  1144. # Convert numpy array into PIL image and save to local
  1145. # Divide the WSI into regular slices to avoid PIL buffer overflow
  1146. slice_size = 30000
  1147. iters = np.uint8(np.ceil([height / slice_size, width / slice_size]))
  1149. save_img = Image.new('L', (width, height))
  1150. for row in range(iters[0]):
  1151. for col in range(iters[1]):
  1152. # Get start and end pixel location
  1153. start_r, start_c = row * slice_size, col * slice_size
  1154. end_r, end_c = min(height, start_r + slice_size), min(width, start_c + slice_size)
  1156. # Paste the slice mask to image mask
  1157. save_slice = Image.fromarray(mask_img3[start_r:end_r, start_c:end_c], 'L')
  1158. save_img.paste(save_slice, (start_c, start_r))
  1159. # Save masked image
  1160. save_img.save('NA3777-02_AB_Eosin_SLICO_0.01_iter=5_40000_finalmask.png')
  1162. print('Saved')
  1163. print("Total: %s" % (time()-startT))
  1165. # %%
  1166. # Convex hull on boundaries
  1167. from skimage.morphology import convex_hull_image
  1168. currentDT = datetime.datetime.now(); print(str(currentDT)); startT = time()
  1169. mask_img3 = np.copy(mask_img2)
  1170. for region in regionprops(segments_slic, img):
  1171. if region.label in labels:
  1172. min_r, min_c, max_r, max_c = region.bbox
  1173. img_tile = mask_img2[min_r:max_r, min_c:max_c]
  1175. # Outlier removal and convex hull
  1176. outlier_hull_img = convex_hull_image(img_tile)
  1177. mask_img3[min_r:max_r, min_c:max_c] = (outlier_hull_img * 255).astype('uint8')
  1179. print("Outlier removal and convex hull: %s" % (time()-startT))
  1181. # Convert numpy array into PIL image and save to local
  1182. # Divide the WSI into regular slices to avoid PIL buffer overflow
  1183. slice_size = 30000
  1184. iters = np.uint8(np.ceil([height / slice_size, width / slice_size]))
  1186. save_img = Image.new('L', (width, height))
  1187. for row in range(iters[0]):
  1188. for col in range(iters[1]):
  1189. # Get start and end pixel location
  1190. start_r, start_c = row * slice_size, col * slice_size
  1191. end_r, end_c = min(height, start_r + slice_size), min(width, start_c + slice_size)
  1193. # Paste the slice mask to image mask
  1194. save_slice = Image.fromarray(mask_img3[start_r:end_r, start_c:end_c], 'L')
  1195. save_img.paste(save_slice, (start_c, start_r))
  1196. # Save masked image
  1197. save_img.save('NA3777-02_AB_Eosin_SLICO_0.01_iter=5_60000_finalmask.png')
  1199. print('Saved')
  1200. print("Total: %s" % (time()-startT))
  1202. # %%
  1203. from skimage.measure import regionprops
  1204. from scipy.ndimage.morphology import binary_fill_holes
  1205. count = 0
  1206. for region in regionprops(segments_slic, img):
  1207. if count < 2000:
  1208. count+=1
  1209. continue
  1210. print(region.label)
  1211. print(region.bbox)
  1212. print(region.max_intensity)
  1213. print(region.mean_intensity)
  1214. print(region.min_intensity)
  1215. print(region.slice[0])
  1216. print(region.slice[1])
  1217. print(region.coords)
  1218. fig = plt.figure(figsize=(15,30))
  1219. ax = fig.add_subplot(1, 2, 1)
  1220. ax.imshow(img[region.bbox[0]:region.bbox[2], region.bbox[1]:region.bbox[3]], cmap='gray')
  1221. ax = fig.add_subplot(1, 2, 2)
  1222. ax.imshow(binary_fill_holes(img[region.bbox[0]:region.bbox[2], region.bbox[1]:region.bbox[3]]), cmap='gray', vmin=0, vmax=1)
  1223. break

3) Analysis - Tissue Separation.ipynb at commit d641f57, under GPL-3.0 · at the source

Overview

Authors: David Garcia1, Shivam Rajendra Rai Sharma1,2, Naomi Saito3, Laurel Beckett3, Louise Nicole C Sevilla1, La Rissa Vasquez1, Charles S DeCarli4, David Gutman5,6,7, Juan Vizcarra5,6,7, David G Coughlin8, Andrew F Teich9,10, Lorena Garcia3, Dan M Mungas4, Chen-Nee Chuah11, Brittany N Dugger1
ORCID iDs: David Garcia, Shivam Rajendra Rai Sharma, Chen-Nee Chuah, Brittany N Dugger
  1. Department of Pathology and Laboratory Medicine, University of California, Davis, Sacramento, CA, United States
  2. Department of Computer Science, University of California, Davis, Davis, CA, United States
  3. Department of Public Health Sciences, University of California, Davis, Davis, CA, United States
  4. Department of Neurology, Alzheimer’s Disease Research Center, University of California, Davis, School of Medicine, Sacramento, CA, United States
  5. Department of Neurology, Emory University School of Medicine, Atlanta, GA, United States
  6. Department of Psychiatry, Emory University School of Medicine, Atlanta, GA, United States
  7. Department of Biomedical Informatics, Emory University School of Medicine, Atlanta, GA, United States
  8. Department of Neurology, University of California, San Diego, La Jolla, CA, United States
  9. Department of Neurology, Taub Institute for Research on Alzheimer’s Disease and Aging Brain, Columbia University Medical Center, New York, NY, United States
  10. Department of Pathology and Cell Biology, Columbia University Medical Center, New York, NY, United States
  11. Department of Electrical and Computer Engineering, University of California, Davis, Davis, CA, United States
Institutions: University of California, Davis (United States); Emory University (United States); University of California San Diego (United States); Columbia University Irving Medical Center (United States); Columbia University (United States)
Journal: Journal of neuropathology and experimental neurology, volume 85, issue 6, pages 537-550
Dates: published online 10 March 2026; in print June 2026
Type: Research article · Language: English
License: CC BY-NC
Identifiers: DOI 10.1093/jnen/nlaf152 · PMID 41806384 · PMCID PMC13197125 · OpenAlex W7134954698
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), Alzheimer's / dementia (population), clinical / translational (subfield)
Methods: Connectivity, Statistics, Machine learning
Keywords: Alzheimer’s disease, amyloid-beta, convolutional neural network, digital pathology, machine learning, neurodegeneration
MeSH: Alzheimer Disease*, Amyloid beta-Peptides*, Brain*, Machine Learning*, Plaque, Amyloid*, Temporal Lobe*, Aged, Aged, 80 and over, Cerebral Amyloid Angiopathy, Female, Humans, Male (* major topic)
Topic: Dementia and Cognitive Impairment Research (Psychiatry and Mental health, Medicine), according to OpenAlex
Funding: U.S. Department of Health & Human Services | NIH | National Institute on Aging (U.S. National Institute on Aging); NIH (R01AG062517, P30AG072972, P30AG062429, P50AG008702, P30AG066462, U24NS133949, U01AG061357-S1, R01AG052132, R01AG056519); Alzheimer’s Disease Research Center (U01AG024904); National Institute of Justice (2014-R2-CX-0012)
Citations: not cited yet (Europe PMC); 76 references in the paper

Abstract

Machine learning enables scalable quantification of neuropathology, offering deeper phenotyping of Alzheimer’s disease (AD). In this validation study, we quantified amyloid-beta (Aβ) deposits, evaluating multiple brain regions across institutions, and evaluated associations with clinical, demographic, and genetic factors in persons pathologically diagnosed with AD. All linear models were adjusted for sex, age of death, ethnicity, and center. We analyzed densities (#/mm2) of cored plaques, diffuse plaques, and cerebral amyloid angiopathy (CAA) in 273 individuals from 3 Alzheimer’s Disease Research Centers. Formalin-fixed paraffin-embedded sections of frontal, temporal, and parietal cortices were immunostained and digitized, generating 799 whole-slide images (WSIs). Following log transformation, mixed-effects modeling revealed the parietal cortex had the highest cored plaque densities (P < .001); the temporal cortex had the highest diffuse plaque (P < .001); CAA showed no regional differences. Wilcoxon rank-sum test, and covariates adjusted linear models showed ApoE ε4− status was associated with higher cored plaque densities in the temporal lobe (P = .04). ApoE ε4+ status was associated with diffuse plaques in the temporal lobe (P = .001), and CAA in the frontal lobe (P = .004). These findings provide further validation and provide exploratory associations advancing deeper phenotyping of AD.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

ucdrubinet/BrainSec

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d641f5702d35be8c1ef234c704875808b9b6c0d1, 9 March 2026
Languages: Jupyter (27), Python (26), C/C++ (20), C++ (4), Shell (1), CUDA (1)
Size: 105 files, 79 scripts
Software Heritage: not archived
Found in: “DATA AVAILABILITY”
Holds: README, license file, environment (docker/Dockerfile, docker/Dockerfile.old, install/requirements.txt), tests, 27 notebooks
Not found: CITATION.cff, continuous integration, documentation
Tools: NumPy (43 files), Matplotlib (28 files), Pillow (28 files), pandas (16 files), Keras (15 files), TensorFlow (15 files), OpenCV (12 files), scikit-image (11 files), SciPy (10 files), PyTorch (9 files), scikit-learn (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
81 files
At the source: github.com/ucdrubinet/BrainSec

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 79 scripts, each with its path and the digest of its content;
  • 1 match between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

The datasets supporting the findings of this study are publicly available on the DRYAD platform. The dataset comprises over 2 TB of WSI of Aβ-immunostained human brain tissue used in this study. Due to file size constraints, the complete WSI dataset is provided in two parts: Part 1 (DOI: 10.5061/dryad.7h44j107j) and Part 2 (DOI: 10.5061/dryad.wstqjq30t). Detailed instructions for downloading the WSIs and reproducing our analysis pipeline are available on DRYAD and in our GitHub repository: https://github.com/ucdrubinet/BrainSec.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 15 authors, 6 keywords, 12 MeSH terms, 4 funders, 71 references.

Cite

This paper

Garcia, D., Rai Sharma, S. R., Saito, N., Beckett, L., Sevilla, L. N. C., Vasquez, L. R., DeCarli, C. S., Gutman, D., Vizcarra, J., Coughlin, D. G., Teich, A. F., Garcia, L., Mungas, D. M., Chuah, C.-N., & Dugger, B. N. (2026). Clinical and pathologic correlations of machine learning quantification of Aβ deposits across 3 brain regions of decedents with Alzheimer disease. Journal of neuropathology and experimental neurology, 85(6), 537-550. https://doi.org/10.1093/jnen/nlaf152

BibTeX

@article{garcia2026clinical,
author = {Garcia, David and Rai Sharma, Shivam Rajendra and Saito, Naomi and Beckett, Laurel and Sevilla, Louise Nicole C and Vasquez, La Rissa and DeCarli, Charles S and Gutman, David and Vizcarra, Juan and Coughlin, David G and Teich, Andrew F and Garcia, Lorena and Mungas, Dan M and Chuah, Chen-Nee and Dugger, Brittany N},
title = {{Clinical and pathologic correlations of machine learning quantification of Aβ deposits across 3 brain regions of decedents with Alzheimer disease}},
journal = {Journal of neuropathology and experimental neurology},
year = {2026},
month = jun,
volume = {85},
number = {6},
pages = {537--550},
publisher = {Oxford University Press},
issn = {0022-3069},
doi = {10.1093/jnen/nlaf152},
url = {https://doi.org/10.1093/jnen/nlaf152},
pmid = {41806384},
pmcid = {PMC13197125}
}

RIS

TY - JOUR
AU - Garcia, David
AU - Rai Sharma, Shivam Rajendra
AU - Saito, Naomi
AU - Beckett, Laurel
AU - Sevilla, Louise Nicole C
AU - Vasquez, La Rissa
AU - DeCarli, Charles S
AU - Gutman, David
AU - Vizcarra, Juan
AU - Coughlin, David G
AU - Teich, Andrew F
AU - Garcia, Lorena
AU - Mungas, Dan M
AU - Chuah, Chen-Nee
AU - Dugger, Brittany N
TI - Clinical and pathologic correlations of machine learning quantification of Aβ deposits across 3 brain regions of decedents with Alzheimer disease
T2 - Journal of neuropathology and experimental neurology
J2 - J Neuropathol Exp Neurol
PY - 2026
DA - 2026/06/01
VL - 85
IS - 6
SP - 537
EP - 550
SN - 0022-3069
PB - Oxford University Press
DO - 10.1093/jnen/nlaf152
UR - https://doi.org/10.1093/jnen/nlaf152
LA - en
ER -

CSL-JSON

{
"id": "10.1093/jnen/nlaf152",
"type": "article-journal",
"title": "Clinical and pathologic correlations of machine learning quantification of Aβ deposits across 3 brain regions of decedents with Alzheimer disease",
"container-title": "Journal of neuropathology and experimental neurology",
"author": [
{
"family": "Garcia",
"given": "David"
},
{
"family": "Rai Sharma",
"given": "Shivam Rajendra"
},
{
"family": "Saito",
"given": "Naomi"
},
{
"family": "Beckett",
"given": "Laurel"
},
{
"family": "Sevilla",
"given": "Louise Nicole C"
},
{
"family": "Vasquez",
"given": "La Rissa"
},
{
"family": "DeCarli",
"given": "Charles S"
},
{
"family": "Gutman",
"given": "David"
},
{
"family": "Vizcarra",
"given": "Juan"
},
{
"family": "Coughlin",
"given": "David G"
},
{
"family": "Teich",
"given": "Andrew F"
},
{
"family": "Garcia",
"given": "Lorena"
},
{
"family": "Mungas",
"given": "Dan M"
},
{
"family": "Chuah",
"given": "Chen-Nee"
},
{
"family": "Dugger",
"given": "Brittany N"
}
],
"container-title-short": "J Neuropathol Exp Neurol",
"volume": "85",
"issue": "6",
"page": "537-550",
"DOI": "10.1093/jnen/nlaf152",
"PMID": "41806384",
"PMCID": "PMC13197125",
"ISSN": "0022-3069",
"publisher": "Oxford University Press",
"URL": "https://doi.org/10.1093/jnen/nlaf152",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41598-026-61605-4 [code]
Learning precise segmentation of neurofibrillary tangles from rapid manual point annotations.
Journal: Scientific reports
In common: OpenCV, scikit-image, Pillow, 5 other tools, Alzheimer's / dementia, 9 references
[2] doi:10.1364/boe.605322 [code]
Generalized plaque digitization framework for multi-dimensional mesoscopic images.
Journal: Biomedical optics express
In common: Keras, TensorFlow, OpenCV, 8 other tools, 1 reference
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Keras, TensorFlow, OpenCV, 8 other tools
[4] doi:10.1371/journal.pcbi.1014571 [code]
SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: Keras, TensorFlow, OpenCV, 8 other tools
[5] doi:10.1126/sciadv.aed3650 [code]
Truthful visualizations for mass spectrometry imaging enable high-spatial-resolution interactive &lt;i&gt;m/z&lt;/i&gt; mapping and exploration.
Journal: Science advances
In common: Keras, TensorFlow, OpenCV, 8 other tools
[6] doi:10.1038/s41467-026-72057-9 [code]
Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.
Journal: Nature communications
In common: Keras, TensorFlow, OpenCV, 8 other tools
[7] doi:10.1002/epi.70296 [code]
Fully automated three-dimensional deep learning-based magnetic resonance imaging segmentation of brain cavities in epilepsy surgery.
Journal: Epilepsia
In common: Keras, TensorFlow, OpenCV, 7 other tools, clinical / translational
[8] doi:10.1186/s12880-026-02481-2 [code]
Deep learning-based neuroanatomical profiling reveals population-specific brain changes in multiple sclerosis: a large-scale Middle Eastern study.
Journal: BMC medical imaging
In common: Keras, TensorFlow, OpenCV, 7 other tools, clinical / translational
[9] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Keras, TensorFlow, OpenCV, 7 other tools
[10] doi:10.1016/j.patter.2026.101538 [code]
A multi-modal foundation model for brain disease diagnosis and medical imaging.
Journal: Patterns (New York, N.Y.)
In common: TensorFlow, OpenCV, scikit-image, 7 other tools, clinical / translational

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.