Cropping a domain

MuSpAn allows you to take a large domain and crop out a subset of it to form a new domain. This may be very useful when, for example, you want to restrict your focus to a specific area of your data, or if you want to see how some spatial metric varies across a large sample.

In this tutorial, we’ll explore how we can trim a large domain down into a smaller one (or many smaller ones).

Let’s start by loading in a domain, and making sure that it has some points and some shapes in it. We’ll make the shapes using alpha shapes, so that our ‘crypts’ of A and B cells are represented with a boundary shape.

import muspan as ms
import matplotlib.pyplot as plt

# Load in a domain and make some ring shapes in it
domain = ms.datasets.load_example_domain('Synthetic-Points-Architecture')
type_A_cells = ms.query.query(domain,('label','Celltype'),'is','A')
type_B_cells = ms.query.query(domain,('label','Celltype'),'is','B')
rings = ms.query.query_container(type_A_cells,'OR',type_B_cells)
new_IDs = domain.convert_objects(rings,collection_name='New object',return_IDs=True,object_type='shape',conversion_method='alpha shape',conversion_method_kwargs=dict(alpha=20))

ms.visualise.visualise(domain,'Celltype',marker_size=5)

MuSpAn domain loaded successfully. Domain summary:
Domain name: Architecture
Number of objects: 5991
Collections: ['Cell centres']
Labels: ['Celltype']
Networks: []
Distance matrices: []

(<Figure size 1000x800 with 2 Axes>, <Axes: >)

Let’s use one of those new IDs to make a new domain containing all the points inside a specific crypt. We’ll pick out one ID (object 5991, the first of the alpha shapes that we’ve just added) and see how it looks in the whole domain.

ms.visualise.visualise(domain,objects_to_plot=[new_IDs[0]],show_boundary=True)
plt.title(f'Object ID {new_IDs[0]}')

Text(0.5, 1.0, 'Object ID 5991')

Great, we know which object ID 5991 is now. Now let’s crop the domain down to that object.

outputs = ms.helpers.crop_domain(domain, [5991])
tiny_domain = outputs[5991]
ms.visualise.visualise(tiny_domain,'Celltype',show_boundary=True)

(<Figure size 1000x800 with 2 Axes>, <Axes: >)

There we have it - a brand new domain whose boundary matches the outside of the shape we just passed it. Notice a couple of things here. First, the labels in tiny_domain have been inherited from the big domain - we were able to colour the points by 'Celltype' without any difficulty.

The other, potentially more misleading thing, is that crop_domain didn’t immediately return a new domain - instead, it returned a dictionary that we called outputs, and we pulled our new domain out using the key 5991. This might seem a strange syntax when you’re cropping out a single domain, but makes much more sense if we crop out many different domains at once:

outputs = ms.helpers.crop_domain(domain, new_IDs)
print(outputs)
#ms.visualise.visualise(tiny_domain,'Celltype',show_boundary=True)

{5991: <muspan.domain.domain object at 0x311cb3ad0>, 5992: <muspan.domain.domain object at 0x3103add50>, 5993: <muspan.domain.domain object at 0x31035ae10>, 5994: <muspan.domain.domain object at 0x311f02590>, 5995: <muspan.domain.domain object at 0x17f1a1290>, 5996: <muspan.domain.domain object at 0x311e01190>, 5997: <muspan.domain.domain object at 0x311c28610>, 5998: <muspan.domain.domain object at 0x3118f38d0>, 5999: <muspan.domain.domain object at 0x311c28290>, 6000: <muspan.domain.domain object at 0x311df5350>, 6001: <muspan.domain.domain object at 0x311db10d0>, 6002: <muspan.domain.domain object at 0x3119768d0>, 6003: <muspan.domain.domain object at 0x312500250>, 6004: <muspan.domain.domain object at 0x312500550>, 6005: <muspan.domain.domain object at 0x310395a50>, 6006: <muspan.domain.domain object at 0x3117168d0>, 6007: <muspan.domain.domain object at 0x31291bed0>, 6008: <muspan.domain.domain object at 0x310358b90>, 6009: <muspan.domain.domain object at 0x3118d4210>, 6010: <muspan.domain.domain object at 0x3126f7fd0>, 6011: <muspan.domain.domain object at 0x17fe4d050>}

outputs now contains a dictionary of 21 domains, corresponding to the 21 IDs in the list of new_IDs. In general, we don’t have to pass specific object IDs here - any query like object that returns some shapes will result in a new domain for each shape.

fig, axes = plt.subplots(nrows=3,ncols=7,figsize=(21,7))
i, j = (0,0)
for key in outputs:
    this_domain = outputs[key]
    ms.visualise.visualise(this_domain,'Celltype',ax=axes[j,i],show_boundary=True)
    axes[j,i].set_title(key)
    i = i + 1
    if i > 6:
        i = 0
        j = j + 1

A situation where we might be interested in cropping the domain is when we’re using a lattice, for example. Let’s add a hexagonal lattice to our original domain, and use it to crop out a single hexagon.

hex_IDs = ms.region_based.generate_hexgrid(domain,150)
ms.visualise.visualise(domain,'ROI')

(<Figure size 1000x800 with 2 Axes>, <Axes: >)

generate_hexgrid has created hexagons with the label ‘ROI’. Let’s use a query to get the hexagon with ROI 14, and crop out a new domain with that as the boundary. Before we do, note that every object in ROI 14 has also been assigned this as a label - if we just get all the objects with ROI equal to ‘ROI_14’, we’ll find a lot of objects.

ms.visualise.visualise(domain,'ROI',objects_to_plot=('ROI','ROI 14'))

(<Figure size 1000x800 with 2 Axes>, <Axes: >)

As a result, we’ll use a more complex query to get the ID of the hexagon we’re interested in.

q14 = ms.query.query(domain,'ROI','is','ROI 14')
qHex = ms.query.query(domain,('collection',),'is','Hexgrid')
our_hexagon = q14 & qHex
ms.visualise.visualise(domain,'ROI',objects_to_plot=our_hexagon,show_boundary=True)

(<Figure size 1000x800 with 2 Axes>, <Axes: >)

When cropping to this hexagon, we have a choice to make - the hexagon intersects with some of our crypt shapes. Do we want these shapes to be excluded from the new domain? Do we want to crop them? Do we want to expand the hexagon so that they’re kept intact? Fortunately, MuSpAn allows us to choose which option we like best.

options = ['exclude','clip','expand']

fig, axes = plt.subplots(ncols=3,figsize=(21,7))
for i, option in enumerate(options):
    new_domain_dict = ms.helpers.crop_domain(domain, our_hexagon, handle_intersecting_objects=option)
    the_key = list(new_domain_dict.keys())[0]
    ms.visualise.visualise(new_domain_dict[the_key],'Celltype',ax=axes[i],show_boundary=True)
    axes[i].set_title(f'Method = {option}')
    print(f'Method {option} - Area: {new_domain_dict[the_key].boundary.area}')
    print(f'Method {option} - Number of objects: {new_domain_dict[the_key].n_objects}')

Method exclude - Area: 58456.71475544251
Method exclude - Number of objects: 252
Method clip - Area: 58456.71475544251
Method clip - Number of objects: 254
Method expand - Area: 69430.46241487304
Method expand - Number of objects: 254

Notice the key differences between the three options. exclude is the strictest, and will completely discard any shape which intersects with the object boundary. By contrast, clip leads to a domain with the same area as exclude, but any objects intersecting with the boundary will be clipped (note that this may turn one object into many, if the clipped shape is no longer connected). Finally, expand is the most tolerant method, ensuring that all our shapes remain unaltered at the expense of increasing the domain area. The best method to choose will depend on your use case.

Finally, there is one more way of cropping to this hexagon. Notice that in the very red plot a few cells up (where we visualised all objects with the label ROI equal to ROI 14), only one of the two crypt shapes shown here was shown. This is because when we generate a hexgrid, objects are assigned a single label based on the region that their centroid is contained in. We can crop the domain to these objects immediately by specifying crop_method='label' in crop_domain.

new_domain_dict = ms.helpers.crop_domain(domain,crop_method='label',label_name='ROI')
print(new_domain_dict)

{'ROI 0': <muspan.domain.domain object at 0x3169fdc50>, 'ROI 1': <muspan.domain.domain object at 0x316a39ad0>, 'ROI 10': <muspan.domain.domain object at 0x311ce6450>, 'ROI 11': <muspan.domain.domain object at 0x17f806f90>, 'ROI 12': <muspan.domain.domain object at 0x3117d6810>, 'ROI 13': <muspan.domain.domain object at 0x311d88110>, 'ROI 14': <muspan.domain.domain object at 0x17f807490>, 'ROI 15': <muspan.domain.domain object at 0x31689d950>, 'ROI 16': <muspan.domain.domain object at 0x3126f7e10>, 'ROI 17': <muspan.domain.domain object at 0x316b3f710>, 'ROI 18': <muspan.domain.domain object at 0x3118af810>, 'ROI 19': <muspan.domain.domain object at 0x316d0d7d0>, 'ROI 2': <muspan.domain.domain object at 0x316d83bd0>, 'ROI 20': <muspan.domain.domain object at 0x31162e650>, 'ROI 22': <muspan.domain.domain object at 0x3163e1f10>, 'ROI 23': <muspan.domain.domain object at 0x31162f3d0>, 'ROI 24': <muspan.domain.domain object at 0x3168d8650>, 'ROI 25': <muspan.domain.domain object at 0x316d86f50>, 'ROI 26': <muspan.domain.domain object at 0x311c3e7d0>, 'ROI 3': <muspan.domain.domain object at 0x316d86610>, 'ROI 4': <muspan.domain.domain object at 0x3164f0ed0>, 'ROI 5': <muspan.domain.domain object at 0x312923a50>, 'ROI 6': <muspan.domain.domain object at 0x311bf2d10>, 'ROI 7': <muspan.domain.domain object at 0x312dd8c10>, 'ROI 8': <muspan.domain.domain object at 0x316408d10>, 'ROI 9': <muspan.domain.domain object at 0x312e21b90>}

This has used the ROI label to crop a new domain based on every possible value of the label ROI. This is in general substantially faster than repeatedly iterating over shapes and cropping to them. Cropping by label will find all shapes with the same value of that label, and create a new domain which encloses all of them. The new domain boundary is the convex hull around the points/shapes contained within that label.

ms.visualise.visualise(new_domain_dict['ROI 14'],'Celltype',show_boundary=True)

(<Figure size 1000x800 with 2 Axes>, <Axes: >)

In this example, because we generated a hexgrid which added new shapes to the domain, we’ve picked up a giant hexagon within our new domain! If we wanted, we could now delete this (remember that the hexagon is in a different collection to the other objects, which would make it straightforward to remove). This cropping functionality is extremely useful in situations where we have local neighbourhoods - for instance, if we had assigned each cell centre a neighbourhood ID (see, for instance, the ‘Paper examples’ tutorial for making Figure 5 of the MuSpAn paper). It means that we can crop objects out of our domain even without having a shape that defines the boundary we wish to crop to. As an example, let’s make new domains which contains only points with the same cell type.

new_domain_dict = ms.helpers.crop_domain(domain,crop_method='label',label_name='Celltype')
print(new_domain_dict)

fig, axes = plt.subplots(ncols=4,figsize=(21,7))
for i, option in enumerate(['A','B','C','D']):
    ms.visualise.visualise(new_domain_dict[option],'Celltype',ax=axes[i],show_boundary=True)

{'A': <muspan.domain.domain object at 0x3123d6250>, 'B': <muspan.domain.domain object at 0x312d39a50>, 'C': <muspan.domain.domain object at 0x316ee0090>, 'D': <muspan.domain.domain object at 0x3122409d0>}

Just for fun, let’s see one of these plots again!

options = ['expand']
for i, option in enumerate(options):
    new_domain_dict = ms.helpers.crop_domain(domain, our_hexagon, handle_intersecting_objects=option)
    the_key = list(new_domain_dict.keys())[0]
    ms.visualise.visualise(new_domain_dict[the_key],'Celltype',show_boundary=True)
    plt.gca().set_title(f'Method = {option}')
    print(f'Method {option} - Area: {new_domain_dict[the_key].boundary.area}')
    print(f'Method {option} - Number of objects: {new_domain_dict[the_key].n_objects}')

Method expand - Area: 69430.46241487304
Method expand - Number of objects: 254