alignment_test.m

## Copyright (C) 2007 P. Cloetens
## 
## 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 2 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, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

## alignment
## function [offset(,tilt)] = alignment(varargin)
##      Tomography alignment
##      determines offset (and tilt) of rotation axis
##      uses two images acquired at 180 degrees apart, a reference image and a dark image
##      assumptions:
##      vertical rotation axis
##      offset motor translates from left to right in the image
##      images are manually given or stored in /data/beamline/inhouse/align
##      with beamline either id19 or set by the global variable FT_BL
##      align0000.edf: image at 0 degrees
##      align0001.edf: image at 180 degrees
##      align0002.edf: reference image
##      align0002.edf: dark image
##      tilt (degrees) is determined only when axisposition method is (choose)highlow
##      offset (pixels) gives current axis position with respect to image middle
##          corrections to be applied are given when possible
##
##      arguments:
##      argument 1: axisposition method ( default: 'highlow' )
##          possible values:
##              highlow: line by line correlation, gives offset and tilt
##              global: 2D correlation, gives offset only
##              accurate: 2D correlation, followed by accurate real space determination
##              choosehighlow: same as highlow, but top and bottom are chosen graphically
##
##      e.g.:   [offset,tilt] = alignment;
##              [offset,tilt] = alignment('highlow');
##              offset = alignment('global');
##              offset = alignment('accurate');
##              [offset,tilt] = alignment('choosehighlow');

## Author: P. Cloetens <cloetens@esrf.fr>
## 
## 2007-07-13 P. Cloetens <cloetens@esrf.fr>
## * Initial revision
## * Based on alignment.m by Wolfgag and Peter

function [offset,tilt] = alignment(axisposition);
    global FT_BL

    if isempty(FT_BL)
        beamline = 'id19';
    else
        beamline = FT_BL;
    end
    if !nargin
        axisposition = [];
    endif
    if isempty(axisposition)
        axisposition = "highlow";    
    endif

    headerlength=1024;
    # defaults if semi-automatic alignment
    align_dir = ['/data/' beamline '/inhouse/align'];
    align_prefix = 'align';
    nbeg = 0; # image at 0 degrees
    nend = 1; # image at 180 degrees
    nref = 2; # reference image
    ndark = 3; # dark image

    #### message to check beamline is right
    printf('You are working on beamline %s\n',beamline)

    ####### semi-automatic acquisition ? ########
    name=sprintf('%s/%s%4.4i.edf',align_dir,align_prefix,nend);
    descript = dir(name);
    if isempty(descript)
        align_auto = 0;
    elseif ((datenum(clock)-datenum(descript.date))*24 > 0.5)
        # check if images for alignment are more than 30 minutes old
        align_auto = 0;
    else
        align_auto = 1;
    end	

    ####### initialisation ########
    if align_auto
        n1 = align_dir;
        n2 = align_prefix;
    else
        # manual initialisation 
        # n1: directory, default is present directory
        disp('You will have to enter values manually')
        disp('Please use alignment in *tomo for semi-automatic alignment')
        # n1: directory, default is present directory
        n1=cleandirectoryname;
        # prefix
        n2 = input('prefix ?        : ','s');
        ndark = input('number of dark image ?       : ');
        nref = input('number of reference image ?  : ');
        nbeg = input('number image 0 degrees?      : ');
        nend = input('number image 180 degrees?    : ');
    end # initialisation

    ##########
    # reading images
    ##########

    name=sprintf('%s/%s%4.4i.edf',n1,n2,ndark);
    d=edfread(name);
    name=sprintf('%s/%s%4.4i.edf',n1,n2,nref);
    f=edfread(name);
    name=sprintf('%s/%s%4.4i.edf',n1,n2,nbeg);
    a0=edfread(name);
    name=sprintf('%s/%s%4.4i.edf',n1,n2,nend);
    a1=edfread(name);
    fp=fopen(name,'rb');
    hd=fscanf(fp,'%c',headerlength);
    fclose(fp);
    pixsize=findheader(hd,'pixel_size','float');
    if isempty(pixsize)
        pixsize=findheader(hd,'optic_used','float');
        if isempty(pixsize)
            disp('pixel_size and optic_used are not in header')
            disp('Pixel size is unknown')
        endif
    endif

    ##########
    # selecting region, flatfield and mirror of image 180 degrees
    ##########


    if strcmp(axisposition,'choosehighlow')
        [rb,re,cb,ce]=imagej_getcorners(a0,'Indicate top of region for alignment','Indicate bottom of region for alignment');
    else
        cb = 1;
        ce = size(a0,2);
    endif
    a0=(a0(:,cb:ce)-d(:,cb:ce))./max((f(:,cb:ce)-d(:,cb:ce)),1);
    a1=(a1(:,cb:ce)-d(:,cb:ce))./max((f(:,cb:ce)-d(:,cb:ce)),1);
    a1=flipud(a1);

    ##########
    # determination tilt and offset
    ##########

    [n,m]=size(a0);
    if strcmp(FT_BL,'id22')
        offset_motor_name = 'ys';
    else
        offset_motor_name = 'pmy';
    endif
    limz=0.5; # maximum vertical drift before warning 
    switch axisposition
    case {'global','accurate'}
        sizepow2=2^floor(log(m)/log(2)); # to use a good size for fft's
        shift=correlate(cut(a0,n,sizepow2),cut(a1,n,sizepow2),0);
        if strcmp(axisposition,'accurate')
            if (abs(shift(2))>10*limz)
                # this is suspicious, we don't trust result of correlate
                printf('Pre-correlation yields %5.2f pixels vertical motion\n',shift(2))
                printf('We do not consider this\n')
                shift=[0 0];
            endif                
            # limit the size of region for comparison to cutsize in both directions
            cutsize = 1024;
            oldshift = round(shift);
            if ((n>cutsize)|(m>cutsize))
                shift = oldshift + findshift(cut(a0,min(n,cutsize),min(m,cutsize)),cut(circshift(a1,oldshift),min(n,cutsize),min(m,cutsize)));
            else
	        shift = findshift(a0,a1,[],oldshift);
	    endif
            if (norm(shift-oldshift)>2)
                printf('Pre-correlation and accurate correlation are not consistent\n')
                printf('!!!Please check or use manual!!!\n')
            endif
        endif
        offset=shift(1)/2;
        if abs(shift(2))>limz
	    disp('Verify alignment or motion sample !!!')
	    printf('Vertical motion: %5.2f pixels\n',shift(2))
        end
    case {'highlow','choosehighlow'}
        c=correlate1(a0,a1);
        plot(c);
        y=1:m;
        dec=polyfit(y,c,1);
        tilt=-dec(1)/2*180/pi;
        offset=(dec(2)+m/2*dec(1))/2;
        printf('Tilt to apply (degrees): mvr rotc %5.3f\n',tilt)
    endswitch
    if !isempty(pixsize)
        printf('Pixel size is %3.1f microns\n',pixsize)
        printf('Offset to apply (mm) = %5.2f pixels: mvr %s %6.4f\n',-offset,offset_motor_name,-offset*pixsize/1000)
    else
        printf('Offset to apply = %5.2f pixels\n',-offset)
    endif
endfunction