function Y = convertPW_MD(X,sn)

% Usage: Y = convertPW_MD(X,sn)
% 
% This function converts coordinates from the Paxinos & Watson (PW) Space
% to the Minimum Deformation (MD) Space and the other way around. It uses
% the SPM8 transformations 'MD2PW.mat' and 'PW2MD.mat'. These were obtained
% using SPM8 "Unified Segmentation" (Ashburner, J., Friston, K.J., 2005.
% Neuroimage 26, 839-51). 
% _________________________________________________________________________
% Case 1: Converting coordinates from MD Space to PW Space
% --------------------------------------------------------
% Y_PW = convertPW_MD(X_MD,'MD2PW.mat');
% 
% X_MD: Coordinates (Nx3 array) in the MD Space
% Y_MD: Converted Coordinates (Nx3 array) in the PW Space
% _________________________________________________________________________
% Case 2: Converting coordinates from PW Space to MD Space
% --------------------------------------------------------
% Y_PW = convertPW_MD(X_MD,'PW2MD.mat');
% 
% X_PW: Coordinates (Nx3 array) in the PW Space
% Y_MD: Converted coordinates (Nx3 array) in the MD Space
%
% Pedro Antonio Valdes-Hernandez

prm = load(sn);
prm.VG(1).dim(1) = prm.VG(1).dim(1);
prm.VG(1).dim(2) = prm.VG(1).dim(2);
prm.VG(1).dim(3) = prm.VG(1).dim(3);
xg = 1:prm.VG(1).dim(1);
yg = 1:prm.VG(1).dim(2);
zg = 1:prm.VG(1).dim(3);
prm.Affine = prm.Affine/prm.VG(1).mat;
prm.VG(1).mat(1:3,4) = prm.VG(1).mat(1:3,4);
prm.Affine = prm.Affine*prm.VG(1).mat;
[Xg,Yg] = ndgrid(xg,yg);
X = applyM(X,inv(prm.VG(1).mat));
xng = X(:,1);
yng = X(:,2);
zng = X(:,3);
BX = spm_dctmtx(prm.VG(1).dim(1),size(prm.Tr,1),xg-1);
BY = spm_dctmtx(prm.VG(1).dim(2),size(prm.Tr,2),yg-1);
BZ = spm_dctmtx(prm.VG(1).dim(3),size(prm.Tr,3),zg-1);
hx = zeros(prm.VG(1).dim(1),prm.VG(1).dim(2),prm.VG(1).dim(3));
hy = zeros(prm.VG(1).dim(1),prm.VG(1).dim(2),prm.VG(1).dim(3));
hz = zeros(prm.VG(1).dim(1),prm.VG(1).dim(2),prm.VG(1).dim(3));
if not(isempty(prm.Tr))
    for j = zg
        tx = get_2Dtrans(prm.Tr(:,:,:,1),BZ,j);
        ty = get_2Dtrans(prm.Tr(:,:,:,2),BZ,j);
        tz = get_2Dtrans(prm.Tr(:,:,:,3),BZ,j);
        hx(:,:,j) = Xg    + BX*tx*BY';
        hy(:,:,j) = Yg    + BX*ty*BY';
        hz(:,:,j) = zg(j) + BX*tz*BY';
    end
    c = spm_bsplinc(hx,[1 1 1 0 0 0]);
    x = spm_bsplins(c,double(xng),double(yng),double(zng),[1 1 1 0 0 0]);
    c = spm_bsplinc(hy,[1 1 1 0 0 0]);
    y = spm_bsplins(c,double(xng),double(yng),double(zng),[1 1 1 0 0 0]);
    c = spm_bsplinc(hz,[1 1 1 0 0 0]);
    z = spm_bsplins(c,double(xng),double(yng),double(zng),[1 1 1 0 0 0]);
else
    x = xng; y = yng; z = zng;
end
Y = applyM([x y z],prm.VF.mat*prm.Affine);
%---

function x = applyM(x,M)

S = size(M);
if all(S == [3 3])
    x = x*M';
elseif all(S == [4 4])
    x = [x ones(size(x,1),1)]*M';
    x(:,4) = [];
else
    error('incorrect dimensions of the transformation matrix');
end

function T2 = get_2Dtrans(T3,B,j)
d   = [size(T3) 1 1 1];
tmp = reshape(T3,d(1)*d(2),d(3));
T2  = reshape(tmp*B(j,:)',d(1),d(2));