function Oblique_Reg(brikdir, func_brik, anat_brik,outname,target_Rx, ...
    input_Rx, target_angle_dir)
% This program registers axial anatomicals to oblique functionals (not vice versa)
% 
% Usage: Oblique_Reg(brikdir, func_brik, anat_brik,outname,[target_Rx],[input_Rx],[target_angle_dir])
% 
% Requried Inputs: 
%   brikdir: brik directory
%   func_brik: the prefix of the functional brik or full name if NIFTI format
%   anat_brik: the prefix of the anatomical brik or full name if NIFTI format 
%   outname:  name of the registerd anatomical dataset 
%              (NIFTI format if contains .nii, else AFNI brik format)
%
% Optional Inputs: (if the Rx information is not in the header, they become required inputs!!!)
%   target_Rx(optional): the Rx corrdinates of the target brik (using this format [-R +L -A +P -I +S])
%   input_rx(optional): the Rx coordinate of the input brik (same format)
%   target_angle_dir(obsolete): 1 or -1, fudge factor when program doesnt work. (1 or -1)
%                               


%**********AFNI BRIK ORIENTATION CONVENTION****************%
%    0 1 2 3 4 5 
%    R L P A I S
%    - + + - - +
%**********************************************************%
%    but Matlab matrix index starts with 1, so:
%    1 2 3 4 5 6 
%    R L P A I S
%    - + + - - +
%**********************************************************%
global out_brik input_brik target_brik outputnifti
outputnifti=0;

slab(1:6)=0;
dir_n=[-1 1 1 -1 -1 1];

func_brik=[brikdir '/' func_brik];
anat_brik=[brikdir '/' anat_brik];
outname  =[brikdir '/' outname];
kk=strfind(func_brik,'.nii');
if kk
    disp('Converting NIFTI to AFNI briks!')
    target_brik=[func_brik(1:(kk(end)-1)) '_tempxyz'];
    ev=sprintf('!3dcopy %s %s',func_brik,target_brik);
    eval(ev);
else
    target_brik=[func_brik '_tempxyz'];

    ev=sprintf('!3dcopy %s+orig %s',func_brik,target_brik);
    eval(ev);
end

kk=strfind(anat_brik,'.nii');
if kk
    disp('Converting NIFTI to AFNI briks!')
    input_brik=[anat_brik(1:(kk(end)-1)) '_tempxyz'];
    ev=sprintf('!3dcopy %s %s',anat_brik,input_brik);
    eval(ev);
else
    input_brik=[anat_brik '_tempxyz'];
    ev=sprintf('!3dcopy %s+orig %s',anat_brik,input_brik);
    eval(ev);
end

kk=strfind(outname,'.nii');
if kk
    outputnifti=1;
    out_brik=[outname(1:kk(end)-1) '_tempxyz'];
else
    outputnifti=0;
    out_brik=outname;
end
    
stat=sprintf('!3dAttribute ORIGIN %s+orig > temp.txt',target_brik);
eval(stat);
origin=load('temp.txt');

stat=sprintf('!3dAttribute ORIENT_SPECIFIC %s+orig > temp.txt',target_brik);
eval(stat);
orient_specific=load('temp.txt');

stat=sprintf('!3dAttribute DELTA %s+orig > temp.txt',target_brik);
eval(stat);
delta=load('temp.txt');

stat=sprintf('!3dAttribute DATASET_DIMENSIONS %s+orig > temp.txt',target_brik);
eval(stat);
dimensions=load('temp.txt');

stat=sprintf('!3dAttribute DATASET_DIMENSIONS %s+orig > temp.txt',input_brik);
eval(stat);
input_dimensions=load('temp.txt');
matrix_size=input_dimensions(1:2);


stat=sprintf('!3dAttribute DELTA %s+orig > temp.txt',input_brik);
eval(stat);
input_delta=load('temp.txt');

stat=sprintf('!3dAttribute ORIENT_SPECIFIC %s+orig > temp.txt',input_brik);
eval(stat);
orient=load('temp.txt');


if ~exist('target_Rx') | isempty(target_Rx)
   target_Rx=getRx(func_brik); 
end
if ~exist('input_Rx') | isempty(input_Rx)
   input_Rx=getRx(anat_brik); 
end

if ~exist('target_angle_dir') | isempty(target_angle_dir)
    if orient_specific(3)>orient_specific(2)
       target_angle_dir=-1;
    else
       target_angle_dir=1;
    end
end
!rm -f temp.txt
        
ind=orient_specific(1);
target_center(ind+1)=origin(1)+(dimensions(1)-1)*delta(1)/2;
target_center(ind+1+(-1)^ind)=origin(1)+(dimensions(1)-1)*delta(1)/2;

ind=orient_specific(2);
target_center(ind+1)=origin(2)+(dimensions(2)-1)*delta(2)/2;
target_center(ind+1+(-1)^ind)=origin(2)+(dimensions(2)-1)*delta(2)/2;

ind=orient_specific(3);
target_center(ind+1)=origin(3)+(dimensions(3)-1)*delta(3)/2;
target_center(ind+1+(-1)^ind)=origin(3)+(dimensions(3)-1)*delta(3)/2;

ind1=orient_specific(2);
ind2=ind1+(-1)^ind1;

if ind1>ind2
    startind1=ind2+1;
    endind1=ind1+1;
else
    startind1=ind1+1;
    endind1=ind2+1;
end

ind1=orient_specific(3);
ind2=ind1+(-1)^ind1;

if ind1>ind2
    startind2=ind2+1;
    endind2=ind1+1;
else
    startind2=ind1+1;
    endind2=ind2+1;
end

target_theta= target_angle_dir*atan((target_Rx(startind1)-...
    target_Rx(endind1))/(target_Rx(startind2)-target_Rx(endind2)));

input_theta=0;

shift(1)=(mean([input_Rx(1) input_Rx(2)])-mean([target_Rx(1) target_Rx(2)]));
shift(2)=shift(1);
    
shift(3)=(mean([input_Rx(3) input_Rx(4)])-mean([target_Rx(3) target_Rx(4)]))*cos(target_theta)+...
    (mean([input_Rx(5) input_Rx(6)])-mean([target_Rx(5) target_Rx(6)]))*sin(target_theta);
shift(4)=shift(3);

shift(5)=(mean([input_Rx(5) input_Rx(6)])-mean([target_Rx(5) target_Rx(6)]))*cos(target_theta)-...
    (mean([input_Rx(3) input_Rx(4)])-mean([target_Rx(3) target_Rx(4)]))*sin(target_theta);
shift(6)=shift(5);

    
ind=orient(1);
input_dim(ind+1)=-(input_dimensions(1)-1)*input_delta(1)/2;
input_dim(ind+1+(-1)^ind)=-input_dim(ind+1);

ind=orient(2);
input_dim(ind+1)=-(input_dimensions(2)-1)*input_delta(2)/2;
input_dim(ind+1+(-1)^ind)=-input_dim(ind+1);

ind=orient(3);
input_dim(ind+1)=-(input_dimensions(3)-1)*input_delta(3)/2;
input_dim(ind+1+(-1)^ind)=-input_dim(ind+1);


for i=1:6
    slab(i)=target_center(i)+shift(i)+input_dim(i);
end

str_x=sprintf('-xorigin_raw %s',num2str(slab(orient(1)+1)));
str_y=sprintf('-yorigin_raw %s',num2str(slab(orient(2)+1)));
str_z=sprintf('-zorigin_raw %s',num2str(slab(orient(3)+1)));

eval(sprintf('!rm -f %s_temp+orig*',input_brik));
stat=sprintf('!3dcopy %s+orig %s_temp', input_brik, input_brik);
fprintf('executed:%s\n',stat)
eval(stat);
 
stat=sprintf('!3drefit %s %s %s %s_temp+orig', str_x, str_y, str_z,input_brik);
fprintf('executed:%s\n',stat)
eval(stat);

rot_angle=(input_theta-target_theta)*180/pi;
if (rot_angle~=0)
    stat=sprintf('!3drotate -verbose -prefix %s -rotate %sR 0A 0S %s_temp+orig', ...
        out_brik,num2str(rot_angle), input_brik);
    fprintf('executed:%s\n',stat)
    eval(stat);
else
    stat=sprintf('!3dcopy %s_temp+orig %s', input_brik, out_brik);
    fprintf('executed:%s\n',stat)
    eval(stat);
end

if outputnifti
     stat=sprintf('!3dAFNItoNIFTI -prefix %s %s+orig',outname, out_brik);
     fprintf('executed:%s\n',stat)
     eval(stat);
end

obclean();
disp('Done!');


%///////////////////////    obclean ()  //////////////////////////////
function obclean()
global out_brik input_brik target_brik outputnifti
if exist([out_brik '+orig.BRIK'])==2 | exist([out_brik '+orig.HEAD'])==2
   if outputnifti
       eval(['!rm -f ' out_brik '+orig.*']);
   end
end
if exist([input_brik '_temp+orig.BRIK'])==2 | exist([input_brik '_temp+orig.HEAD'])==2
   eval(sprintf('!rm -f %s_temp+orig*',input_brik));
end

if exist([input_brik '+orig.BRIK'])==2 | exist([input_brik '+orig.HEAD'])==2
   eval(sprintf('!rm -f %s+orig*',input_brik));
end

if exist([target_brik '+orig.BRIK'])==2 | exist([target_brik '+orig.HEAD'])==2
   eval(sprintf('!rm -f %s+orig*',target_brik));
end

if exist('temp.txt')==2 
   eval('!rm -f temp.txt');
end
if exist('AFNINotes')==2 
   eval('!rm -f AFNINotes');
end

%///////////////////////  getRx  //////////////////////////////
function Rxx=getRx(brikname)
global out_brik input_brik target_brik outputnifti

RxStr='';
sk=strfind(brikname,'.nii');

if (sk) %NIFTI format
    command=['nifti_tool -disp_exts -infile ', brikname, ' > AFNINotes' ];
    system(command);
else
    command = ['3dNotes ',brikname,'+orig > AFNINotes'];
    system(command);
end
RxStart='';
RxEnd='';
fid = fopen('AFNINotes');
readline=1;
while readline
    ss = fgetl(fid);
    if ss==-1 readline=0;
    elseif (~isempty(ss) & length(ss) >=8)
        if strfind(ss,'Rx info:')
        %if strcmp(ss(1:8),'Rx info:');
            ss = fgetl(fid);
            [RxStart] = strread(ss, '%s');        % Rx info: start values here,
            ss = fgetl(fid);                                % end values on following line
            [RxEnd] = strread(ss, '%s');
            readline=0;
        end
    end
end
fclose(fid);
if isempty(RxStart) | isempty(RxEnd)
    obclean();
    error(['No Rx info found in ' brikname ' header! Please provide Rx as inputs.']) ;
end

Rx={RxStart{2},RxEnd{2},RxStart{3},RxEnd{3},RxStart{4},RxEnd{4}};

eval('!rm -f AFNINotes');

tmpRx=zeros(1,6);
for cntRx=1:6
    if strfind(Rx{cntRx},'L')
        tmpRx(cntRx) = str2num(Rx{cntRx}(strfind(Rx{cntRx},'L')+1:end));
    elseif strfind(Rx{cntRx},'P')
        tmpRx(cntRx) = str2num(Rx{cntRx}(strfind(Rx{cntRx},'P')+1:end));
    elseif strfind(Rx{cntRx},'S')
        tmpRx(cntRx) = str2num(Rx{cntRx}(strfind(Rx{cntRx},'S')+1:end));
    elseif strfind(Rx{cntRx},'R')
        tmpRx(cntRx) = -str2num(Rx{cntRx}(strfind(Rx{cntRx},'R')+1:end));
    elseif strfind(Rx{cntRx},'A')
        tmpRx(cntRx) = -str2num(Rx{cntRx}(strfind(Rx{cntRx},'A')+1:end));
    elseif strfind(Rx{cntRx},'I')
        tmpRx(cntRx) = -str2num(Rx{cntRx}(strfind(Rx{cntRx},'I')+1:end));
    end
end
Rxx=tmpRx;