<jupyter_start><jupyter_text>3 - Model loading and NM size/yield analysisIn this notebook we will:1. Load an image for analysis.2. Load our previously-trained model.3. Use model to label the SEM image.4. Perform post-processing on model output to learn about our NM characteristics Note: A GPU instance is not necessary for this notebook as we will only be performing inference which is not as computationally-expensive as training.  Install detectron2Again, we will be using Facebook's [detectron2](https://github.com/facebookresearch/detectron2) library to run the interence on our images to let's install it.<jupyter_code># install dependencies: (use cu101 because colab has CUDA 10.1)
!pip install -U torch==1.5 torchvision==0.6 -f https://download.pytorch.org/whl/cu101/torch_stable.html 
!pip install cython pyyaml==5.1
!pip install -U 'git+https://github.com/cocodataset/cocoapi.git#subdirectory=PythonAPI'
import torch, torchvision
print(torch.__version__, torch.cuda.is_available())
!gcc --version
# install detectron2:
!pip install detectron2==0.1.2 -f https://dl.fbaipublicfiles.com/detectron2/wheels/cu101/index.html
# Setup detectron2 logger
import detectron2
from detectron2.utils.logger import setup_logger
setup_logger('logs')

# import some common libraries
import numpy as np
import pandas as pd
import seaborn as sns
import os, cv2, random, tifffile, json, datetime, time, urllib
from glob import glob
from google.colab.patches import cv2_imshow
from PIL import Image
from pathlib import Path
from matplotlib import pyplot as plt
plt.rc('axes', axisbelow=True)

# import some common detectron2 utilities
from detectron2 import model_zoo
from detectron2.engine import DefaultPredictor
from detectron2.config import get_cfg
from detectron2.utils.visualizer import Visualizer
from detectron2.data import DatasetCatalog, MetadataCatalog
# Define the classes in our dataset
class_dict = {'slit': '1',
              'nanomembrane': '2',
              'parasitic': '3',
              'bottom_nucleus': '4',
              'side_nucleus': '5',
              'nanowire': '6',
              'overgrowth': '7'}

# Define a dict that maps class number back to class name too
inv_class_dict = dict(map(reversed, class_dict.items()))

# Define some paths/constants that will be useful later
desired_mag = 50000  # Used to filter the TIFF input files

root = Path('./DeepSEM/segmentation-NMs/')
dataset_dir = root.joinpath('datasets')
output_dir = root.joinpath('output')
models_dir = root.joinpath('trained_models')

github_url = 'https://github.com/Martin09/' + str(root).replace('/','/trunk/')

imgs_zip = dataset_dir.joinpath('Nick_NMs_allrawimgs.zip')
imgs_dir = dataset_dir.joinpath(imgs_zip.stem)
imgs_google_drive_id = '1M2_0GLScsNY53w8hU2xJdXtisESfkOqI'

test_dir = imgs_dir.joinpath('test')
train_dir = imgs_dir.joinpath('train')

dataset_root_name = 'nm_masks'
train_name = dataset_root_name + '_train'
test_name = dataset_root_name + '_test'

model_path = models_dir.joinpath('nm_seg_it20k_loss0.028.yaml')
weights_path = models_dir.joinpath('nm_seg_it20k_loss0.028.pth')

weights_google_drive_id = '1btMy-EyU2sTSSPQk8kf663DYgO-sD3kR'<jupyter_output><empty_output><jupyter_text>3.1 - Unpack and load our images<jupyter_code># # Optional: Save everything to your own GoogleDrive
# from google.colab import drive
# drive.mount('/content/gdrive/')
# %cd "/content/gdrive/My Drive/path/to/save/location"

# Clone just the relevant folder from the DeepSEM repo
!rm -rf $root
!apt install subversion
!svn checkout $github_url $root

# # Alternative: Clone whole DeepSEM repository
# !rm -rf DeepSEM  # Remove folder if it already exists
# !git clone https://github.com/Martin09/DeepSEM<jupyter_output><empty_output><jupyter_text>For simplicity, I will use our previous training images for inference. However these could be replaced with any similar un-labelled SEM images.<jupyter_code># Check if .zip file exists, if not, download it from Google Drive
if imgs_zip.exists():
  print('Dataset already exists. Skipping download!')
else:
  print('Dataset does not exist... Downloading!')
  !gdown --id $imgs_google_drive_id -O $imgs_zip

# Unzip raw dataset
!rm -rf $imgs_dir
!unzip -o $imgs_zip -d $imgs_dir<jupyter_output><empty_output><jupyter_text>Now we will sort the input files which have many different magnifications into images that only have the desired magnification (50k in this case).<jupyter_code>in_files = list(imgs_dir.rglob('*.tif'))

images = []
# Start to loop over all TIFF files
for file in in_files:
    # Open each file using the TiffFile library
    with tifffile.TiffFile(file) as tif:
        
        # Extract magnification data
        mag = tif.sem_metadata['ap_mag'][1] 
        if type(mag) is str:  # Apply correction for "k" ex: mag = "50 k"
            mag = float(mag.split(' ')[0]) * 1000
        else:
            mag = float(mag)

        # Only filter the images that have the magnification that we are interested in
        if not mag == desired_mag:
          continue

    images.append(file)<jupyter_output><empty_output><jupyter_text>Load a random image and show it.<jupyter_code>im_path = random.sample(images,1)[0]
im = cv2.imread(str(im_path), cv2.IMREAD_GRAYSCALE)
print(im.shape)
cv2_imshow(im)<jupyter_output><empty_output><jupyter_text>Do some pre-processing to get it ready to feed into our model.<jupyter_code># Model expects an RGB image, so copy the greyscale data into other 2 channels
im_RGB = np.repeat(im[:, :, np.newaxis], 3, axis=2)
print(im_RGB.shape)
cv2_imshow(im_RGB)

# For use later
img_h = im_RGB.shape[0]
img_w = im_RGB.shape[1]<jupyter_output><empty_output><jupyter_text>3.2 - Load our modelNow we will load a trained model and use it to label the above image. First we load a default config with `get_cfg()` and we then overwrite some of its parameters with our saved YAML configuration file. One important point is that we need to have `cfg.MODEL.WEIGHTS` set to point to the weights file. As this file can be quite big (>300MB) and since Github isn't designed to host big binary files, I have saved the weights for this model on my Google Drive instead. However, if you have your weights saved locally (ex: on your Google Drive), you can skip this download.<jupyter_code># Check if .zip file exists, if not, download it from Google Drive
if weights_path.exists():
  print('Dataset already exists. Skipping download!')
else:
  print('Dataset does not exist... Downloading!')
  !gdown --id $weights_google_drive_id -O $weights_path<jupyter_output><empty_output><jupyter_text>Now we can go ahead with the rest of the configuration of the model.<jupyter_code>cfg = get_cfg()
cfg.merge_from_file(model_path)
cfg.MODEL.WEIGHTS = str(weights_path)
cfg.MODEL.DEVICE = 'cpu'  # CPU is enough for inference, no need for GPU

# If we have a lot of objects to detect, need to set higher # of proposals here:
cfg.MODEL.RPN.POST_NMS_TOPK_TEST = 1000
cfg.MODEL.RPN.PRE_NMS_TOPK_TEST = 1000
cfg.TEST.DETECTIONS_PER_IMAGE = 200

cfg.MODEL.ROI_HEADS.SCORE_THRESH_TEST = 0.5   # Set the testing threshold for this model
cfg.MODEL.ROI_HEADS.NMS_THRESH_TEST = 0.2     # Non-max supression threshold
cfg.MODEL.ROI_HEADS.NUM_CLASSES = len(class_dict) # We have three classification classes 

# Setting allowed input sizes (avoid scaling)
cfg.INPUT.MIN_SIZE_TEST = 0
cfg.INPUT.MAX_SIZE_TEST = 99999


# A bit of a hacky way to be able to use the DefaultPredictor:
# Register a "fake" dataset to then set the 'thing_classes' metadata
# (there is probably a better way to do this...)
cfg.DATASETS.TEST = ('placeholder')
DatasetCatalog.clear()
DatasetCatalog.register("placeholder", lambda _: None)
MetadataCatalog.get("placeholder").set(thing_classes=list(class_dict))
predictor = DefaultPredictor(cfg)
outputs = predictor(im_RGB)
print('Number of detected objects = {}'.format(len(outputs["instances"])))
# Verify outputs manually
# outputs["instances"].pred_classes
# outputs["instances"].pred_boxes
# outputs["instances"].scores
# We can use Visualizer to draw the predictions on the image.
v = Visualizer(im_RGB[:, :, ::-1], MetadataCatalog.get(cfg.DATASETS.TEST[0]), scale=1.5)
v = v.draw_instance_predictions(outputs["instances"].to("cpu"))
cv2_imshow(v.get_image()[:, :, ::-1])<jupyter_output><empty_output><jupyter_text>3.4 - Post-processing model outputHowever, just getting the output from the model isn't enough. Now we have to do bit more work to post-process the output and extract things like nanomembrane yield, sizes and other interesting data!First lets divide up the output of the neural net for further processing:<jupyter_code>cl = np.array(outputs["instances"].pred_classes.cpu())  # Classes
s = np.array(outputs["instances"].scores.cpu()) # Prediction scores
b =  np.array([x.numpy() for x in outputs["instances"].pred_boxes])  # Bounding boxes
c = np.array(outputs["instances"].pred_boxes.get_centers())  # Bounding box centres
m =  np.array([x.numpy() for x in outputs["instances"].pred_masks])  # Segmentation masks<jupyter_output><empty_output><jupyter_text>Now we can loop over all the possible classes and display images with segmentation masks of each class individually.<jupyter_code>for clas in range(len(class_dict)):
  i_filt = list(np.argwhere(cl==clas).flatten()) # Choose only the indixes with specific class

  print(f"{inv_class_dict[str(clas+1)]}:")

  # We can use Visualizer to draw the predictions on the image.
  v = Visualizer(im_RGB[:, :, ::-1], MetadataCatalog.get(cfg.DATASETS.TEST[0]), scale=1.0)
  v = v.draw_instance_predictions(outputs["instances"][[i_filt]].to("cpu"))
  cv2_imshow(v.get_image()[:, :, ::-1])<jupyter_output><empty_output><jupyter_text>Now, before we can start to mess around with dimensional analysis we first need to extract the pixel size from the raw TIF image:<jupyter_code>with tifffile.TiffFile(im_path) as tif:
    
    # Extract magnification data
    mag = tif.sem_metadata['ap_mag'][1] 
    if type(mag) is str:  # Apply correction for "k" ex: mag = "50 k"
        mag = float(mag.split(' ')[0]) * 1000
    else:
        mag = float(mag)

    # Extract pixel size data
    pixel_size = float(tif.sem_metadata['ap_pixel_size'][1])  # nm
    if 'µm' in tif.sem_metadata['ap_pixel_size'][2]: # Correction for um
        pixel_size *= 1000

    # Extract tilt data
    tilt = tif.sem_metadata['ap_stage_at_t'][1]  # degrees tilt

pixel_size_x = pixel_size  # nm
pixel_size_y = pixel_size / np.cos(np.deg2rad(tilt))  # nm<jupyter_output><empty_output><jupyter_text>Let's put the output into a handy Pandas Dataframe before any more processing:<jupyter_code># Define data structure
data = { 'class':[],
         'class_num':[],
         'score':[],
         'bbox':[],
         'bbox_centre':[],
         'height':[],
         'width':[],
         'mask':[],
         'area':[],
         'area_bbox':[]}

# Loop over all objects
for i in range(len(outputs["instances"])):
  
  data['class'].append(inv_class_dict[str(cl[i]+1)])
  data['class_num'].append(cl[i])

  data['score'].append(s[i])

  data['bbox'].append(b[i])
  data['bbox_centre'].append(c[i])

  h = (b[i,3] - b[i,1]) * pixel_size_y
  data['height'].append(h)

  w = (b[i,2] - b[i,0]) * pixel_size_x  
  data['width'].append(w)

  data['area_bbox'].append(w*h)

  data['mask'].append(m[i])

  data['area'].append(m[i].astype(int).sum() * pixel_size_x * pixel_size_y)

df = pd.DataFrame.from_dict(data)<jupyter_output><empty_output><jupyter_text>Let's plot a simple bar graph of the number of objects in each class.<jupyter_code>fig_size = (8,6)

fig = plt.figure(figsize=fig_size)

counts = df['class'].value_counts()
total = sum(counts)
values = [c/total*100 for c in counts]
labels = list(counts.index)

sns.barplot(x=labels, y=values) # height should be three times width);

for i, v in enumerate(values):
    plt.text(i, v + np.max(values)*0.01, f'{v:.1f}% ({counts[i]:.0f})', color='k', ha='center')

plt.ylim([0, np.max(values)*(1.1)])
plt.title('Growth Structure Yield (count)')
plt.ylabel('Yield (%)')
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Let's start with slit length/width.<jupyter_code>df_slits = df[df['class']=='slit']

print(f"Mean slit width: {df_slits['height'].mean():.0f} +/- {df_slits['height'].std():.0f} nm")
print(f"Mean slit length: {df_slits['width'].mean():.0f} +/- {df_slits['width'].std():.0f} nm")<jupyter_output><empty_output><jupyter_text>Can use the slit bounding boxes to plot a comparison of total area fraction:<jupyter_code># Calculate slit area using the bounding box height (slit width) and image width
# -> Avoids problems of weirdly-shaped slit segmentation masks
slit_area = 0
for h in df_slits['height']:
  slit_area += h * (img_w * pixel_size_x)

print(f"slit area (bbox): {slit_area:.0f} nm2")

class_areas = []

# Loop over all classes
for clas in range(len(class_dict)):
  
  df_filt = df[df['class_num']==clas]

  # There isn't any object with this class in the current image
  if df_filt.size == 0:
      print(f"{inv_class_dict[str(clas+1)]} area: 0 nm2")
      class_areas.append(0)
      continue

  # Stack all masks in this class together
  overall_mask = df_filt['mask'].sum()

  # This approach avoids double counting pixels if masks overlap with eachother
  overall_mask = overall_mask.astype(bool).astype(int)

  # Calculate area
  area = overall_mask.sum() * pixel_size_x * pixel_size_y
  print(f"{inv_class_dict[str(clas+1)]} area: {area:.0f} nm2")
  
  # Add area to the areas list
  class_areas.append(area)

# Remove the "slit" class entry of the areas array since we will use the bbox value from above
class_areas = class_areas[1:]<jupyter_output><empty_output><jupyter_text>Now we can plot the total area of each class:<jupyter_code>fig = plt.figure(figsize=fig_size)

values = [area/slit_area*100 for area in class_areas]
labels = [inv_class_dict[str(clas+1)] for clas in range(1, len(class_dict))]

sns.barplot(x=labels, y=values)

for i, v in enumerate(values):
    plt.text(i, v + np.max(values)*0.01, f'{v:.1f}%\n({class_areas[i]:.0f} nm$^2$)', color='k', ha='center')

plt.ylim([0, np.max(values)*(1.1)])
plt.title('Growth Structure Yield (% of slit area)')
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Then we can also look at nanomembrane dimensions:<jupyter_code>df_nms = df[df['class']=='nanomembrane']
print(f"Mean nanomembrane width: {df_nms['height'].mean():.0f} +/- {df_nms['height'].std():.0f} nm")
print(f"Mean nanomembrane length: {df_nms['width'].mean():.0f} +/- {df_nms['width'].std():.0f} nm")
print(f"Max nanomembrane length: {df_nms['width'].max():.0f} nm")<jupyter_output><empty_output><jupyter_text>If we want to compare the dimensions of all classes we can define a handy plotting function that we then use to plot different values:<jupyter_code># Define a helper function for easy box plotting
def make_box_plot(dat, x, x_lab):
  # Create figure
  fig = plt.figure(figsize=fig_size)

  # Plot the orbital period with horizontal boxes
  ax = sns.boxplot(x=x, y='class', data=dat, whis=[0, 100], palette="vlag")
  ax.set(xlim=(0, dat[x].max()*1.2))
  ax.xaxis.grid(True)

  # Add in points to show each observation
  g = sns.swarmplot(x=x, y='class', data=dat, size=5, color=".3", linewidth=0)

  # Tweak the visual presentation
  plt.xlabel(x_lab)
  plt.ylabel("")
  title = x_lab.split(" ")[0]
  plt.title(f"{title} Distributions")
# Filter dataframe to exclude the slits, plot the data
df_no_slit = df[df['class']!='slit']
make_box_plot(df_no_slit,'height','Width (nm)');
# Filter dataframe to exclude the slits, plot the data
df_no_slit = df[df['class']!='slit']
make_box_plot(df_no_slit,'width','Length (nm)')
# Filter dataframe to exclude the slits, plot the data
df_no_slit = df[df['class']!='slit']
make_box_plot(df_no_slit,'area','Area (nm$^2$)')
# Plot a scatter plot of growth structure length vs width to try and see some trends
fig = plt.figure(figsize=fig_size)
ax = sns.scatterplot(x="width", y="height", hue="class", style="class", s=50, data=df_no_slit)
ax.grid()
ax.set(xlabel='Length (nm)')
ax.set(ylabel='Width (nm)')
ax.set(title='Length vs Width Plot')

# Add a horizontal line to represent median slit width
ax.axhline(df_slits['height'].median(), ls='--', color='r')
ax.text(0, df_slits['height'].median()*1.02, "Median Slit Width", color='r');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1296863513844646, 'fraction_numerical': 0.014457240872335212, 'mean_word_length': 5.028434268833087, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 75, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '211805', 'license': 'MIT', 'n_tokens_mistral': 5843, 'n_tokens_neox': 5614, 'n_words': 1560, 'path': 'segmentation-NMs/3_nm_seg_inference.ipynb', 'repo_name': 'Martin09/DeepSEM', 'size': '16324'}