refFile
and preablationFile
.Eprofile.txt
and Eprofile.zip
, generated using Fiji plugin Zebra_ablate.py
.zStep
defined in the Fiji plugin.pixelSize
.codeDir
containing the repository files.FIJIexecutable
.
In [3]:
expDir = 'Y:/Nikita/Janelia_SPIM_selected-MVregistration/2015-07-15fish5dpf(H2B_ablation_RoL)/'
#Reference file, in which ablation ROIs were selected using Fiji plugin Zebra_ablate
refFile = expDir + 'dOMR0_20150715_143228/proc/ref.tif'
#Pre-ablation file name, presumably after some drift/rotation relative to the Reference file
preablationFile = expDir + 'preablation_20150715_151113/TM00000_CM0_CHN00.tif'
#existing ablation profile generated by Fiji plugin Zebra_ablate
ablationCoordFile = expDir + 'Eprofile.txt'
ablationCoordROIFile = ablationCoordFile + '.zip'
#Set image pixel size (x, y) in um
pixelSize = (0.406, 0.406)
In [4]:
regDir = expDir + 'MVregistration/'
ablationCoordFileRegistered = regDir + 'Eprofile_Corrected.txt'
ablationCoordROIFileRegistered = regDir + 'Eprofile_Corrected.txt.zip'
In [5]:
#Software folders
#codeDir = 'C:/Users/AhrensLab/Documents/GitHub/zebrascope_targets-master/MultiviewRegistration/'
#FIJIexecutable = 'C:/Users/AhrensLab/Desktop/Fiji.app.newest/ImageJ-win64.exe'
codeDir = 'C:/Users/nvladim/Documents/Nikita-Janelia/mapping_perturbation_data_ms/code/zebrascope_targets/MultiviewRegistration/'
FIJIexecutable = 'C:/Users/nvladim/Desktop/Fiji.app/ImageJ-win64.exe'
In [6]:
#Check if all files exist
import os
if not os.path.exists(expDir):
print("folder not found: ", expDir)
if not os.path.exists(refFile):
print("file not found: ", refFile)
if not os.path.exists(preablationFile):
print("folder not found: ", preablationFile)
if not os.path.exists(ablationCoordFile):
print("file not found: ", ablationCoordFile)
if not os.path.exists(ablationCoordROIFile):
print("file not found: ", ablationCoordROIFile)
if not os.path.exists(FIJIexecutable):
print("file not found: ", FIJIexecutable)
if not os.path.exists(codeDir):
print("folder not found: ", codeDir)
if not os.path.exists(regDir):
os.mkdir(regDir)
print("created output folder: ", regDir)
In [7]:
%%time
#Prepare files
from shutil import copyfile
import re
copyfile(codeDir + '/dataset_default.xml', regDir + '/dataset.xml') #copy XML file template, overwrite existing
#copy files to the registration folder
copyfile(refFile, regDir + 'stack0.tif')
copyfile(preablationFile, regDir + 'stack1.tif')
copyfile(ablationCoordFile, regDir + 'Eprofile_NoCorrection.txt')
copyfile(ablationCoordROIFile, regDir + 'Eprofile_NoCorrection.txt.zip')
#Modify the data path in `macro_detectPoints_register.ijm` file:
with open(codeDir + 'macro_detectPoints_register_template.ijm', 'rt') as f:
with open(codeDir + 'macro_detectPoints_register.ijm', 'wt') as f1:
for line in f:
if(re.search('[regDir]', line)):
line1 = line.replace('[regDir]/',regDir)
else:
line1 = line
f1.writelines(line1)
In [8]:
%%time
#Run the `Multiview Plugin` in batch mode:
#1. Find Interest Points in both images
#2. Register second stack relative to the first stack.
import subprocess
args = FIJIexecutable + ' -macro ' + codeDir + 'macro_detectPoints_register.ijm'
logStr = subprocess.check_output(args)
In [15]:
#Load the transformation matrix and apply to the ablation coordinates:
import xml.etree.ElementTree as ET
import numpy as np
import re
np.set_printoptions(precision=2, suppress=True)
tree = ET.parse(regDir + 'dataset.xml')
root = tree.getroot()
transformList = []
print('Found the following tranforms in ' + regDir + 'dataset.xml' ':\n')
for registration in root.iter('ViewRegistration'):
if registration.get('setup')=="1":
for transform in registration.iter('ViewTransform'):
name = transform.find('Name').text
print('Transform: ' + name + '\n')
#print transform.find('affine').text + '\n'
transformList.append(transform.find('affine').text)
#Parse the XML file into transform matrices, starting with `calibration` first
TransformMatrices = np.zeros((len(transformList),3,4))
for i in range(len(transformList)):
TransformMatrices[i,:,:] = np.fromstring(transformList[-i-1], sep = ' ').reshape(3,4)
print('Transforms have been parsed')
# Combine all registration transformations into single matrix
TransformMatrixAll = np.identity(4)
for i in range(1,len(transformList)):
TransformMatrixAll = np.dot(TransformMatrixAll,np.vstack((TransformMatrices[i,:,:], np.array([0,0,0,1]))))
print('Transforms have been combined\n')
# Inverse transf. matrix, to tranform stack0.tif ablation coordinates to stack1.tif space
TransformMatrixInv = np.linalg.inv(TransformMatrixAll)
print('Transform matrix has been inverted (px):\n')
print(TransformMatrixInv)
print('\n')
print('Translation terms (um):')
print('Drift along X: ' + '{:.1f}'.format(TransformMatrixInv[0,3]*pixelSize[0]))
print('Drift along Y: ' + '{:.1f}'.format(TransformMatrixInv[1,3]*pixelSize[0]))
print('Drift along Z: ' + '{:.1f}'.format(TransformMatrixInv[2,3]*pixelSize[0]))
print('\n')
#Parse ablation coordinates
xyz1List = []
with open(ablationCoordFile, 'rt') as f:
for line in f:
if(re.search('COORDINATES(.+)', line)):
coords = re.findall('\d+', line)
xyz1List.append((float(coords[0]), float(coords[1]), float(coords[2]), 1))
xyz1array = np.asarray(xyz1List)
print('Original ablation points have been parsed')
# convert Z-coordinates (um) into the same units as XY voxels
xyz1array_Zscaled = xyz1array.copy()
xyz1array_Zscaled[:,2] = xyz1array[:,2]/pixelSize[0]
#Apply transforms
xyzArrayRegistered_Zscaled = np.zeros(xyz1array.shape)
for i in range(xyz1array.shape[0]):
xyzArrayRegistered_Zscaled[i,:] = np.dot(TransformMatrixInv,xyz1array_Zscaled[i,:])
print('Transform has been applied')
# Convert Z-coordinates (px) back into microns
xyzArrayRegistered = xyzArrayRegistered_Zscaled.copy()
xyzArrayRegistered[:,2] = xyzArrayRegistered_Zscaled[:,2]*pixelSize[0]
# Write corrected ablation profile
with open(ablationCoordFile, 'rt') as f:
with open(ablationCoordFileRegistered, 'wt') as f1:
i_cell = 0
for line in f:
if(re.search('COORDINATES(.+)', line)):
line1 = ['COORDINATES('+str(np.round(xyzArrayRegistered[i_cell,0]))+',' \
+str(np.round(xyzArrayRegistered[i_cell,1]))+',' \
+str(np.round(xyzArrayRegistered[i_cell,2]))+')'+'\n']
i_cell+=1
else:
line1 = line
f1.writelines(line1)
print('\nCorrected ablation coordinate file is written into:\n' + ablationCoordFileRegistered)
#Create ZIP file with corrected rois
#Default values:
zStart = None
zStep = None
xOffset = None
yOffset = None
radiusROI = None
#Read actual values from e-profile
with open(ablationCoordFile, 'rt') as f:
for line in f:
if(re.search('<zStart', line)):
strVar = re.findall('\d+\.\d+', line)
zStart = float(strVar[0])
if(re.search('<zStep', line)):
strVar = re.findall('\d+\.\d+', line)
zStep = float(strVar[0])
elif(re.search('<xOffset', line)):
strVar = re.findall('\d+\.\d+', line)
xOffset = float(strVar[0])
elif(re.search('<yOffset', line)):
strVar = re.findall('\d+\.\d+', line)
yOffset = float(strVar[0])
elif(re.search('<radiusROI', line)):
strVar = re.findall('\d+\.\d+', line)
radiusROI = float(strVar[0])
print('\nParameters parsed from the coordinate file:')
print('zStart(um) = ' + str(zStart))
print('zStep(um) = ' + str(zStep))
print('xOffset(px) = ' + str(xOffset))
print('yOffset(px) = ' + str(yOffset))
# Compress all .roi files into one ZIP archine, readable by FIJI
import sys, os
sys.path.append(codeDir + '/PymageJ-devel')
from pymagej.roi import ROIEncoder, ROIRect, ROIOval
import zipfile
zf = zipfile.ZipFile(ablationCoordROIFileRegistered, mode = "w", compression=zipfile.ZIP_DEFLATED)
for i in range(xyzArrayRegistered.shape[0]):
roi_obj = ROIOval(int(-yOffset + xyzArrayRegistered[i,1]), \
int(-xOffset + xyzArrayRegistered[i,0]), \
int(-yOffset + xyzArrayRegistered[i,1] + 2*radiusROI), \
int(-xOffset + xyzArrayRegistered[i,0] + 2*radiusROI))
z_plane = int(np.round(xyzArrayRegistered[i,2]/zStep)) + 1
roi_obj.header = {'POSITION': z_plane}
ROIfileName = regDir + '/z' + str(z_plane) + 'cell' + str(i) + '.roi'
with ROIEncoder(ROIfileName, roi_obj) as roi:
roi.write()
zf.write(ROIfileName, 'z' + str(z_plane) + 'cell' + str(i) + '.roi')
os.remove(ROIfileName)
zf.close()
print('\nCorrected ZIP file with Fiji ROIs is written into:\n' + ablationCoordROIFileRegistered)
In [ ]:
In [ ]: