<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Nanocrystal Segmentation
This notebook demonstrates two approaches to nanocrystal segmentation
Step1:  Load demonstration data
Step2:  Plot data to inspect
Step3:  Remove the background
Step4:  Plot the background subtracted data
Step5:  Find the position of the direct beam in the background subtracted data.
Step6:  Apply the same shifts to the raw data.
Step7:  Set calibrations
Step8:  <a id='vdf'></a>
2. Virtual Image Based Segmentation
2.1. Peak Finding & Refinement
Find all diffraction peaks for all PED patterns. 
The parameters were found by interactive peak finding
Step9:  Visualise the number of diffraction peaks found at each probe position
Step10:  Exclude peaks too close to the detector edge for sub-pixel refinement.
Step11:  Refine the peak positions using center of mass
Step12:  2.2. Determine Unique Peaks
Determine the unique diffraction peaks by clustering
Step13:  Visualise the detected unique peaks by plotting them on the maximum of the signal.
Step14:  Visualise both the clusters and the unique peaks obtained after DBSCAN clustering. 
NB The cluster colors are randomly generated, so run it again if it is hard to discern two close clusters.
Step15:  Filter the unique vectors by magnitude in order to exclude the direct beam from the following analysis
Step16:  Plot the unique vectors
Step17:  Optionally save and load the unique peaks
Step18:  2.3. Virtual Imaging & Segmentation
Calculate VDF images for all unique peaks
Step19:  Plot the VDF images for inspection
Step20:  First find adequate parameters by looking at watershed segmentation of a single VDF image.
Step21:  Perform segmentation on all the VDF images
Step22:  Plot the segments for inspection
Step23:  Calculate normalised cross-correlations between all VDF image segments to identify those that are related to the same crystal.
Step24:  If the correlation value exceeds corr_threshold for certain segments, those segments are summed. These segments are discarded if the number of these segments are below vector_threshold, as this number corresponds to the number of detected diffraction peaks associated with the single crystal. The vector_threshold criteria is included to avoid including segment images resulting from noise or incorrect segmentation.
Step25:  Simulate virtual diffraction patterns for each summed segment
Step26:  Plot the final results from the VDF image-based segmentation
Step27:  <a id='nmf'></a> 
3. NMF Based Segmentation
3.1. NMF Decomposition
Create a signal mask so that the region in the centre of each PED pattern, including the direct beam, can be excluded in the machine learning.
Step28:  Perform single value decomposition (SVD)
Step29:  Investigate the scree plot and use it as a guide to determine the number of components
Step30:  Perform NMF decomposition with specified number of components
Step31:  Plot the NMF results
Step32:  Discard the components related to background (#0) and to the carbon film (#4)
Step33:  3.2. Correlate NMF Loading Maps
NMF often leads to splitting of some crystals into several components. Therefore the correlation between loadings and between component patterns are calculated, and if both the correlation values for loadings and factors exceed threshold values, those loadings and factors are summed. 
Calculate the matrix of normalised cross-correlation for both the loadings and patterns first, to find suitable correlation threshold values.
Step34:  Perform correlation and summation of the factors and loadings
Step35:  Plot the NMF reuslts after correlation and summation
Step36:  First investigate how the parameters influence the segmentation on
one single loading map.
Step37:  Set a threshold for the minimum intensity value that a loading segment must contain in order to be kept.
Step38:  Plot the final results from the NMF-based segmentation
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import hyperspy.api as hs
import matplotlib.pyplot as plt
import pyxem as pxm
Explanation: Nanocrystal Segmentation
This notebook demonstrates two approaches to nanocrystal segmentation:
1. Virtual dark-field (VDF) imaging-based segmentation
2. Non-negative matrix factorisation (NMF)-based segmentation
The segmentation is demonstrated on a SPED dataset of partly overlapping MgO nanoparticles, where some of the particles share the same orientation.
This functionality has been checked to run with pyxem-0.11.0 (May 2020). Bugs are always possible, do not trust the code blindly, and if you experience any issues please report them here: https://github.com/pyxem/pyxem-demos/issues
Contents
<a href='#gen'> Setting up, Loading Data, Pre-processing</a>
<a href='#vdf'> Virtual Image Based Segmentation</a>
<a href='#nmf'> NMF Based Segmentation</a>
<a id='gen'></a>
1. Setting up, Loading Data, Pre-processing
Import pyxem and other required libraries
End of explanation
dp = hs.load('./data/06/mgo_nanoparticles.hdf5')
dp
Explanation: Load demonstration data
End of explanation
dp.plot(cmap='magma_r')
Explanation: Plot data to inspect
End of explanation
sigma_min = 1.7
sigma_max = 13.2
dp_rb= dp.subtract_diffraction_background('difference of gaussians', 
                             min_sigma=sigma_min, 
                             max_sigma=sigma_max)
dp.data
dp_rb.data.compute()
Explanation: Remove the background
End of explanation
dp_rb.plot(cmap='magma_r')
Explanation: Plot the background subtracted data
End of explanation
shifts = dp_rb.center_direct_beam(method='cross_correlate',
                                  half_square_width=15,
                                  return_shifts=True,
                                  radius_start=2,
                                  radius_finish=6)
# this step is likely to be quiet slow
shifts.compute()
Explanation: Find the position of the direct beam in the background subtracted data.
End of explanation
dp_rb.align2D(shifts=shifts, crop=False)
Explanation: Apply the same shifts to the raw data.
End of explanation
scale = 0.03246*2 
scale_real = 3.03*2
dp.set_diffraction_calibration(scale)
dp.set_scan_calibration(scale_real)
# see https://github.com/pyxem/pyxem/issues/717 for why have to do this
dp_rb.compute()
dp_rb = pxm.signals.ElectronDiffraction2D(dp_rb)
dp_rb.set_diffraction_calibration(scale)
dp_rb.set_scan_calibration(scale_real)
Explanation: Set calibrations
End of explanation
import hyperspy
hyperspy.__version__
peaks = dp_rb.find_peaks(method='minmax', 
                         interactive=False)
Explanation: <a id='vdf'></a>
2. Virtual Image Based Segmentation
2.1. Peak Finding & Refinement
Find all diffraction peaks for all PED patterns. 
The parameters were found by interactive peak finding:
peaks = dp_rb.find_peaks_interactive(imshow_kwargs={'cmap': 'magma_r'})
End of explanation
from pyxem.signals.diffraction_vectors import DiffractionVectors
peaks = DiffractionVectors.from_peaks(peaks,center=(36,36),calibration=scale)
diff_map = peaks.get_diffracting_pixels_map()
diff_map.plot()
Explanation: Visualise the number of diffraction peaks found at each probe position
End of explanation
# Currently a Bug.  Need to Fix
peaks.detector_shape= (72,72)
peaks.pixel_calibration=scale
peaks_filtered = peaks.filter_detector_edge(exclude_width=2)
Explanation: Exclude peaks too close to the detector edge for sub-pixel refinement.
End of explanation
from pyxem.generators.subpixelrefinement_generator import SubpixelrefinementGenerator
from pyxem.signals.diffraction_vectors import DiffractionVectors
refine_gen = SubpixelrefinementGenerator(dp_rb, peaks_filtered)
peaks_refined = DiffractionVectors(refine_gen.center_of_mass_method(square_size=4))
peaks_refined.axes_manager.set_signal_dimension(0)
Explanation: Refine the peak positions using center of mass
End of explanation
def replace_nan(data):
    new_data = data[0]
    new_data[np.isnan(data[0])]=0
    return np.array(new_data)
peaks_refined = peaks_refined.map(replace_nan, inplace=False, ragged=True)
peaks_refined.set_signal_type("diffraction_vectors")
distance_threshold = scale*0.89
min_samples = 10
unique_peaks = peaks_refined.get_unique_vectors(method='DBSCAN',
                                                distance_threshold=distance_threshold,
                                                min_samples=min_samples)
print(np.shape(unique_peaks.data)[0], ' unique vectors were found.')
Explanation: 2.2. Determine Unique Peaks
Determine the unique diffraction peaks by clustering
End of explanation
radius_px = dp_rb.axes_manager.signal_shape[0]/2
reciprocal_radius = radius_px * scale
unique_peaks.plot_diffraction_vectors(
    method='DBSCAN',
    unique_vectors=unique_peaks,
    distance_threshold=distance_threshold,
    xlim=reciprocal_radius,
    ylim=reciprocal_radius,
    min_samples=min_samples,
    image_to_plot_on=dp_rb.max(),
    image_cmap='magma_r',
    plot_label_colors=False)
Explanation: Visualise the detected unique peaks by plotting them on the maximum of the signal.
End of explanation
peaks_refined.plot_diffraction_vectors(
    method='DBSCAN',
    xlim=reciprocal_radius, 
    ylim=reciprocal_radius,
    unique_vectors=unique_peaks, 
    distance_threshold=distance_threshold,
    min_samples=min_samples, 
    image_to_plot_on=dp_rb.max(), 
    image_cmap='gray_r',
    plot_label_colors=True, 
    distance_threshold_all=scale*0.1)
Explanation: Visualise both the clusters and the unique peaks obtained after DBSCAN clustering. 
NB The cluster colors are randomly generated, so run it again if it is hard to discern two close clusters.
End of explanation
Gs = unique_peaks.filter_magnitude(min_magnitude=10*scale,
                                   max_magnitude=np.inf)
print(np.shape(Gs)[0], ' unique vectors.')
Explanation: Filter the unique vectors by magnitude in order to exclude the direct beam from the following analysis
End of explanation
Gs.plot_diffraction_vectors(unique_vectors=Gs,
                            distance_threshold=distance_threshold,
                            xlim=reciprocal_radius,
                            ylim=reciprocal_radius,
                            min_samples=min_samples,
                            image_to_plot_on=dp_rb.max(),
                            image_cmap='magma',
                            plot_label_colors=False)
Explanation: Plot the unique vectors
End of explanation
np.save('peaks.npy', Gs.data)
Gs = np.load('peaks.npy', allow_pickle=True)
Gs = pxm.signals.DiffractionVectors(Gs)
Gs.axes_manager.set_signal_dimension(0)
Explanation: Optionally save and load the unique peaks
End of explanation
from pyxem.generators import VirtualDarkFieldGenerator
radius=scale*2
vdfgen = VirtualDarkFieldGenerator(dp_rb, Gs)
VDFs = vdfgen.get_virtual_dark_field_images(radius=radius)
Explanation: 2.3. Virtual Imaging & Segmentation
Calculate VDF images for all unique peaks
End of explanation
#%matplotlib notebook
VDFs.plot(cmap='magma', scalebar=False)
Explanation: Plot the VDF images for inspection
End of explanation
from pyxem.utils.segment_utils import separate_watershed
min_distance = 5.5
min_size = 10
max_size = 1000
max_number_of_grains = 1000
marker_radius = 2
exclude_border = 2
i = 25
sep_i = separate_watershed(
    VDFs.inav[i].data, min_distance=min_distance, min_size=min_size,
    max_size=max_size, max_number_of_grains=max_number_of_grains,
    exclude_border=exclude_border, marker_radius=marker_radius,
    threshold=True, plot_on=True)
Explanation: First find adequate parameters by looking at watershed segmentation of a single VDF image.
End of explanation
segs = VDFs.get_vdf_segments(min_distance=min_distance,
                             min_size=min_size,
                             max_size = max_size,
                             max_number_of_grains = max_number_of_grains,
                             exclude_border=exclude_border,
                             marker_radius=marker_radius,
                             threshold=True)
print(np.shape(segs.segments)[0],' segments were found.')
Explanation: Perform segmentation on all the VDF images
End of explanation
segs.segments.plot(cmap='magma_r')
Explanation: Plot the segments for inspection
End of explanation
ncc_vdf = segs.get_ncc_matrix()
ncc_vdf.plot(scalebar=False, cmap='RdBu')
Explanation: Calculate normalised cross-correlations between all VDF image segments to identify those that are related to the same crystal.
End of explanation
corr_threshold=0.7
vector_threshold=5
segment_threshold=4
corrsegs = segs.correlate_vdf_segments(
    corr_threshold=corr_threshold, vector_threshold=vector_threshold,
    segment_threshold=segment_threshold)
print(np.shape(corrsegs.segments)[0],' correlated segments were found.')
Explanation: If the correlation value exceeds corr_threshold for certain segments, those segments are summed. These segments are discarded if the number of these segments are below vector_threshold, as this number corresponds to the number of detected diffraction peaks associated with the single crystal. The vector_threshold criteria is included to avoid including segment images resulting from noise or incorrect segmentation.
End of explanation
sigma = scale*1.5
virtual_sig = corrsegs.get_virtual_electron_diffraction(
    calibration=scale, shape=(int(radius_px*2), int(radius_px*2)), sigma=sigma)
virtual_sig.set_diffraction_calibration(scale)
Explanation: Simulate virtual diffraction patterns for each summed segment
End of explanation
hs.plot.plot_images(corrsegs.segments, cmap='magma_r', axes_decor='off',
                    per_row=np.shape(corrsegs.segments)[0],
                    suptitle='', scalebar=False, scalebar_color='white',
                    colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right':0.78})
hs.plot.plot_images(virtual_sig, cmap='magma_r', axes_decor='off',
                    per_row=np.shape(corrsegs.segments)[0],
                    suptitle='', scalebar=False, scalebar_color='white',
                    colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right': 0.78})
Explanation: Plot the final results from the VDF image-based segmentation
End of explanation
dpm = pxm.signals.Diffraction2D(dp.inav[0,0])
signal_mask = dpm.get_direct_beam_mask(radius=10)
signal_mask.plot()
Explanation: <a id='nmf'></a> 
3. NMF Based Segmentation
3.1. NMF Decomposition
Create a signal mask so that the region in the centre of each PED pattern, including the direct beam, can be excluded in the machine learning.
End of explanation
dp.change_dtype('float32')
dp.decomposition(algorithm='svd',
                 normalize_poissonian_noise=True,
                 centre=None,
                 signal_mask=signal_mask.data)
dp.plot_decomposition_results()
Explanation: Perform single value decomposition (SVD)
End of explanation
num_comp=11
ax = dp.plot_explained_variance_ratio(n=200, threshold=num_comp,
                                      hline=True, xaxis_labeling='ordinal',
                                      signal_fmt={'color':'k', 'marker':'.'}, 
                                      noise_fmt={'color':'gray', 'marker':'.'})
Explanation: Investigate the scree plot and use it as a guide to determine the number of components
End of explanation
dp.decomposition(normalize_poissonian_noise=True,
                 algorithm='nmf',
                 output_dimension=num_comp,
                 signal_mask=signal_mask.data)
dp_nmf = dp.get_decomposition_model(components=np.arange(num_comp))
factors = dp_nmf.get_decomposition_factors()
loadings = dp_nmf.get_decomposition_loadings()
Explanation: Perform NMF decomposition with specified number of components
End of explanation
hs.plot.plot_images(loadings, cmap='magma_r', axes_decor='off', per_row=11,
             suptitle='', scalebar=False, scalebar_color='white', colorbar=False,
             padding={'top': 0.95, 'bottom': 0.05,
                      'left': 0.05, 'right':0.78})
hs.plot.plot_images(factors, cmap='magma_r', axes_decor='off', per_row=11,
             suptitle='', scalebar=False, scalebar_color='white', colorbar=False,
             padding={'top': 0.95, 'bottom': 0.05,
                      'left': 0.05, 'right':0.78})
Explanation: Plot the NMF results
End of explanation
from hyperspy.signals import Signal2D
factors = Signal2D(np.delete(factors.data, [0, 4], axis = 0))
loadings = Signal2D(np.delete(loadings.data, [0, 4], axis = 0))
hs.plot.plot_images(factors, cmap='magma_r', axes_decor='off',
                    per_row=9, suptitle='', scalebar=False,
                    scalebar_color='white', colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right':0.78})
hs.plot.plot_images(loadings, cmap='magma_r', axes_decor='off',
                    per_row=9, suptitle='', scalebar=False,
                    scalebar_color='white', colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right':0.78})
Explanation: Discard the components related to background (#0) and to the carbon film (#4)
End of explanation
from pyxem.signals.segments import LearningSegment
learn = LearningSegment(factors=factors, loadings=loadings)
ncc_nmf = learn.get_ncc_matrix()
ncc_nmf.plot(scalebar=False, cmap='RdBu')
corr_th_factors = 0.45
corr_th_loadings = 0.3
Explanation: 3.2. Correlate NMF Loading Maps
NMF often leads to splitting of some crystals into several components. Therefore the correlation between loadings and between component patterns are calculated, and if both the correlation values for loadings and factors exceed threshold values, those loadings and factors are summed. 
Calculate the matrix of normalised cross-correlation for both the loadings and patterns first, to find suitable correlation threshold values.
End of explanation
learn_corr = learn.correlate_learning_segments(corr_th_factors=corr_th_factors,
                                               corr_th_loadings=corr_th_loadings)
Explanation: Perform correlation and summation of the factors and loadings
End of explanation
hs.plot.plot_images(learn_corr.loadings, cmap='magma_r', axes_decor='off',
                    per_row=7, suptitle='', scalebar=False,
                    scalebar_color='white', colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right':0.78})
hs.plot.plot_images(learn_corr.factors, cmap='magma_r', axes_decor='off',
                    per_row=7, suptitle='', scalebar=False,
                    scalebar_color='white', colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right':0.78})
Explanation: Plot the NMF reuslts after correlation and summation
End of explanation
from pyxem.utils.segment_utils import separate_watershed
min_distance = 10
min_size = 50
max_size = 100000
max_number_of_grains = 100000
marker_radius = 2
exclude_border = 1
threshold = True
i =1
sep_i = separate_watershed(
    learn_corr.loadings.data[i], min_distance=min_distance,
    min_size=min_size, max_size=max_size, 
    max_number_of_grains=max_number_of_grains,
    exclude_border=exclude_border, 
    marker_radius=marker_radius, threshold=True, plot_on=True)
Explanation: First investigate how the parameters influence the segmentation on
one single loading map.
End of explanation
min_intensity_threshold = 10000
learn_corr_seg = learn_corr.separate_learning_segments(
    min_intensity_threshold=min_intensity_threshold,
    min_distance = min_distance, min_size = min_size,
    max_size = max_size, 
    max_number_of_grains = max_number_of_grains,
    exclude_border = exclude_border,
    marker_radius = marker_radius, threshold = True)
Explanation: Set a threshold for the minimum intensity value that a loading segment must contain in order to be kept.
End of explanation
hs.plot.plot_images(learn_corr_seg.loadings, 
                    cmap='magma_r', axes_decor='off',
                    per_row=10, suptitle='', scalebar=False,
                    scalebar_color='white', colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right':0.78})
hs.plot.plot_images(learn_corr_seg.factors, 
                    cmap='magma_r', axes_decor='off',
                    per_row=10, suptitle='', scalebar=False,
                    scalebar_color='white', colorbar=False,
                    padding={'top': 0.95, 'bottom': 0.05,
                             'left': 0.05, 'right':0.78})
Explanation: Plot the final results from the NMF-based segmentation
End of explanation
<END_TASK>