Aligning partial coronal brain sections with the Allen Brain Altas
In this notebook, we align a STARmap dataset of a the right hemisphere of a mouse brain to the Allen Brain Atlas.
We will use STalign to achieve this alignment. We will first load the relevant code libraries.
[1]:
#Import dependencies
import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt
import pandas as pd # for csv.
from matplotlib import cm
from matplotlib.lines import Line2D
import os
from os.path import exists,split,join,splitext
from os import makedirs
import glob
import requests
from collections import defaultdict
import nrrd
import torch
from torch.nn.functional import grid_sample
import tornado
from STalign import STalign
import copy
import pandas as pd
We have already downloaded STARmap datasets and stored them in a folder called starmap_data.
Now, we can read in the cell infomation using pandas as pd.
We must load all cell x and y positions into the variable x and y respectively. Importantly, the scale of the loaded data must approximately match the Allen Brain Atlas, where ~1 pixel unit is ~1 um
[2]:
#Load file
filename = '__'
df = pd.read_csv(r"/home/manju/Documents/STalign-master/STalign_build/docs/starmap_data/well11_spatial.csv.gz")
#Load x position
x = np.array(df['X'])[1:].astype(float) #change to x positions of cells
#Load y position
y = np.array(df['Y'])[1:].astype(float) #change to column y positions of cells
[3]:
dx=50
blur = 1
Now we can load in data from the Allen Brain Altas.
We need to load in and save files containing the ontology of regions, a 3D cell density atlas, and annotations of brain regions in the aforementioned 3D atlas.
Depending on which Allen Atlas you would like to use, you can change the url that you download these files from. In this example, we are aligning the adult mouse (P56) Allen Atlas.
[4]:
url = 'http://api.brain-map.org/api/v2/data/query.csv?criteria=model::Structure,rma::criteria,[ontology_id$eq1],rma::options[order$eq%27structures.graph_order%27][num_rows$eqall]'
ontology_name,namesdict = STalign.download_aba_ontology(url, 'allen_ontology.csv') #url for adult mouse
<Response [200]>
[5]:
imageurl = 'http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/ara_nissl/ara_nissl_50.nrrd'
labelurl = 'http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation/ccf_2017/annotation_50.nrrd'
imagefile, labelfile = STalign.download_aba_image_labels(imageurl, labelurl, 'aba_nissl.nrrd', 'aba_annotation.nrrd')
Our first step in this alignment will be to rasterize the single cell centroid positions into an image. This converts our x and y arrays of cells into a pixelated image. For this example, we choose to rasterize at a resolution of 15 (um).
[6]:
#Rasterize Image
X_,Y_,W = STalign.rasterize(x,y,dx=dx, blur = blur,draw=False)
[7]:
X_
[7]:
array([-1.42879227e+03, -1.37879227e+03, -1.32879227e+03, -1.27879227e+03,
-1.22879227e+03, -1.17879227e+03, -1.12879227e+03, -1.07879227e+03,
-1.02879227e+03, -9.78792270e+02, -9.28792270e+02, -8.78792270e+02,
-8.28792270e+02, -7.78792270e+02, -7.28792270e+02, -6.78792270e+02,
-6.28792270e+02, -5.78792270e+02, -5.28792270e+02, -4.78792270e+02,
-4.28792270e+02, -3.78792270e+02, -3.28792270e+02, -2.78792270e+02,
-2.28792270e+02, -1.78792270e+02, -1.28792270e+02, -7.87922701e+01,
-2.87922701e+01, 2.12077299e+01, 7.12077299e+01, 1.21207730e+02,
1.71207730e+02, 2.21207730e+02, 2.71207730e+02, 3.21207730e+02,
3.71207730e+02, 4.21207730e+02, 4.71207730e+02, 5.21207730e+02,
5.71207730e+02, 6.21207730e+02, 6.71207730e+02, 7.21207730e+02,
7.71207730e+02, 8.21207730e+02, 8.71207730e+02, 9.21207730e+02,
9.71207730e+02, 1.02120773e+03, 1.07120773e+03, 1.12120773e+03,
1.17120773e+03, 1.22120773e+03, 1.27120773e+03, 1.32120773e+03,
1.37120773e+03, 1.42120773e+03, 1.47120773e+03, 1.52120773e+03,
1.57120773e+03, 1.62120773e+03, 1.67120773e+03, 1.72120773e+03,
1.77120773e+03, 1.82120773e+03, 1.87120773e+03, 1.92120773e+03,
1.97120773e+03, 2.02120773e+03, 2.07120773e+03, 2.12120773e+03,
2.17120773e+03, 2.22120773e+03, 2.27120773e+03, 2.32120773e+03,
2.37120773e+03, 2.42120773e+03, 2.47120773e+03, 2.52120773e+03,
2.57120773e+03, 2.62120773e+03, 2.67120773e+03, 2.72120773e+03,
2.77120773e+03, 2.82120773e+03, 2.87120773e+03, 2.92120773e+03,
2.97120773e+03, 3.02120773e+03, 3.07120773e+03, 3.12120773e+03,
3.17120773e+03, 3.22120773e+03, 3.27120773e+03, 3.32120773e+03,
3.37120773e+03, 3.42120773e+03, 3.47120773e+03, 3.52120773e+03,
3.57120773e+03, 3.62120773e+03, 3.67120773e+03, 3.72120773e+03,
3.77120773e+03, 3.82120773e+03, 3.87120773e+03, 3.92120773e+03,
3.97120773e+03, 4.02120773e+03, 4.07120773e+03, 4.12120773e+03,
4.17120773e+03, 4.22120773e+03, 4.27120773e+03, 4.32120773e+03,
4.37120773e+03, 4.42120773e+03, 4.47120773e+03, 4.52120773e+03,
4.57120773e+03, 4.62120773e+03, 4.67120773e+03, 4.72120773e+03,
4.77120773e+03, 4.82120773e+03, 4.87120773e+03, 4.92120773e+03,
4.97120773e+03, 5.02120773e+03, 5.07120773e+03, 5.12120773e+03,
5.17120773e+03, 5.22120773e+03, 5.27120773e+03, 5.32120773e+03,
5.37120773e+03, 5.42120773e+03, 5.47120773e+03, 5.52120773e+03,
5.57120773e+03, 5.62120773e+03, 5.67120773e+03, 5.72120773e+03,
5.77120773e+03, 5.82120773e+03, 5.87120773e+03, 5.92120773e+03,
5.97120773e+03, 6.02120773e+03, 6.07120773e+03, 6.12120773e+03,
6.17120773e+03, 6.22120773e+03, 6.27120773e+03, 6.32120773e+03,
6.37120773e+03, 6.42120773e+03, 6.47120773e+03, 6.52120773e+03,
6.57120773e+03, 6.62120773e+03, 6.67120773e+03, 6.72120773e+03,
6.77120773e+03, 6.82120773e+03, 6.87120773e+03, 6.92120773e+03,
6.97120773e+03, 7.02120773e+03, 7.07120773e+03, 7.12120773e+03,
7.17120773e+03, 7.22120773e+03, 7.27120773e+03, 7.32120773e+03,
7.37120773e+03, 7.42120773e+03, 7.47120773e+03, 7.52120773e+03,
7.57120773e+03, 7.62120773e+03, 7.67120773e+03, 7.72120773e+03,
7.77120773e+03, 7.82120773e+03, 7.87120773e+03, 7.92120773e+03,
7.97120773e+03, 8.02120773e+03, 8.07120773e+03, 8.12120773e+03,
8.17120773e+03, 8.22120773e+03, 8.27120773e+03, 8.32120773e+03,
8.37120773e+03, 8.42120773e+03, 8.47120773e+03, 8.52120773e+03,
8.57120773e+03, 8.62120773e+03, 8.67120773e+03, 8.72120773e+03,
8.77120773e+03, 8.82120773e+03, 8.87120773e+03, 8.92120773e+03,
8.97120773e+03, 9.02120773e+03, 9.07120773e+03, 9.12120773e+03,
9.17120773e+03, 9.22120773e+03, 9.27120773e+03, 9.32120773e+03,
9.37120773e+03, 9.42120773e+03, 9.47120773e+03, 9.52120773e+03,
9.57120773e+03, 9.62120773e+03, 9.67120773e+03, 9.72120773e+03,
9.77120773e+03, 9.82120773e+03, 9.87120773e+03, 9.92120773e+03,
9.97120773e+03, 1.00212077e+04, 1.00712077e+04, 1.01212077e+04,
1.01712077e+04, 1.02212077e+04, 1.02712077e+04, 1.03212077e+04,
1.03712077e+04, 1.04212077e+04, 1.04712077e+04, 1.05212077e+04,
1.05712077e+04, 1.06212077e+04, 1.06712077e+04, 1.07212077e+04,
1.07712077e+04, 1.08212077e+04, 1.08712077e+04, 1.09212077e+04,
1.09712077e+04, 1.10212077e+04, 1.10712077e+04, 1.11212077e+04,
1.11712077e+04, 1.12212077e+04, 1.12712077e+04, 1.13212077e+04,
1.13712077e+04, 1.14212077e+04, 1.14712077e+04, 1.15212077e+04,
1.15712077e+04, 1.16212077e+04, 1.16712077e+04, 1.17212077e+04,
1.17712077e+04, 1.18212077e+04, 1.18712077e+04, 1.19212077e+04,
1.19712077e+04, 1.20212077e+04, 1.20712077e+04, 1.21212077e+04,
1.21712077e+04, 1.22212077e+04, 1.22712077e+04, 1.23212077e+04,
1.23712077e+04, 1.24212077e+04, 1.24712077e+04, 1.25212077e+04,
1.25712077e+04, 1.26212077e+04, 1.26712077e+04, 1.27212077e+04,
1.27712077e+04, 1.28212077e+04, 1.28712077e+04, 1.29212077e+04,
1.29712077e+04, 1.30212077e+04, 1.30712077e+04, 1.31212077e+04,
1.31712077e+04, 1.32212077e+04, 1.32712077e+04, 1.33212077e+04,
1.33712077e+04, 1.34212077e+04, 1.34712077e+04, 1.35212077e+04,
1.35712077e+04, 1.36212077e+04, 1.36712077e+04, 1.37212077e+04,
1.37712077e+04, 1.38212077e+04, 1.38712077e+04, 1.39212077e+04,
1.39712077e+04, 1.40212077e+04, 1.40712077e+04, 1.41212077e+04,
1.41712077e+04, 1.42212077e+04, 1.42712077e+04, 1.43212077e+04,
1.43712077e+04, 1.44212077e+04, 1.44712077e+04, 1.45212077e+04,
1.45712077e+04, 1.46212077e+04, 1.46712077e+04, 1.47212077e+04,
1.47712077e+04, 1.48212077e+04, 1.48712077e+04, 1.49212077e+04,
1.49712077e+04, 1.50212077e+04, 1.50712077e+04, 1.51212077e+04,
1.51712077e+04, 1.52212077e+04, 1.52712077e+04, 1.53212077e+04,
1.53712077e+04, 1.54212077e+04, 1.54712077e+04, 1.55212077e+04,
1.55712077e+04, 1.56212077e+04, 1.56712077e+04, 1.57212077e+04,
1.57712077e+04, 1.58212077e+04, 1.58712077e+04, 1.59212077e+04,
1.59712077e+04, 1.60212077e+04, 1.60712077e+04, 1.61212077e+04,
1.61712077e+04, 1.62212077e+04, 1.62712077e+04, 1.63212077e+04,
1.63712077e+04, 1.64212077e+04, 1.64712077e+04, 1.65212077e+04,
1.65712077e+04, 1.66212077e+04, 1.66712077e+04, 1.67212077e+04,
1.67712077e+04, 1.68212077e+04, 1.68712077e+04, 1.69212077e+04,
1.69712077e+04, 1.70212077e+04, 1.70712077e+04, 1.71212077e+04,
1.71712077e+04, 1.72212077e+04, 1.72712077e+04, 1.73212077e+04,
1.73712077e+04, 1.74212077e+04, 1.74712077e+04, 1.75212077e+04,
1.75712077e+04, 1.76212077e+04, 1.76712077e+04, 1.77212077e+04,
1.77712077e+04, 1.78212077e+04, 1.78712077e+04, 1.79212077e+04,
1.79712077e+04, 1.80212077e+04, 1.80712077e+04, 1.81212077e+04,
1.81712077e+04, 1.82212077e+04, 1.82712077e+04, 1.83212077e+04,
1.83712077e+04, 1.84212077e+04, 1.84712077e+04, 1.85212077e+04,
1.85712077e+04, 1.86212077e+04, 1.86712077e+04, 1.87212077e+04,
1.87712077e+04, 1.88212077e+04, 1.88712077e+04, 1.89212077e+04,
1.89712077e+04, 1.90212077e+04, 1.90712077e+04, 1.91212077e+04,
1.91712077e+04, 1.92212077e+04, 1.92712077e+04, 1.93212077e+04,
1.93712077e+04, 1.94212077e+04, 1.94712077e+04, 1.95212077e+04,
1.95712077e+04, 1.96212077e+04, 1.96712077e+04, 1.97212077e+04,
1.97712077e+04, 1.98212077e+04, 1.98712077e+04, 1.99212077e+04,
1.99712077e+04, 2.00212077e+04, 2.00712077e+04, 2.01212077e+04,
2.01712077e+04, 2.02212077e+04, 2.02712077e+04, 2.03212077e+04,
2.03712077e+04, 2.04212077e+04, 2.04712077e+04, 2.05212077e+04,
2.05712077e+04, 2.06212077e+04, 2.06712077e+04, 2.07212077e+04,
2.07712077e+04, 2.08212077e+04, 2.08712077e+04, 2.09212077e+04,
2.09712077e+04, 2.10212077e+04, 2.10712077e+04, 2.11212077e+04,
2.11712077e+04, 2.12212077e+04, 2.12712077e+04, 2.13212077e+04,
2.13712077e+04, 2.14212077e+04, 2.14712077e+04, 2.15212077e+04,
2.15712077e+04, 2.16212077e+04, 2.16712077e+04, 2.17212077e+04,
2.17712077e+04, 2.18212077e+04, 2.18712077e+04, 2.19212077e+04,
2.19712077e+04, 2.20212077e+04, 2.20712077e+04, 2.21212077e+04,
2.21712077e+04, 2.22212077e+04, 2.22712077e+04, 2.23212077e+04,
2.23712077e+04, 2.24212077e+04, 2.24712077e+04, 2.25212077e+04,
2.25712077e+04, 2.26212077e+04, 2.26712077e+04, 2.27212077e+04,
2.27712077e+04, 2.28212077e+04, 2.28712077e+04, 2.29212077e+04,
2.29712077e+04, 2.30212077e+04, 2.30712077e+04, 2.31212077e+04,
2.31712077e+04, 2.32212077e+04, 2.32712077e+04, 2.33212077e+04,
2.33712077e+04, 2.34212077e+04, 2.34712077e+04, 2.35212077e+04,
2.35712077e+04, 2.36212077e+04, 2.36712077e+04, 2.37212077e+04,
2.37712077e+04, 2.38212077e+04, 2.38712077e+04, 2.39212077e+04,
2.39712077e+04, 2.40212077e+04, 2.40712077e+04, 2.41212077e+04,
2.41712077e+04, 2.42212077e+04, 2.42712077e+04, 2.43212077e+04,
2.43712077e+04, 2.44212077e+04, 2.44712077e+04, 2.45212077e+04,
2.45712077e+04, 2.46212077e+04, 2.46712077e+04, 2.47212077e+04,
2.47712077e+04, 2.48212077e+04, 2.48712077e+04, 2.49212077e+04,
2.49712077e+04, 2.50212077e+04, 2.50712077e+04, 2.51212077e+04,
2.51712077e+04, 2.52212077e+04, 2.52712077e+04, 2.53212077e+04,
2.53712077e+04, 2.54212077e+04, 2.54712077e+04, 2.55212077e+04,
2.55712077e+04, 2.56212077e+04, 2.56712077e+04, 2.57212077e+04,
2.57712077e+04, 2.58212077e+04, 2.58712077e+04, 2.59212077e+04,
2.59712077e+04, 2.60212077e+04, 2.60712077e+04, 2.61212077e+04,
2.61712077e+04, 2.62212077e+04, 2.62712077e+04, 2.63212077e+04,
2.63712077e+04, 2.64212077e+04, 2.64712077e+04, 2.65212077e+04,
2.65712077e+04, 2.66212077e+04, 2.66712077e+04, 2.67212077e+04,
2.67712077e+04, 2.68212077e+04, 2.68712077e+04, 2.69212077e+04,
2.69712077e+04, 2.70212077e+04, 2.70712077e+04, 2.71212077e+04,
2.71712077e+04, 2.72212077e+04, 2.72712077e+04, 2.73212077e+04,
2.73712077e+04, 2.74212077e+04, 2.74712077e+04, 2.75212077e+04,
2.75712077e+04, 2.76212077e+04, 2.76712077e+04, 2.77212077e+04,
2.77712077e+04, 2.78212077e+04, 2.78712077e+04, 2.79212077e+04,
2.79712077e+04, 2.80212077e+04, 2.80712077e+04, 2.81212077e+04,
2.81712077e+04, 2.82212077e+04, 2.82712077e+04, 2.83212077e+04,
2.83712077e+04, 2.84212077e+04, 2.84712077e+04, 2.85212077e+04,
2.85712077e+04, 2.86212077e+04, 2.86712077e+04, 2.87212077e+04,
2.87712077e+04, 2.88212077e+04, 2.88712077e+04, 2.89212077e+04,
2.89712077e+04, 2.90212077e+04, 2.90712077e+04, 2.91212077e+04,
2.91712077e+04, 2.92212077e+04, 2.92712077e+04, 2.93212077e+04,
2.93712077e+04, 2.94212077e+04, 2.94712077e+04, 2.95212077e+04,
2.95712077e+04, 2.96212077e+04, 2.96712077e+04, 2.97212077e+04,
2.97712077e+04, 2.98212077e+04, 2.98712077e+04, 2.99212077e+04,
2.99712077e+04, 3.00212077e+04, 3.00712077e+04, 3.01212077e+04,
3.01712077e+04, 3.02212077e+04, 3.02712077e+04, 3.03212077e+04,
3.03712077e+04, 3.04212077e+04, 3.04712077e+04, 3.05212077e+04,
3.05712077e+04, 3.06212077e+04, 3.06712077e+04, 3.07212077e+04,
3.07712077e+04, 3.08212077e+04, 3.08712077e+04, 3.09212077e+04,
3.09712077e+04, 3.10212077e+04, 3.10712077e+04, 3.11212077e+04,
3.11712077e+04, 3.12212077e+04, 3.12712077e+04, 3.13212077e+04,
3.13712077e+04, 3.14212077e+04, 3.14712077e+04, 3.15212077e+04,
3.15712077e+04, 3.16212077e+04, 3.16712077e+04, 3.17212077e+04,
3.17712077e+04, 3.18212077e+04, 3.18712077e+04, 3.19212077e+04,
3.19712077e+04, 3.20212077e+04, 3.20712077e+04, 3.21212077e+04,
3.21712077e+04, 3.22212077e+04, 3.22712077e+04, 3.23212077e+04,
3.23712077e+04, 3.24212077e+04, 3.24712077e+04, 3.25212077e+04,
3.25712077e+04, 3.26212077e+04, 3.26712077e+04, 3.27212077e+04,
3.27712077e+04, 3.28212077e+04, 3.28712077e+04, 3.29212077e+04,
3.29712077e+04, 3.30212077e+04, 3.30712077e+04, 3.31212077e+04,
3.31712077e+04, 3.32212077e+04, 3.32712077e+04, 3.33212077e+04,
3.33712077e+04, 3.34212077e+04, 3.34712077e+04, 3.35212077e+04,
3.35712077e+04, 3.36212077e+04, 3.36712077e+04, 3.37212077e+04,
3.37712077e+04, 3.38212077e+04, 3.38712077e+04, 3.39212077e+04,
3.39712077e+04, 3.40212077e+04, 3.40712077e+04, 3.41212077e+04,
3.41712077e+04, 3.42212077e+04, 3.42712077e+04, 3.43212077e+04,
3.43712077e+04, 3.44212077e+04, 3.44712077e+04, 3.45212077e+04,
3.45712077e+04, 3.46212077e+04, 3.46712077e+04, 3.47212077e+04,
3.47712077e+04, 3.48212077e+04, 3.48712077e+04, 3.49212077e+04,
3.49712077e+04, 3.50212077e+04, 3.50712077e+04, 3.51212077e+04,
3.51712077e+04, 3.52212077e+04, 3.52712077e+04, 3.53212077e+04,
3.53712077e+04, 3.54212077e+04, 3.54712077e+04, 3.55212077e+04,
3.55712077e+04, 3.56212077e+04, 3.56712077e+04, 3.57212077e+04,
3.57712077e+04, 3.58212077e+04, 3.58712077e+04, 3.59212077e+04,
3.59712077e+04, 3.60212077e+04, 3.60712077e+04, 3.61212077e+04,
3.61712077e+04, 3.62212077e+04, 3.62712077e+04, 3.63212077e+04,
3.63712077e+04, 3.64212077e+04, 3.64712077e+04, 3.65212077e+04])
We can visualize the resulting rasterized image.
[8]:
#Plot unrasterized/rasterized images
fig,ax = plt.subplots(1,2)
ax[0].scatter(x,y,s=0.5,alpha=0.25)
ax[0].invert_yaxis()
ax[0].set_title('List of cells')
ax[0].set_aspect('equal')
W = W[0]
extent = (X_[0],X_[-1],Y_[0],Y_[-1])
ax[1].imshow(W, origin='lower')
ax[1].invert_yaxis()
ax[1].set_title('Rasterized')
# save figure
#fig.canvas.draw()
#fig.savefig(outname[:-4]+'_image.png')
[8]:
Text(0.5, 1.0, 'Rasterized')
Now we need to find an approximate slice number in the Allen Atlas to initialize the alignment. Evaluate whether the value of slice is similar to the target image by viewing a side by side comparison of the ABA slice and the target image. If not, change the value of slice. In this example, slice = 0 is anterior and slice = 264 is posterior.
[9]:
#find slice
#peruse through images in atlas
# Loading the atlas
slice = 140
vol,hdr = nrrd.read(imagefile)
A = vol
vol,hdr = nrrd.read(labelfile)
L = vol
dxA = np.diag(hdr['space directions'])
nxA = A.shape
xA = [np.arange(n)*d - (n-1)*d/2.0 for n,d in zip(nxA,dxA)]
XA = np.meshgrid(*xA,indexing='ij')
fig,ax = plt.subplots(1,2)
extentA = STalign.extent_from_x(xA[1:])
ax[0].imshow(A[slice],extent=extentA)
ax[0].set_title('Atlas Slice')
ax[1].imshow(W,extent=extentA)
ax[1].set_title('Target Image')
fig.canvas.draw()
[10]:
len(xA[2])
[10]:
228
Next, we need to find an approximate rotation angle of the Allen Atlas to initialize the alignment. Evaluate whether the rotation of the atlas is similar to the target image by viewing a side by side comparison of the ABA slice and the target image. If not, change the value of theta_deg. Note: the rotation here is defined in degrees.
[11]:
from scipy.ndimage import rotate
theta_deg = 90
fig,ax = plt.subplots(1,2)
extentA = STalign.extent_from_x(xA[1:])
ax[0].imshow(rotate(A[slice], angle=theta_deg),extent=extentA)
ax[0].set_title('Atlas Slice')
ax[1].imshow(W,extent=extentA)
ax[1].set_title('Target Image')
fig.canvas.draw()
Next, specify if there are any points on the image above that can be approximately matched with eachother and add them in the following format. This helps speed up alignment. (OPTIONAL)
[12]:
points_atlas = np.array([[0,0]]) #(-y,-x)
points_target = np.array([[-3700,0]]) # (-y,-x)
Li,Ti = STalign.L_T_from_points(points_atlas,points_target)
We can define atlas and target points, xI and xJ, as well as atlas and target images, I and J.
[13]:
xJ = [Y_,X_]
J = W[None]/np.mean(np.abs(W))
xI = xA
I = A[None] / np.mean(np.abs(A),keepdims=True)
I = np.concatenate((I,(I-np.mean(I))**2))
Inorm = STalign.normalize(I, t_min=0, t_max=1)
Jnorm = STalign.normalize(J, t_min=0, t_max=1)
Now that we can chosen the slice number and rotation angle, we are ready to initialize parameters related to our alignment.
[14]:
sigmaA = 0.1 #standard deviation of artifact intensities
sigmaB = 0.1 #standard deviation of background intensities
sigmaM = 0.1 #standard deviation of matching tissue intenities
muA = torch.tensor([0.7,0.7,0.7],device='cpu') #average of artifact intensities
muB = torch.tensor([0,0,0],device='cpu') #average of background intensities
We can change the parameters above by looking at the intensity histogram of our target image (below). We need to consider intensities of artifacts (tissue that is present in the target image and absent in the atlas), which is usually in the upper range of intensity values. We also need to find the intensity of the background values that does not correspond to any tissue, usually around 0, and standard deviations for these values. The matching tissue intensities are regions that should have tissue that can be aligned in both the atlas and the target image.
[15]:
fig,ax = plt.subplots()
ax.hist(Jnorm.ravel())
plt.xlabel('Intensity')
plt.ylabel('Number of Pixels')
plt.title('Intensity Histogram of Target Image')
[15]:
Text(0.5, 1.0, 'Intensity Histogram of Target Image')
The following parameters vary depending on your target image and the duration & accuracy of the alignment. The scale parameters refer to a scaling of the atlas in (x, y, z) to match the size of the target image. The nt parameter ensures smoothness and invertibility of the image transformation. The niter refers to the number of iteration steps in gradient descent to achieve minimum error between the altas and the target image.
[16]:
# initialize variables
scale_x = 4 #default = 0.9
scale_y = 4 #default = 0.9
scale_z = 0.9 #default = 0.9
theta0 = (np.pi/180)*theta_deg
# get an initial guess
if 'Ti' in locals():
T = np.array([-xI[0][slice],np.mean(xJ[0])-(Ti[0]*scale_y),np.mean(xJ[1])-(Ti[1]*scale_x)])
else:
T = np.array([-xI[0][slice],np.mean(xJ[0]),np.mean(xJ[1])])
#T = np.array([-xI[0][slice],0,0])
scale_atlas = np.array([[scale_z,0,0],
[0,scale_x,0],
[0,0,scale_y]])
L = np.array([[1.0,0.0,0.0],
[0.0,np.cos(theta0),-np.sin(theta0)],
[0.0,np.sin(theta0),np.cos(theta0)]])
L = np.matmul(L,scale_atlas)#np.identity(3)
Now, we can perform alignment of the atlas to our target slice.
[24]:
%%time
#returns mat = affine transform, v = velocity, xv = pixel locations of velocity points
transform = STalign.LDDMM_3D_to_slice(
xI,Inorm,xJ,Jnorm,
T=T,L=L,
nt=4,niter=800,
a=250,
device='cpu',
sigmaA = sigmaA, #standard deviation of artifact intensities
sigmaB = sigmaB, #standard deviation of background intensities
sigmaM = sigmaM, #standard deviation of matching tissue intenities
muA = muA, #average of artifact intensities
muB = muB #average of background intensities
)
/home/manju/.local/lib/python3.8/site-packages/STalign/STalign.py:1617: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
axE[2].set_yscale('log')
CPU times: user 2h 36min 26s, sys: 1h 16min 5s, total: 3h 52min 32s
Wall time: 13min 27s
[25]:
A = transform['A']
v = transform['v']
xv = transform['xv']
Xs = transform['Xs']
[26]:
print(A)
tensor([[ 1.2073e+00, 3.4987e-02, -1.3848e-02, -6.5878e+02],
[-6.2258e-02, 2.3100e-02, -4.8150e+00, 2.6193e+04],
[ 1.2410e-02, 4.6378e+00, -4.3118e-02, 1.8045e+04],
[ 0.0000e+00, 0.0000e+00, 0.0000e+00, 1.0000e+00]],
dtype=torch.float64)
Once we have aligned the atlas to our target image, we can project the brain regions defined in the atlas to our target image and export these annotations, along with information about the positions of the atlas that align with the slice.
[27]:
df = STalign.analyze3Dalign(labelfile, xv,v,A, xJ, dx, scale_x=scale_x, scale_y=scale_y,x=x,y=y, X_=X_, Y_=Y_, namesdict=namesdict,device='cpu')
/home/manju/.local/lib/python3.8/site-packages/STalign/STalign.py:1726: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).
A = torch.tensor(A)
/home/manju/.local/lib/python3.8/site-packages/STalign/STalign.py:1727: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).
v = torch.tensor(v)
/home/manju/.local/lib/python3.8/site-packages/STalign/STalign.py:1739: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).
XJ = torch.tensor(XJ)
[28]:
df
[28]:
| coord0 | coord1 | coord2 | x | y | struct_id | acronym | |
|---|---|---|---|---|---|---|---|
| 0 | 714.048218 | -3618.158516 | 5537.465920 | 10834.272727 | 1488.454545 | 0 | bg |
| 1 | 712.695551 | -3561.576081 | 5562.228887 | 12139.406780 | 864.788136 | 0 | bg |
| 2 | 712.907165 | -3574.675730 | 5547.685416 | 11848.132075 | 1225.679245 | 527 | AUDd1 |
| 3 | 713.804790 | -3607.279554 | 5543.688647 | 11090.980392 | 1342.539216 | 0 | bg |
| 4 | 713.301162 | -3589.913500 | 5543.613236 | 11517.905660 | 1351.735849 | 0 | bg |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 43336 | 667.020365 | -2502.125095 | 4416.027593 | 34741.963636 | 25621.704545 | 0 | bg |
| 43337 | 668.244750 | -2520.296934 | 4459.336172 | 34380.794118 | 24686.573529 | 0 | bg |
| 43338 | 668.459205 | -2534.351518 | 4446.841324 | 34106.068182 | 25049.780303 | 0 | bg |
| 43339 | 667.999837 | -2508.449947 | 4465.799706 | 34648.200000 | 24484.694118 | 0 | bg |
| 43340 | 667.127079 | -2502.288521 | 4422.370802 | 34753.990385 | 25492.134615 | 0 | bg |
43341 rows × 7 columns
Now, we can explore our alignment and brain region annotations!
We can visualize the MERFISH slice overlayed with the matched section in the Allen Brain Atlas.
[29]:
It = torch.tensor(I,device='cpu',dtype=torch.float64)
AI = STalign.interp3D(xI,It,Xs.permute(3,0,1,2),padding_mode="border")
Ishow_source = ((AI-torch.amin(AI,(1,2,3))[...,None,None])/(torch.amax(AI,(1,2,3))-torch.amin(AI,(1,2,3)))[...,None,None,None]).permute(1,2,3,0).clone().detach().cpu()
Jt = torch.tensor(J,device='cpu',dtype=torch.float64)
Ishow_target = Jt.permute(1,2,0).cpu()/torch.max(Jt).item()
import matplotlib as mpl
fig,ax = plt.subplots(1,3, figsize=(15,5))
ax0 = ax[0].imshow(Ishow_target, cmap = mpl.cm.Blues,alpha=0.9)
ax[0].set_title('STARmap Slice')
ax1 = ax[1].imshow(Ishow_source[0,:,:,0], cmap = mpl.cm.Reds,alpha=0.2)
ax[1].set_title('z=0 slice of Aligned 3D Allen Brain Atlas')
ax2 = ax[2].imshow(Ishow_target, cmap = mpl.cm.Blues,alpha=0.9)
ax2 = ax[2].imshow(Ishow_source[0,:,:,0], cmap = mpl.cm.Reds,alpha=0.3)
ax[2].set_title('Overlayed')
plt.show()
[30]:
from mpl_toolkits import mplot3d
x=df['coord0']
y=-df['coord1']
z=df['coord2']
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter3D(x,y,z,s=0.05)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
[30]:
Text(0.5, 0, 'z')
We can also plot all of the brain regions that are captured in the slice, as defined by our alignment.
[31]:
STalign.plot_brain_regions(df)
If we also interested in only a few regions, in this example let’s take VISp4 and VISp5, we can plot those regions using the following function.
[32]:
brain_regions = ['CA1']
STalign.plot_subset_brain_regions(df, brain_regions)
[ ]: