Add some code for finding He bubbles

%% Cell type:code id:0e02bb89-8003-47a2-839f-58536b067490 tags:

``` python

# Some parameters for the behavior of the extraction

# Specifies the spacing of the sampling grid.

lattice_parameter = 3.177

# The radius to use for a He splat

He_radius = 0.01

# Number of He atoms that have to escape for a bubble to be

# considered burst

burst_size = 100

# Directory containing VTK files holding atoms

#data_path = '/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm/vtk-files'

data_path = '/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files'

```

%% Cell type:code id:21fe4234-a349-402d-afdf-c9cb9a04f36a tags:

``` python

import sys

sys.path.append('/Users/4d5/builds/VTK/Debug/lib/python3.11/site-packages')

import vtk

```

%% Cell type:code id:9eaffe57-7c2c-401b-a253-7ee2628e1160 tags:

``` python

def read_he_atoms(filename):

    '''Given a vtk file containing points from a LAMMPS atom file, creates a

point cloud containing a point and vertex cell for each He atom. All W

atoms are removed. This is returned as a `vtkUnstructuredGrid`.'''

    reader = vtk.vtkDataSetReader()

    reader.SetFileName(filename)

    tocloud = vtk.vtkConvertToPointCloud()

    tocloud.SetInputConnection(reader.GetOutputPort())

    tocloud.SetCellGenerationMode(tocloud.VERTEX_CELLS)

    extract_he = vtk.vtkThreshold()

    extract_he.SetInputConnection(tocloud.GetOutputPort())

    extract_he.SetInputArrayToProcess(

        0, 0, 0, 'vtkDataObject::FIELD_ASSOCIATION_POINTS', 'type')

    extract_he.SetLowerThreshold(1.9)

    extract_he.SetUpperThreshold(2.1)

    extract_he.Update()

    return extract_he.GetOutput()

```

%% Cell type:code id:07203a73-966e-4e20-a766-b5cf27f43f59 tags:

``` python

def find_bubbles(atoms):

    '''Given a VTK data object containing a point cloud of atoms, identify

which atoms together form a bubble by clustering nearby atoms. Each

bubble is given a unique identifier. This function returns a VTK data

object of the same structure as the input with a `RegionId` field added

that identifies which bubble each atom belongs to.'''

    # If there are no atoms, then there is nothing to find.

    if atoms.GetNumberOfPoints() < 1:

        return atoms

    # Create a density field of atoms by splatting the atoms into a grid.

    # Here we use `vtkGaussianSplatter`. If this is taking too long, it

    # is possible to speed things up using a `vtkFastSplatter`.

    bounds = atoms.GetBounds()

    sampling_resolution = (

        int((bounds[1] - bounds[0]) / (lattice_parameter / 2)),

        int((bounds[3] - bounds[2]) / (lattice_parameter / 2)),

        int((bounds[5] - bounds[4]) / (lattice_parameter / 2)))

    splatter = vtk.vtkGaussianSplatter()

    splatter.SetInputDataObject(atoms)

    splatter.SetSampleDimensions(sampling_resolution)

    splatter.SetRadius(He_radius)

    splatter.SetNormalWarping(0)

    splatter.SetAccumulationModeToSum()

    # Use an image connectivity filter to find regions with nearby atoms.

    # The filter will assign each connected region a unique ID and write

    # that to point data.

    connectivity = vtk.vtkImageConnectivityFilter()

    connectivity.SetInputConnection(splatter.GetOutputPort())

    connectivity.SetLabelScalarTypeToInt()

    #Lowering threshold to include regions to prevent missing atoms at edge

    # of bubble.

    connectivity.SetScalarRange(0.1, 100000000000)

    # Frustratingly, there appears to be a bug in `vtkImageConnectivityFilter`

    # that does not properly pass the origin and spacing of the image, which

    # means that it will no longer align with the atom data. Restore it to

    # the proper values.

    splatter.Update()

    splatter_output = splatter.GetOutput()

    reset_position = vtk.vtkImageChangeInformation()

    reset_position.SetInputConnection(connectivity.GetOutputPort())

    reset_position.SetOutputOrigin(splatter_output.GetOrigin())

    reset_position.SetOutputSpacing(splatter_output.GetSpacing())

    # Sample the connectivity field back onto the atoms that generated them.

    # This will add a `RegionId` field to the atoms that identifes which

    # cluster each one belongs to.

    probe = vtk.vtkProbeFilter()

    probe.SetInputDataObject(0, atoms)

    probe.SetInputConnection(1, reset_position.GetOutputPort())

    probe.CategoricalDataOn()

    probe.PassPointArraysOn()

    probe.Update()

    return probe.GetOutput()

```

%% Cell type:code id:0c64f950-83bb-4c26-9f13-87aa550ffc3f tags:

``` python

def make_bubble_dict(atoms):

    '''Given a VTK data object containing a point cloud of atoms and a

`RegionId` field identifying a containing bubble for each atom (e.g.

atoms read with `read_he_atoms` and then processed with `find_bubbles`),

returns a dictionary with a key for each bubble id. Each value is a set

of the ids of atoms in that bubble.'''

    if atoms.GetNumberOfPoints() < 1:

        return {}

    bubble_ids = atoms.GetPointData().GetArray('RegionId')

    atom_ids = atoms.GetPointData().GetArray('id')

    bubbles = {}

    for bubble_id in range(1, bubble_ids.GetValueRange()[1] + 1):

        bubbles[bubble_id] = set()

    for index in range(bubble_ids.GetNumberOfValues()):

        bubble_id = bubble_ids.GetValue(index)

        # Note: if bubble_id == 0, then the atom is not part of a bubble

        if bubble_id > 0:

            atom_id = int(atom_ids.GetTuple1(index))

            bubbles[bubble_id].add(atom_id)

    return bubbles

```

%% Cell type:code id:e11aa19a-4946-4468-80dd-13546d7a2b68 tags:

``` python

def find_burst_bubbles(bubbles, next_atoms):

    '''Given a dictionary of bubbles (created by `make_bubble_dict` and a VTK

data object containing the atoms for the next timestep, find bubbles that

appear to have burst. When a bubble bursts, the atoms become evacuated

from their cavity into a vacuum and disappear from the simulation. When

a large number of atoms in a bubble are not in the next timestep, that is

indicitive of the bubble bursting.

The function returns a dictionary with the keys being ids of bubbles that

have burst and the values being the number of atoms that went missing.'''

    # Create a set of all the atom ids in the next timestep for easy finding.

    atom_ids = next_atoms.GetPointData().GetArray('id')

    next_atom_set = set()

    for index in range(atom_ids.GetNumberOfValues()):

        next_atom_set.add(int(atom_ids.GetTuple1(index)))

    # Iterate over all bubbles from the previous timestep and count how many

    # have suddenly disappeared. If the count is above the `burst_size`

    # threshold, record it.

    burst_bubbles = {}

    for bubble_id, bubble_set in bubbles.items():

        if len(bubble_set) < burst_size:

            continue

        num_gone = 0

        for atom_id in bubble_set:

            if not atom_id in next_atom_set:

                num_gone += 1

        if num_gone >= burst_size:

            burst_bubbles[bubble_id] = num_gone

    return burst_bubbles

```

%% Cell type:code id:0f551511-82e4-4e98-8557-6315b7beacb3 tags:

``` python

# Iterate over all the files in the data path and find potential bubble bursts.

bubble_dict = None

max_loss = 0

max_loss_ratio = 0.0

import glob

for filename in sorted(glob.glob(data_path + '/*.vtk')):

    he_atoms = read_he_atoms(filename)

    if bubble_dict:

        burst_bubbles = find_burst_bubbles(bubble_dict, he_atoms)

        if len(burst_bubbles) > 0:

            print(filename)

            print(burst_bubbles)

        for bubble_id, loss in burst_bubbles.items():

            max_loss = max(max_loss, loss)

            loss_ratio = float(loss)/float(len(bubble_dict[bubble_id]))

            max_loss_ratio = max(max_loss_ratio, loss_ratio)

    bubbles = find_bubbles(he_atoms)

    bubble_dict = make_bubble_dict(bubbles)

print(max_loss)

print(max_loss_ratio)

```

%% Output

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-084475.vtk

    {100: 167}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-084476.vtk

    {21: 151}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-091829.vtk

    {64: 164}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-091830.vtk

    {26: 359}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-091831.vtk

    {20: 185}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-091832.vtk

    {13: 113}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-094432.vtk

    {52: 299}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-094433.vtk

    {22: 205}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-094434.vtk

    {61: 110}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-095042.vtk

    {58: 102}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-098328.vtk

    {126: 351}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-098329.vtk

    {33: 380}

    /Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm_detailed/vtk-files/110-25x25nm-detailed-098330.vtk

    {32: 160}

    380

    0.8429752066115702

%% Cell type:code id:8c9e9b81-4315-49d6-b409-a9e45c8ffbc6 tags:

``` python

```

%% Cell type:markdown id:a16f8887-aad2-4c98-84b3-77ee765c93b0 tags:

Some parameters for the behavior of the extraction.

%% Cell type:code id:663024d6-3321-41f6-86d8-cf56ea22b290 tags:

``` python

# Specifies the spacing of the sampling grid.

lattice_parameter = 3.177

# The radius to use for a He splat

He_radius = 0.01

# Number of He atoms that have to escape for a bubble to be

# considered burst

burst_size = 1

```

%% Cell type:code id:994f58bf-ce63-4c35-918f-69dd2fb784d6 tags:

``` python

import sys

sys.path.append('/Users/4d5/builds/VTK/Debug/lib/python3.11/site-packages')

import vtk

```

%% Cell type:markdown id:a36e89f3-fe92-44e1-96ac-b828494fea58 tags:

Read in a data file of atoms. This is an unstructured grid, but it has no cells. Just points for each atom.

%% Cell type:code id:8af6940d-a1a5-451a-9bd1-2fb967fa28e7 tags:

``` python

reader = vtk.vtkDataSetReader()

reader.SetFileName('/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm/vtk-files/110-25x25nm100000.vtk')

```

%% Cell type:markdown id:86d36655-7c2a-4459-b9e5-e58734a93638 tags:

Convert the data to a point cloud, which will add a vertex cell for each point. This will allow us to apply filters like Threshold.

%% Cell type:code id:0266fccc-3ad7-4bd5-9c30-1beb1badb20f tags:

``` python

tocloud = vtk.vtkConvertToPointCloud()

tocloud.SetInputConnection(reader.GetOutputPort())

tocloud.SetCellGenerationMode(tocloud.VERTEX_CELLS)

```

%% Cell type:markdown id:0025f269-3c5a-4b68-9112-09b1393b8f0b tags:

Use a threshold to get just the He atoms (type = 2).

%% Cell type:code id:e6470202-0b28-41a5-838a-f3c2ffadb5f2 tags:

``` python

extract_he = vtk.vtkThreshold()

extract_he.SetInputConnection(tocloud.GetOutputPort())

extract_he.SetInputArrayToProcess(0, 0, 0, 'vtkDataObject::FIELD_ASSOCIATION_POINTS', 'type')

extract_he.SetLowerThreshold(1.9)

extract_he.SetUpperThreshold(2.1)

```

%% Cell type:markdown id:304b629d-8a99-4a80-8068-93753720864b tags:

Create a density field of He atoms by splatting the atoms into a grid. Here we use `vtkGaussianSplatter`. Some previous code by Ollie has a faster version that uses `vtkFastSplatter`. If we want to speed things up, that could be an easy way to go.

%% Cell type:code id:a44f124c-f9be-422c-b813-5c529d5a3c35 tags:

``` python

tocloud.Update()

bounds = tocloud.GetOutput().GetBounds()

sampling_resolution = (int((bounds[1] - bounds[0]) / (lattice_parameter / 2)),

                       int((bounds[3] - bounds[2]) / (lattice_parameter / 2)),

                       int((bounds[5] - bounds[4]) / (lattice_parameter / 2)))

```

%% Cell type:code id:3768f452-c720-4ef1-b8ae-966c1b449a3a tags:

``` python

splatter = vtk.vtkGaussianSplatter()

splatter.SetInputConnection(extract_he.GetOutputPort())

splatter.SetSampleDimensions(sampling_resolution)

splatter.SetRadius(He_radius)

splatter.SetNormalWarping(0)

splatter.SetAccumulationModeToSum()

```

%% Cell type:markdown id:91106a3a-fbbe-4c4c-ac9d-4b0fda722553 tags:

Use an image connectivity filter to find regions with nearby He atoms. The filter will assign each connected region a unique ID and write that to point data.

%% Cell type:code id:0abf1673-cf56-4287-ac85-7cc4346e4bb8 tags:

``` python

connectivity = vtk.vtkImageConnectivityFilter()

connectivity.SetInputConnection(splatter.GetOutputPort())

connectivity.SetLabelScalarTypeToInt()

#Lowering threshold to include regions to prevent missing atoms at edge of bubble.

connectivity.SetScalarRange(0.1, 100000000000)

```

%% Cell type:markdown id:bb71178d-5ba7-4120-a89c-cebd8df0a5bd tags:

Frustratingly, there appears to be a bug in `vtkImageConnectivityFilter` that does not properly pass the origin and spacing of the image, which means it will no longer align with the atom data. Restore it to the proper values.

%% Cell type:code id:7617df5d-0d51-446d-aed5-ee5693d4990b tags:

``` python

splatter.Update()

splatter_output = splatter.GetOutput()

reset_position = vtk.vtkImageChangeInformation()

reset_position.SetInputConnection(connectivity.GetOutputPort())

reset_position.SetOutputOrigin(splatter_output.GetOrigin())

reset_position.SetOutputSpacing(splatter_output.GetSpacing())

```

%% Cell type:markdown id:7227ea72-f31b-437a-894f-8bc328e25d71 tags:

Sample the connectivity field back onto the atoms that generated them. This will add a `RegionId` field to the atoms that identifies which cluster each one belongs to.

%% Cell type:code id:a5fa732f-ddc6-4eda-be34-adb036417ab1 tags:

``` python

probe = vtk.vtkProbeFilter()

probe.SetInputConnection(0, extract_he.GetOutputPort())

probe.SetInputConnection(1, reset_position.GetOutputPort())

probe.CategoricalDataOn()

probe.PassPointArraysOn()

```

%% Cell type:code id:d05c411a-ae1f-4b4f-90cd-0693809f4d3d tags:

``` python

writer = vtk.vtkDataSetWriter()

writer.SetInputConnection(probe.GetOutputPort())

writer.SetFileTypeToBinary()

writer.SetFileName('/Users/4d5/Downloads/test.vtk')

writer.Write()

```

%% Output

    1

%% Cell type:code id:eb3af2c0-1e42-4085-834d-b5c5fe427028 tags:

``` python

# Create a dictionary with a key for each bubble. Each value is a set

# of the atoms in that bubble (by their ids).

probe.Update()

clustered_he = probe.GetOutput()

bubble_ids = clustered_he.GetPointData().GetArray('RegionId')

he_ids = clustered_he.GetPointData().GetArray('id')

bubbles = {}

for bubble_id in range(1, bubble_ids.GetValueRange()[1] + 1):

    bubbles[bubble_id] = set()

for atom_id in range(bubble_ids.GetNumberOfValues()):

    bubble_id = bubble_ids.GetValue(atom_id)

    if bubble_id > 0:

        he_id = int(he_ids.GetTuple1(atom_id))

        bubbles[bubble_id].add(he_id)

```

%% Cell type:code id:98ac76bf-27b6-4f8c-952d-ac4246026262 tags:

``` python

# Read in the atoms for the next timestep

reader = vtk.vtkDataSetReader()

reader.SetFileName('/Users/4d5/data/PSI/karl-hammond-2023/110-25x25nm/vtk-files/110-25x25nm100250.vtk')

tocloud = vtk.vtkConvertToPointCloud()

tocloud.SetInputConnection(reader.GetOutputPort())

tocloud.SetCellGenerationMode(tocloud.VERTEX_CELLS)

extract_he = vtk.vtkThreshold()

extract_he.SetInputConnection(tocloud.GetOutputPort())

extract_he.SetInputArrayToProcess(0, 0, 0, 'vtkDataObject::FIELD_ASSOCIATION_POINTS', 'type')

extract_he.SetLowerThreshold(1.9)

extract_he.SetUpperThreshold(2.1)

extract_he.Update()

next_data = extract_he.GetOutput()

```

%% Cell type:code id:2dd6660b-3929-42aa-bc75-36afff54a49d tags:

``` python

# Create a set of all the He ids in the next timestep

he_ids = next_data.GetPointData().GetArray('id')

he_set = set()

for atom_id in range(he_ids.GetNumberOfValues()):

    he_set.add(int(he_ids.GetTuple1(atom_id)))

```

%% Cell type:code id:cfae5559-5259-40a5-8c01-7c3747ca4c51 tags:

``` python

# Iterate over all bubbles from the previous timestep and see if there

# are any where a lot of atoms have suddenly disappeared.

for bubble_id, bubble_set in bubbles.items():

    if len(bubble_set) < burst_size:

        continue

    num_gone = 0

    for he_id in bubble_set:

        if not he_id in he_set:

            num_gone += 1

    if (num_gone > 0):

        print(bubble_id, num_gone)

```

%% Output

    33 1

    58 1

    62 1

    83 1

    102 3

    132 4

    149 32

    234 1

    245 10

    309 2

    357 3

    384 2

    396 1

    430 1

    447 1

    460 2

    477 10

    518 2

%% Cell type:code id:0a835232-7356-4c0e-b05b-b776ffe66683 tags:

``` python

```