%% bflyMotionEstimate.m script
%%  This is an example of estimating linear motion from 
%%  3D Butterfly navigator data. 
%%
%%  This example uses the function transEstAll3(). There
%%  is a completely matlab-based implementation of the 
%%  function and a compilable mex version (with omp/acml) 
%%  that runs significantly faster.
%%
%% (c) Joseph Y Cheng (jycheng@mrsrl.stanford.edu) 2012

%% SVN info: 
%%   Date:     $Date: 2014-10-17 17:10:46 -0700 (Fri, 17 Oct 2014) $
%%   Revision: $Revision: 1 $
%%   Author:   $Author: jycheng $
%%   Id:       $Id: runBflyMotionEstimate.m 1 2014-10-18 00:10:46Z jycheng $


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load Data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;
    dbdisp('Loading example data...');
    load example1.mat;
toc
% The following data structures are stored in the .mat file:
%   DATAA          -- raw k-space data with size nx*ny*nz*nc:
%                       nx -- readout length
%                       ny -- number of phase-encodes
%                       nz -- number of slice-encodes
%                       nc -- number of coils (receivers)
%   NAVA           -- raw navigator data with size nnav*nall:
%                       nnav -- navigator length
%                       nall -- total number of readouts x nc (= ny*nz*nc)
%   TR             -- [us] repitition time
%   k              -- [1/cm] k-space trajectory for navigator data
%   nFRead         -- reconstruction dimension for x 
%                     (for partial k-space imaging in x)
%   res            -- [mm] image resolution
%   yorder, zorder -- y-phase encode and z-phase encode ordering

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Determine dimensions from input data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nx,ny,nz,nc] = size(DATAA);
[nnav,nall]   = size(NAVA);
[npe]         = length(yorder); % same length as zorder

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters that can be adjusted
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NREF  = round(npe/2); % reference point for motion estimation
FILTW = floor(nz/4);  % filter length to filter motion estimation
PFIT  = 1;            % polynomial-fitting order to remove any
                      %   motion estimation drifiting
DEBUG_MOTION = 200;   % flag to plot motion in figure #200

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate Motion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%addpath([cd '/offreconmod']);
tic;
    dbdisp('Estimating motion...');
    [DDX,DDY,DDZ] = transEstAll(NAVA, k, res, ...
                                nc, npe, ...
                                NREF, PFIT, FILTW,...
                                DEBUG_MOTION);
toc