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)


    

    




Wall time: 22.6 s

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)


    

    




Wall time: 43 s

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)


    

    




Found the following tranforms in Y:/Nikita/Janelia_SPIM_selected-MVregistration/2015-07-15fish5dpf(H2B_ablation_RoL)/MVregistration/dataset.xml:

Transform: Redundant geometric local descriptor matching (translation invariant), RigidModel3D on [cells (c=0)]

Transform: calibration

Transforms have been parsed
Transforms have been combined

Transform matrix has been inverted (px):

[[ 1.   -0.    0.   -3.93]
 [ 0.    1.   -0.   -1.81]
 [-0.    0.    1.    5.16]
 [ 0.    0.    0.    1.  ]]


Translation terms (um):
Drift along X: -1.6
Drift along Y: -0.7
Drift along Z: 2.1


Original ablation points have been parsed
Transform has been applied

Corrected ablation coordinate file is written into:
Y:/Nikita/Janelia_SPIM_selected-MVregistration/2015-07-15fish5dpf(H2B_ablation_RoL)/MVregistration/Eprofile_Corrected.txt

Parameters parsed from the coordinate file:
zStart(um) = 0.0
zStep(um) = 5.0
xOffset(px) = 0.0
yOffset(px) = 512.0

Corrected ZIP file with Fiji ROIs is written into:
Y:/Nikita/Janelia_SPIM_selected-MVregistration/2015-07-15fish5dpf(H2B_ablation_RoL)/MVregistration/Eprofile_Corrected.txt.zip

In [ ]:

    




    

In [ ]: