Correct only z drift with Correct 3D

registration
Tags: #<Tag:0x00007fd6a17aef98>

#1

I have been using Correct 3D Drift with great results to process multi channel XYZT datasets.
Does anyone know if I can use it to correct only the drift in z without changing xy?
If not possible with Correct 3D Drift, is there another solution out there?

Thank you in advance for the help

Rafael


#2

Dear @Buono,

What exactly do you mean by that? Translating the z-slices along the z-axis over time while keeping xy constant? If so, would such a translation move individual slices or the complete stack for a time point?

Best,
Stefan


#3

Dear @stelfrich,

Thanks for the answer!

Yes

I would like to translate the complete stack for a time point to compensate for a physical drift in z. That z drift is a product of my sample moving away from the original focal point between time points.

The current challenge that I am trying to tackle is the following:
I work with plant roots. They literally grow away of the field of view during a time series.
I have dual channel 3D time series. Channel 2 is great for registration of xy translation that occurs due to growth of my samples. Channel 1 is great for seeing drifts in z.

My reasoning is to first register based on channel 2 and then register again, but just in z this time, based on channel 1.

Best,

Rafael


#4

How about reslicing the data to show xz by y (third dimension). In theory, if the data is aligned really nicely in xy, the xz images should just need registering with a translation transformation and then reslice the stack back to xy by z.

Whatever you are using to align in xy, just use it on the resliced data also.

What do you think?

Rob


#5

I would never have thought of that, but it should work!

The reason I was asking so specifically what @Buono wants to achieve is that I was sceptical about ImageJ’s ability to properly handle transformed stacks in a time series: ImageJ1 only supports fixed distances in z (also in xy) for the whole series. Hence, I don’t know how well the reslicing works and what the output will look like in that case. How many slices will be in the series afterwards, i.e. will a transformation be applied to the whole stack at a time point or individual slices?

Have you tried your suggestion, @7rebor? I’d really be interested in the outcome…


#6

I’m not entirely sure what you mean by “fixed distances”, do you mean fixed distance of transformations or slice depths/pixel sizes? Are you saying that an individual stack in a time series couldn’t be entirely transformed relative to the other stacks?

The 3D time series I deal with are over months and only have 9 timepoints, so I usually keep the timepoints separate and use max projections (for example) to align them to each other in xy (using MultiStackReg and using the transform files to apply to the entire stack). Separately I manually choose the correct z slice for each timepoint and crop them using a script comparing stack sizes and positions of reference slices to finally concatenate them to a 4D stack. Obviously the manual method is not possible if you have 100 timepoints per series!

I did have a go with my example (above). What I did was use MultiStackReg to align the two xzy stacks to each other, but on second thoughts I would use the max projections of those stacks to apply a transform to the whole stack to speed it up (although this might not work, I’ll have a go now and see what happens, could be useful for me - eliminates my manual step).

EDIT: Neither of them really work for my data, I guess the problem is that each z slice is reduced to 1 pixel in the resliced version and the registration isn’t sensitive enough to register it perfectly. It was worth a try!


#7

I mean fixed “slice depth” (although that assumes there is an extent of a sample point, which there isn’t). Imagine the following situation (sorry for the sloppy artwork):

How will the information of the translation (along z-axis) be encoded in the time series? This only gets worth if the “slice depth” were changed by the translation. I fear I am missing something obvious here…


#8

Nice artwork!

I see what you mean, the best possibility is correcting the drift by the value of the slice depth then? So if it drifts by 1 um down, you can shift all the slice positions up by one for that timepoint. An easy way to do this is to remove/add slices as necessary from the start/end of the individual timepoints to change the slice values in that stack relative to the other timepoints. This is what I do for my small time series, where I don’t need very exact drift correction (just within 1 um, which is a single slice for me).

An excerpt from my script:

imgList = getList("image.titles"); // list of all the headers of open image windows

slice = newArray(9);

selectWindow(imgList[0]); // selects first window
Stack.getDimensions(w, h, channelCount, sliceCount, frameCount); // gets dimensions of the stack
run("Tile"); // tiles the open image windows
waitForUser("Select the same slice in all stacks"); // user selects the closest slice in subsequent timepoints to the reference slice
selectWindow(imgList[0]); // selects first window
Stack.getPosition(channel, slice[0], frame); // gets slice user has chosen

for (k=1; k < imgList.length; k++) { // for all but the first image in the image list
		selectWindow(imgList[k]);
		Stack.getPosition(channel, slice[k], frame); // gets slice user has chosen
		if (slice[k] < slice[0]) { // if it is closer to the start of the stack additional slices are added to make the slice at the same slice# as the reference stack
			sliceAdditions = (slice[0] - slice[k]); // # additions is calculated by the diff. between ref and current stack
			selectWindow(imgList[k]);
			for (n=0; n < sliceAdditions; n++) {
			    Stack.setSlice(1); // set the current slice to be the first slice
				run("Add Slice", "add=slice prepend"); // add the slice before the current position (add to beginning of stack)
			}
		}
		if (slice[k] > slice[0]) { // if the selected slice is closer to the end than the reference
			sliceDeletions = (slice[k] - slice[0]); // #slice to delete from the start to make the slice in same position as ref slice
			selectWindow(imgList[k]);
			for (n=0; n < sliceDeletions; n++) {
				Stack.setSlice(1);
				run("Delete Slice", "delete=slice");
			}
		}
	}

After that I also match the stack sizes by adding empty slices to the start/end.

Something else I thought of a while ago is to choose a slice in the stack and then using the SNR plugin from EPFL to work out the same slice in the next timepoint that corresponds to that slice. This would be done by taking the highest SNR from comparisons (this assumes the slices are aligned perfectly in xy). Then repeat sequentially through the timeseries to keep the slice in the same position in each timeseries and essentially adding/removing slices in the stacks to align the entire time series in z. Although I don’t know how to read out the table values to do this.

P.S. I tried correct 3D drift and it seems to be unreliable compared to my approach (which does most of my stacks almost perfectly), the only problem is that my approach requires some manual input for the z alignment. Been wondering how to automate that for a long time (this is why I’m interested in this thread).


#9

I think it is the only reasonable option.


Maybe @Christian_Tischer (who has taken over Lead on Correct 3D Drift) input on that topic as well.

Also, @Buono, you can give the workflow

  1. Image > Stacks > Reslice
  2. Correct 3D drift
  3. Image > Stacks > Reslice

a try.


#10

@stelfrich

I was just trying it. That is a very nice idea!

It did a nice job at correcting z-drift, but in the process it did change the previously performed xy correction. My guess is that it is due to the fact that the sample (a root tip) is growing longer and longer over the time points. So there is always a portion of the image that will have moving features, no matter how much I register it. I am adding an example below of what I mean:

Unregistered stack:

Registered using Correct 3D Drift and an ROI at the TIP region to get that region steady. In a perfect situation this would mean that the root would grow away of the TIP on the registered stack.

I imagine that Correct 3D Drift tries to correct that moving portion at the left of the image. I have tried asking Correct 3D Drift to only take into account the somewhat steady region of the tip on the Reslice-Correct-Reslice strategy, but it also did not work well.

My original idea was to figure out how to get Correct 3D Drift to not apply whatever corrections it finds for xy and only apply the z corrections. In this scenario I would do a first round of registration that should take care of xy and then a second one that would only take care of z. I tried going through the code [here] (https://github.com/fiji/Correct_3D_Drift/blob/master/src/main/resources/scripts/Plugins/Registration/Correct_3D_drift.py) but could not find how to do it.


#11

Is it just stacks of brightfield/widefield images that you are registering, or also confocal images?

I am asking because the extended depth of field plugin will find the most in focus parts of the image and use that to make a projection of the stack. This can be applied to each timepoint using some scripting. It might solve your problem if you’re not using confocal microscopy and therefore only interested in one z-plane.

Correct 3D Drift could still be used for xy correction, but with the projections from extended DoF plugin.


#12

@7rebor

The brightfield was mostly for example. The data is all confocal. I tend to collect brightfield just as a reference.
The current datasets are nuclei in channel 1 and brightfield (T-PMT) on channel 2. I am using the brightfield to register xy. But then I need the nuclei (channel 1) to stop drifting up and down in z.


#13

Since my last reply I tried the following strategy:

-I got the Correct 3D Drift code from here
-Opened it in FIJI and started deleting/substituting with zeroes whatever looked to me like a line that was calculating or applying xy shifts.
-Since I did not really know what I was doing, that took a few tries
-Ended up with a few “modifications” to the original script that mostly do what my current dataset needs (translation only in z)
This modified version has modifications in lines 77, 262, and 476
Below is the modified version in its entirety:

# #%L
# Script to register time frames (stacks) to each other.
# %%
# Copyright (C) 2010 - 2016 Fiji development team
# %%
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public
# License along with this program.  If not, see
# <http://www.gnu.org/licenses/gpl-3.0.html>.
# #L%
###

# Robert Bryson-Richardson and Albert Cardona 2010-10-08 at Estoril, Portugal
# EMBO Developmental Imaging course by Gabriel Martins
#
# Register time frames (stacks) to each other using Stitching_3D library
# to compute translations only, in all 3 spatial axes.
# Operates on a virtual stack.
# 23/1/13 -
# added user dialog to make use of virtual stack an option
# 10/01/16 -
# Christian Tischer (tischitischer@gmail.com)
# major changes and additions:
# - it now also works for 2D time-series (used to be 3D only)
# - option: measure drift on multiple timescales (this allows to also find slow drift components of less than 1 pixel per frame)
# - option: correct sub-pixel drift computing the shifted images using TransformJ
# - option: if a ROI is put on the image, only this part of the image is considered for drift computation
#           the ROI is moved along with the detected drift thereby tracking the structure of interest
# - macro recording is compatible with previous version
# 21/11/16
# - fixed a bug related to hyperstack conversion

from ij import VirtualStack, IJ, CompositeImage, ImageStack, ImagePlus
from ij.process import ColorProcessor
from ij.plugin import HyperStackConverter
from ij.io import DirectoryChooser, FileSaver, SaveDialog
from ij.gui import GenericDialog, YesNoCancelDialog, Roi
from mpicbg.imglib.image import ImagePlusAdapter
from mpicbg.imglib.algorithm.fft import PhaseCorrelation
from org.scijava.vecmath import Point3i  #from javax.vecmath import Point3i # java6
from org.scijava.vecmath import Point3f  #from javax.vecmath import Point3f # java6
from java.io import File, FilenameFilter
from java.lang import Integer
import math, os, os.path

# sub-pixel translation using imglib2
from net.imagej.axis import Axes
from net.imglib2.img.display.imagej import ImageJFunctions
from net.imglib2.realtransform import RealViews, Translation3D, Translation2D
from net.imglib2.view import Views
from net.imglib2.img.imageplus import ImagePlusImgs
from net.imglib2.converter import Converters
from net.imglib2.converter.readwrite import RealFloatSamplerConverter
from net.imglib2.interpolation.randomaccess import NLinearInterpolatorFactory

def translate_single_stack_using_imglib2(imp, dx, dy, dz):
  # wrap into a float imglib2 and translate
  #   conversion into float is necessary due to "overflow of n-linear interpolation due to accuracy limits of unsigned bytes"
  #   see: https://github.com/fiji/fiji/issues/136#issuecomment-173831951
  img = ImagePlusImgs.from(imp.duplicate())
  extended = Views.extendZero(img)
  converted = Converters.convert(extended, RealFloatSamplerConverter())
  interpolant = Views.interpolate(converted, NLinearInterpolatorFactory())
  
  # translate
  if imp.getNDimensions()==3:
    transformed = RealViews.affine(interpolant, Translation3D(0, 0, dz)) ###changed dx and dy to zero
  elif imp.getNDimensions()==2:
    transformed = RealViews.affine(interpolant, Translation2D(dx, dy))
  else:
    IJ.log("Can only work on 2D or 3D stacks")
    return None
  
  cropped = Views.interval(transformed, img)
  # wrap back into bit depth of input image and return
  bd = imp.getBitDepth()
  if bd==8:
    return(ImageJFunctions.wrapUnsignedByte(cropped,"imglib2"))
  elif bd == 16:
    return(ImageJFunctions.wrapUnsignedShort(cropped,"imglib2"))
  elif bd == 32:
    return(ImageJFunctions.wrapFloat(cropped,"imglib2"))
  else:
    return None    

'''
def translate_single_stack_using_imagescience(imp, dx, dy, dz):
  translator = Translate()
  output = translator.run(Image.wrap(imp), dx, dy, dz, Translate.LINEAR)
  return output.imageplus()
'''

def compute_stitch(imp1, imp2):
  """ Compute a Point3i that expressed the translation of imp2 relative to imp1."""
  phc = PhaseCorrelation(ImagePlusAdapter.wrap(imp1), ImagePlusAdapter.wrap(imp2), 5, True)
  phc.process()
  p = phc.getShift().getPosition()
  if len(p)==3: # 3D data
    p3 = p
  elif len(p)==2: # 2D data: add zero shift
    p3 = [p[0],p[1],0]
  return Point3i(p3)

def extract_frame(imp, frame, channel, z_min, z_max):
  """ From a VirtualStack that is a hyperstack, contained in imp,
  extract the timepoint frame as an ImageStack, and return it.
  It will do so only for the given channel. """
  stack = imp.getStack() # multi-time point virtual stack
  stack2 = ImageStack(imp.width, imp.height, None)
  for s in range(int(z_min), int(z_max)+1):
    i = imp.getStackIndex(channel, s, frame)  
    stack2.addSlice(str(s), stack.getProcessor(i))
  return stack2


def extract_frame_process_roi(imp, frame, channel, process, background, roi, z_min, z_max):
  # extract frame and channel 
  imp_frame = ImagePlus("", extract_frame(imp, frame, channel, z_min, z_max)).duplicate()
  # check for roi and crop
  if roi != None:
    #print roi.getBounds()
    imp_frame.setRoi(roi)
    IJ.run(imp_frame, "Crop", "")
  # subtract background  
  if background > 0:
    #IJ.log("Subtracting "+str(background));
    IJ.run(imp_frame, "Subtract...", "value="+str(background)+" stack");
  # enhance edges  
  if process:
    IJ.run(imp_frame, "Mean 3D...", "x=1 y=1 z=0");
    IJ.run(imp_frame, "Find Edges", "stack");

  # return
  return imp_frame

def add_Point3f(p1, p2):
  p3 = Point3f(0,0,0)
  p3.x = p1.x + p2.x
  p3.y = p1.y + p2.y
  p3.z = p1.z + p2.z
  return p3

def subtract_Point3f(p1, p2):
  p3 = Point3f(0,0,0)
  p3.x = p1.x - p2.x
  p3.y = p1.y - p2.y
  p3.z = p1.z - p2.z
  return p3


def shift_between_rois(roi2, roi1):
  """ computes the relative xy shift between two rois 
  """ 
  dr = Point3f(0,0,0)
  dr.x = roi2.getBounds().x - roi1.getBounds().x
  dr.y = roi2.getBounds().y - roi1.getBounds().y
  dr.z = 0
  return dr

def shift_roi(imp, roi, dr):
  """ shifts a roi in x,y by dr.x and dr.y
  if the shift would cause the roi to be outside the imp,
  it only shifts as much as possible maintaining the width and height
  of the input roi
  """ 
  if roi == None:
    return roi
  else:
    r = roi.getBounds()
    # init x,y coordinates of new shifted roi
    sx = 0
    sy = 0
    # x shift
    if (r.x + dr.x) < 0:
      sx = 0
    elif (r.x + dr.x + r.width) > imp.width: 
      sx = int(imp.width-r.width)
    else:
      sx = r.x + int(dr.x)
    # y shift
    if (r.y + dr.y) < 0:
      sy = 0
    elif (r.y + dr.y + r.height) > imp.height: 
      sy = int(imp.height-r.height)
    else:
      sy = r.y + int(dr.y)
    # return shifted roi
    shifted_roi = Roi(sx, sy, r.width, r.height)
    return shifted_roi   
  
def compute_and_update_frame_translations_dt(imp, channel, dt, process, background, z_min, z_max, shifts = None):
  """ imp contains a hyper virtual stack, and we want to compute
  the X,Y,Z translation between every t and t+dt time points in it
  using the given preferred channel. 
  if shifts were already determined at other (lower) dt 
  they will be used and updated.
  """
  nt = imp.getNFrames()
  # get roi (could be None)
  roi = imp.getRoi()
  #if roi:
  #  print "ROI is at", roi.getBounds()   
  # init shifts
  if shifts == None:
    shifts = []
    for t in range(nt):
      shifts.append(Point3f(0,0,0))
  # compute shifts
  IJ.showProgress(0)
  for t in range(dt, nt+dt, dt):
    if t > nt-1: # together with above range till nt+dt this ensures that the last data points are not missed out
      t = nt-1 # nt-1 is the last shift (0-based)
    IJ.log("      between frames "+str(t-dt+1)+" and "+str(t+1))      
    # get (cropped and processed) image at t-dt
    roi1 = shift_roi(imp, roi, shifts[t-dt])
    imp1 = extract_frame_process_roi(imp, t+1-dt, channel, process, background, roi1, z_min, z_max)
    # get (cropped and processed) image at t-dt
    roi2 = shift_roi(imp, roi, shifts[t])
    imp2 = extract_frame_process_roi(imp, t+1, channel, process, background, roi2, z_min, z_max)
    #if roi:
    #  print "ROI at frame",t-dt+1,"is",roi1.getBounds()   
    #  print "ROI at frame",t+1,"is",roi2.getBounds()   
    # compute shift
    local_new_shift = compute_stitch(imp2, imp1)
    if roi: # total shift is shift of rois plus measured drift
      #print "correcting measured drift of",local_new_shift,"for roi shift:",shift_between_rois(roi2, roi1)
      local_new_shift = add_Point3f(local_new_shift, shift_between_rois(roi2, roi1))
    # determine the shift that we knew alrady
    local_shift = subtract_Point3f(shifts[t],shifts[t-dt])
    # compute difference between new and old measurement (which come from different dt)   
    add_shift = subtract_Point3f(local_new_shift,local_shift)
    #print "++ old shift between %s and %s: dx=%s, dy=%s, dz=%s" % (int(t-dt+1),int(t+1),local_shift.x,local_shift.y,local_shift.z)
    #print "++ add shift between %s and %s: dx=%s, dy=%s, dz=%s" % (int(t-dt+1),int(t+1),add_shift.x,add_shift.y,add_shift.z)
    # update shifts from t-dt to the end (assuming that the measured local shift will presist till the end)
    for i,tt in enumerate(range(t-dt,nt)):
      # for i>dt below expression basically is a linear drift predicition for the frames at tt>t
      # this is only important for predicting the best shift of the ROI 
      # the drifts for i>dt will be corrected by the next measurements
      shifts[tt].x += 1.0*i/dt * add_shift.x
      shifts[tt].y += 1.0*i/dt * add_shift.y
      shifts[tt].z += 1.0*i/dt * add_shift.z
      #print "updated shift till frame",tt+1,"is",shifts[tt].x,shifts[tt].y,shifts[tt].z
    IJ.showProgress(1.0*t/(nt+1))
  
  IJ.showProgress(1)
  return shifts


def convert_shifts_to_integer(shifts):
  int_shifts = []
  for shift in shifts: 
    int_shifts.append(Point3i(int(round(0)),int(round(0)),int(round(shift.z)))) ###changed x and y to zero
  return int_shifts

def compute_min_max(shifts):
  """ Find out the top left up corner, and the right bottom down corner,
  namely the bounds of the new virtual stack to create.
  Expects absolute shifts. """
  minx = Integer.MAX_VALUE
  miny = Integer.MAX_VALUE
  minz = Integer.MAX_VALUE
  maxx = -Integer.MAX_VALUE
  maxy = -Integer.MAX_VALUE
  maxz = -Integer.MAX_VALUE
  for shift in shifts:
    minx = min(minx, shift.x)
    miny = min(miny, shift.y)
    minz = min(minz, shift.z)
    maxx = max(maxx, shift.x)
    maxy = max(maxy, shift.y)
    maxz = max(maxz, shift.z)  
  return minx, miny, minz, maxx, maxy, maxz

def zero_pad(num, digits):
  """ for 34, 4 --> '0034' """
  str_num = str(num)
  while (len(str_num) < digits):
    str_num = '0' + str_num
  return str_num

def invert_shifts(shifts):
  """ invert shifts such that they can be used for correction.
  """
  for shift in shifts:
    shift.x *= -1
    shift.y *= -1
    shift.z *= -1
  return shifts


def register_hyperstack(imp, channel, shifts, target_folder, virtual):
  """ Takes the imp, determines the x,y,z drift for each pair of time points, using the preferred given channel,
  and outputs as a hyperstack."""
  # Compute bounds of the new volume,
  # which accounts for all translations:
  minx, miny, minz, maxx, maxy, maxz = compute_min_max(shifts)
  # Make shifts relative to new canvas dimensions
  # so that the min values become 0,0,0
  for shift in shifts:
    shift.x -= minx
    shift.y -= miny
    shift.z -= minz
  #print "shifts relative to new dimensions:"
  #for s in shifts:
  #  print s.x, s.y, s.z
  # new canvas dimensions:r
  width = imp.width + maxx - minx
  height = maxy - miny + imp.height
  slices = maxz - minz + imp.getNSlices()

  print "New dimensions:", width, height, slices
  # Prepare empty slice to pad in Z when necessary
  empty = imp.getProcessor().createProcessor(width, height)

  # if it's RGB, fill the empty slice with blackness
  if isinstance(empty, ColorProcessor):
    empty.setValue(0)
    empty.fill()
  # Write all slices to files:
  stack = imp.getStack()

  if virtual is False:
  	registeredstack = ImageStack(width, height, imp.getProcessor().getColorModel())
  names = []
  
  for frame in range(1, imp.getNFrames()+1):
 
    shift = shifts[frame-1]
    
    #print "frame",frame,"correcting drift",-shift.x-minx,-shift.y-miny,-shift.z-minz
    IJ.log("    frame "+str(frame)+" correcting drift "+str(-shift.x-minx)+","+str(-shift.y-miny)+","+str(-shift.z-minz))
    
    fr = "t" + zero_pad(frame, len(str(imp.getNFrames())))
    # Pad with empty slices before reaching the first slice
    for s in range(shift.z):
      ss = "_z" + zero_pad(s + 1, len(str(slices))) # slices start at 1
      for ch in range(1, imp.getNChannels()+1):
        name = fr + ss + "_c" + zero_pad(ch, len(str(imp.getNChannels()))) +".tif"
        names.append(name)

        if virtual is True:
          currentslice = ImagePlus("", empty)
          currentslice.setCalibration(imp.getCalibration().copy())
          currentslice.setProperty("Info", imp.getProperty("Info"))
          FileSaver(currentslice).saveAsTiff(target_folder + "/" + name)
        else:
          empty = imp.getProcessor().createProcessor(width, height)
          registeredstack.addSlice(str(name), empty)
    
    
    # Add all proper slices
    stack = imp.getStack()
    for s in range(1, imp.getNSlices()+1):
      ss = "_z" + zero_pad(s + shift.z, len(str(slices)))
      for ch in range(1, imp.getNChannels()+1):
         ip = stack.getProcessor(imp.getStackIndex(ch, s, frame))
         ip2 = ip.createProcessor(width, height) # potentially larger
         ip2.insert(ip, shift.x, shift.y)
         name = fr + ss + "_c" + zero_pad(ch, len(str(imp.getNChannels()))) +".tif"
         names.append(name)

         if virtual is True:
           currentslice = ImagePlus("", ip2)
           currentslice.setCalibration(imp.getCalibration().copy())
           currentslice.setProperty("Info", imp.getProperty("Info"));
           FileSaver(currentslice).saveAsTiff(target_folder + "/" + name)
         else:
           registeredstack.addSlice(str(name), ip2)

    # Pad the end
    for s in range(shift.z + imp.getNSlices(), slices):
      ss = "_z" + zero_pad(s + 1, len(str(slices)))
      for ch in range(1, imp.getNChannels()+1):
        name = fr + ss + "_c" + zero_pad(ch, len(str(imp.getNChannels()))) +".tif"
        names.append(name)

        if virtual is True:
          currentslice = ImagePlus("", empty)
          currentslice.setCalibration(imp.getCalibration().copy())
          currentslice.setProperty("Info", imp.getProperty("Info"))
          FileSaver(currentslice).saveAsTiff(target_folder + "/" + name)
        else:
          registeredstack.addSlice(str(name), empty)
 
  if virtual is True:
    # Create virtual hyper stack
    registeredstack = VirtualStack(width, height, None, target_folder)
    for name in names:
      registeredstack.addSlice(name)
  
  registeredstack_imp = ImagePlus("registered time points", registeredstack)
  registeredstack_imp.setCalibration(imp.getCalibration().copy())
  registeredstack_imp.setProperty("Info", imp.getProperty("Info"))
  registeredstack_imp = HyperStackConverter.toHyperStack(registeredstack_imp, imp.getNChannels(), len(names) / (imp.getNChannels() * imp.getNFrames()), imp.getNFrames(), "xyczt", "Composite");    
  
  return registeredstack_imp

  
def register_hyperstack_subpixel(imp, channel, shifts, target_folder, virtual):
  """ Takes the imp, determines the x,y,z drift for each pair of time points, using the preferred given channel,
  and outputs as a hyperstack.
  The shifted image is computed using TransformJ allowing for sub-pixel shifts using interpolation.
  This is quite a bit slower than just shifting the image by full pixels as done in above function register_hyperstack().
  However it significantly improves the result by removing pixel jitter.
  """
  # Compute bounds of the new volume,
  # which accounts for all translations:
  minx, miny, minz, maxx, maxy, maxz = compute_min_max(shifts)
  # Make shifts relative to new canvas dimensions
  # so that the min values become 0,0,0
  for shift in shifts:
    shift.x -= minx
    shift.y -= miny
    shift.z -= minz
  # new canvas dimensions:
  width = int(imp.width + maxx - minx)
  height = int(maxy - miny + imp.height)
  slices = int(maxz - minz + imp.getNSlices())

  #print "New dimensions:", width, height, slices
    
  # prepare stack for final results
  stack = imp.getStack()
  if virtual is True: 
    names = []
  else:
    registeredstack = ImageStack(width, height, imp.getProcessor().getColorModel())
  
  # prepare empty slice for padding
  empty = imp.getProcessor().createProcessor(width, height)

  IJ.showProgress(0)

  # get raw data as stack
  stack = imp.getStack()

  # loop across frames
  for frame in range(1, imp.getNFrames()+1):
      
    IJ.showProgress(frame / float(imp.getNFrames()+1))
    fr = "t" + zero_pad(frame, len(str(imp.getNFrames()))) # for saving files in a virtual stack
    
    # get and report current shift
    shift = shifts[frame-1]
    #print "frame",frame,"correcting drift",-shift.x-minx,-shift.y-miny,-shift.z-minz
    IJ.log("    frame "+str(frame)+" correcting drift "+str(round(-shift.x-minx,2))+","+str(round(-shift.y-miny,2))+","+str(round(-shift.z-minz,2)))

    # loop across channels
    for ch in range(1, imp.getNChannels()+1):      
      
      tmpstack = ImageStack(width, height, imp.getProcessor().getColorModel())

      # get all slices of this channel and frame
      for s in range(1, imp.getNSlices()+1):
        ip = stack.getProcessor(imp.getStackIndex(ch, s, frame))
        ip2 = ip.createProcessor(width, height) # potentially larger
        ip2.insert(ip, 0, 0)
        tmpstack.addSlice("", ip2)

      # Pad the end (in z) of this channel and frame
      for s in range(imp.getNSlices(), slices):
        tmpstack.addSlice("", empty)

      # subpixel translation
      imp_tmpstack = ImagePlus("", tmpstack)
      imp_translated = translate_single_stack_using_imglib2(imp_tmpstack, 0, 0, shift.z) ###Changed shift.x and shift.y to zero
      
      # add translated stack to final time-series
      translated_stack = imp_translated.getStack()
      for s in range(1, translated_stack.getSize()+1):
        ss = "_z" + zero_pad(s, len(str(slices)))
        ip = translated_stack.getProcessor(s).duplicate() # duplicate is important as otherwise it will only be a reference that can change its content  
        if virtual is True:
          name = fr + ss + "_c" + zero_pad(ch, len(str(imp.getNChannels()))) +".tif"
          names.append(name)
          currentslice = ImagePlus("", ip)
          currentslice.setCalibration(imp.getCalibration().copy())
          currentslice.setProperty("Info", imp.getProperty("Info"));
          FileSaver(currentslice).saveAsTiff(target_folder + "/" + name)
        else:
          registeredstack.addSlice("", ip)    

  IJ.showProgress(1)
    
  if virtual is True:
    # Create virtual hyper stack
    registeredstack = VirtualStack(width, height, None, target_folder)
    for name in names:
      registeredstack.addSlice(name)
  
  registeredstack_imp = ImagePlus("registered time points", registeredstack)
  registeredstack_imp.setCalibration(imp.getCalibration().copy())
  registeredstack_imp.setProperty("Info", imp.getProperty("Info"))
  registeredstack_imp = HyperStackConverter.toHyperStack(registeredstack_imp, imp.getNChannels(), slices, imp.getNFrames(), "xyzct", "Composite");    
  
  return registeredstack_imp
  

class Filter(FilenameFilter):
  def accept(self, folder, name):
    return not File(folder.getAbsolutePath() + "/" + name).isHidden()

def validate(target_folder):
  f = File(target_folder)
  if len(File(target_folder).list(Filter())) > 0:
    yn = YesNoCancelDialog(IJ.getInstance(), "Warning!", "Target folder is not empty! May overwrite files! Continue?")
    if yn.yesPressed():
      return True
    else:
      return False
  return True

def getOptions(imp):
  gd = GenericDialog("Correct 3D Drift Options")
  channels = []
  for ch in range(1, imp.getNChannels()+1 ):
    channels.append(str(ch))
  gd.addChoice("Channel for registration:", channels, channels[0])
  gd.addCheckbox("Multi_time_scale computation for enhanced detection of slow drifts?", False)
  gd.addCheckbox("Sub_pixel drift correction (possibly needed for slow drifts)?", False)
  gd.addCheckbox("Edge_enhance images for possibly improved drift detection?", False)
  gd.addNumericField("Only consider pixels with values larger than:", 0, 0)
  gd.addNumericField("Lowest z plane to take into account:", 1, 0)
  gd.addNumericField("Highest z plane to take into account:", imp.getNSlices(), 0)
  gd.addCheckbox("Use virtualstack for saving the results to disk to save RAM?", False)
  gd.addCheckbox("Only compute drift vectors?", False)
  gd.addMessage("If you put a ROI, drift will only be computed in this region;\n the ROI will be moved along with the drift to follow your structure of interest.")
  gd.showDialog()
  if gd.wasCanceled():
    return
  channel = gd.getNextChoiceIndex() + 1  # zero-based
  multi_time_scale = gd.getNextBoolean()
  subpixel = gd.getNextBoolean()
  process = gd.getNextBoolean()
  background = gd.getNextNumber()
  z_min = gd.getNextNumber()
  z_max = gd.getNextNumber()
  virtual = gd.getNextBoolean()
  only_compute = gd.getNextBoolean()
  return channel, virtual, multi_time_scale, subpixel, process, background, z_min, z_max, only_compute

def save_shifts(shifts, roi):
  sd = SaveDialog('please select shift file for saving', 'shifts', '.txt')
  fp = os.path.join(sd.getDirectory(),sd.getFileName())
  f = open(fp, 'w')
  txt = []
  txt.append("ROI zero-based")
  txt.append("\nx_min\ty_min\tz_min\tx_max\ty_max\tz_max")
  txt.append("\n"+str(roi[0])+"\t"+str(roi[1])+"\t"+str(roi[2])+"\t"+str(roi[3])+"\t"+str(roi[4])+"\t"+str(roi[5]))
  txt.append("\nShifts")
  txt.append("\ndx\tdy\tdz")  
  for shift in shifts:
    txt.append("\n"+str(shift.x)+"\t"+str(shift.y)+"\t"+str(shift.z))
  f.writelines(txt)
  f.close()


def run():

  IJ.log("Correct_3D_Drift")
    
  imp = IJ.getImage()
  if imp is None:
    return
  if 1 == imp.getNFrames():
    IJ.showMessage("Cannot register because there is only one time frame.\nPlease check [Image > Properties...].")
    return

  options = getOptions(imp)
  if options is not None:
    channel, virtual, multi_time_scale, subpixel, process, background, z_min, z_max, only_compute = options
  else:
    return # user pressed Cancel

  if z_min < 1:
    IJ.showMessage("The minimal z plane must be >=1.")
    return
    
  if z_max > imp.getNSlices():
    IJ.showMessage("Your image only has "+str(imp.getNSlices())+" z-planes, please adapt your z-range.")
    return
  
  if virtual is True:
    dc = DirectoryChooser("Choose target folder to save image sequence")
    target_folder = dc.getDirectory()
    if target_folder is None:
      return # user canceled the dialog
    if not validate(target_folder):
      return
  else:
    target_folder = None 

  # compute shifts
  IJ.log("  computing drifts..."); #print("\nCOMPUTING SHIFTS:")

  IJ.log("    at frame shifts of 1"); 
  dt = 1; shifts = compute_and_update_frame_translations_dt(imp, channel, dt, process, background, z_min, z_max)
  
  # multi-time-scale computation
  if multi_time_scale is True:
    dt_max = imp.getNFrames()-1
    # computing drifts on exponentially increasing time scales 3^i up to 3^6
    # ..one could also do this with 2^i or 4^i
    # ..maybe make this a user choice? did not do this to keep it simple.
    dts = [3,9,27,81,243,729,dt_max] 
    for dt in dts:
      if dt < dt_max:
        IJ.log("    at frame shifts of "+str(dt)) 
        shifts = compute_and_update_frame_translations_dt(imp, channel, dt, process, background, z_min, z_max, shifts)
      else: 
        IJ.log("    at frame shifts of "+str(dt_max));
        shifts = compute_and_update_frame_translations_dt(imp, channel, dt_max, process, background, z_min, z_max, shifts)
        break

  # invert measured shifts to make them the correction
  shifts = invert_shifts(shifts)
  #print(shifts)
  
  
  # apply shifts
  if not only_compute:
    
    IJ.log("  applying shifts..."); #print("\nAPPLYING SHIFTS:")
    
    if subpixel:
      registered_imp = register_hyperstack_subpixel(imp, channel, shifts, target_folder, virtual)
    else:
      shifts = convert_shifts_to_integer(shifts)
      registered_imp = register_hyperstack(imp, channel, shifts, target_folder, virtual)
    
    
    if virtual is True:
      if 1 == imp.getNChannels():
        ip=imp.getProcessor()
        ip2=registered_imp.getProcessor()
        ip2.setColorModel(ip.getCurrentColorModel())
        registered_imp.show()
      else:
    	registered_imp.copyLuts(imp)
    	registered_imp.show()
    else:
      if 1 ==imp.getNChannels():
        registered_imp.show()
      else:
        registered_imp.copyLuts(imp)
        registered_imp.show()
  
    registered_imp.show()
  
  else:
   
    if imp.getRoi(): 
      xmin = imp.getRoi().getBounds().x
      ymin = imp.getRoi().getBounds().y
      zmin = 0
      xmax = xmin + imp.getRoi().getBounds().width - 1
      ymax = ymin + imp.getRoi().getBounds().height - 1
      zmax = imp.getNSlices()-1  
    else:
      xmin = 0
      ymin = 0
      zmin = 0
      xmax = imp.getWidth() - 1
      ymax = imp.getHeight() - 1
      zmax = imp.getNSlices() - 1  
    
    save_shifts(shifts, [xmin, ymin, zmin, xmax, ymax, zmax])
    IJ.log("  saving shifts...")
    
run()

It seems to work well if not doing sub pixel drift correction.
If I toggle the box for sub pixel drift correction, I still get the xyz shifts reported in the log window (this was on purpose). However, the registered stack has extra pixels to the bottom right (example below). I am assuming that I failed to stop the script from creating a bigger image that takes into account the shifts in xy. Since I generally don’ t use the sub pixel option, this is not currently a real problem for me.

I might have broken more things than what I have noticed so far in this process. It is also not the most efficient way since lots of shifts that are not going to be used are still being calculated. But it might get the job done for now.


#14

Did you compile the plugin again after this? I have no idea how to put it back together? I’m keen to use your edit as then I might be able to remove all user input to my processing, a great success!


#15

@7rebor
I never compiled it.
I tested by just pasting it to File>New>Script… then set the language to python and hit "Run"
When I was satisfied with it, I followed the instructions at the address below to have it in my Plugin menu

Hope it works for you!