function [] = write_Mat2mhd(Matrice, Path_sortie, Pixel_spacing)

%% Exports a 1D, 2D, 3D, 4D Matlab object
%% as a VTK/ITK readable '.mhd' file.
%% (actually : n files for 4D objects)

%% Matrice       : 1D, 2D, 3D, 4D
%% Path_sortie   : Output file name, *no* extention ( ex : myFile, *not* myFile.mhd )
%% Pixel_spacing : pixel sizes, in real world ( ex : [1.5 1.5 4] )

    nbPixelSpacing = size(Pixel_spacing);

    ndim = numel(size(Matrice));

%% WTF ? JPR
    Path_sortie = char(Path_sortie);

%% /// \TODO : throw an exception 'dim > 4 not yet dealt with' ?
 if ndim > 4
    err = MException('ResultChk:OutOfRange','Object dim> 4');
    errCause = MException('ResultChk:BadInput','Object dim > 4');
    err = addCause(err, errCause);
    throw(err);
 end

%% /// \TODO
%% throw exception if ndim vs nbPixelSpacing are unconsistent?
%% or apply an heuristics : a single value means pixels are isotropics?
%%                          no value means pixel Size is [1.0 1.0 1.0]?

Pixel_spacing = num2str(Pixel_spacing);

%% /// \TODO : Check if we *actually* get a vector of strings

%% ---
%% --- Check Pixel Type ---
%% ---

    tab_matlab_types = {'int8',    'uint8',    'int16',    'uint16',    'int32',  'uint32',  'int64',   'uint64',   'single',   'double'    };
    tab_mhd_types    = {'MET_CHAR','MET_UCHAR','MET_SHORT','MET_USHORT','MET_INT','MET_UINT','MET_LONG','MET_ULONG','MET_FLOAT','MET_DOUBLE'};
    tab_pixel_sizes  = { 1,         1,          2,          2,           4,        4,         8,         8,          4,          8          };
  
    [a,b] = size(tab_matlab_types);
 %{
    for i=1:b
        K = isa(Matrice, tab_matlab_types{i});
        if K==1
            mhd_type   = tab_mhd_types{i};
            pixel_size = tab_pixel_sizes{i};
            break;
        end
    end
 %}

    k = class(Matrice);

    for i = 1:b
        test(i)= strcmp(k,tab_matlab_types(i));
    end

    choice = find(test==1);

    mhd_type   = tab_mhd_types{choice};
    pixel_size = tab_pixel_sizes{choice};

%% ---
%% --- Transpose the matrix ---
%% --- 
 
    if ndim == 1
        %% nothing to do ?
%%	/// \TODO : Voir pb vecteur ligne / vecteur colonne
        Matrice_a = Matrice;
        %%[Largeur] = size(Matrice_a);
        %%C2 = num2str([Largeur]);
       
    elseif ndim == 2
        %% matlab stores 2D objects column per column, C/C++ stores row per row
        Matrice_a = transpose(Matrice);
        %%[Hauteur,Largeur] = size(Matrice_a);
        %%C2 = num2str([Hauteur,Largeur]);

    elseif ndim == 3
        [haut,largeur,prof] = size(Matrice);
        %% /// \TODO : check the *actual* order of [haut,largeur,prof) !      
        for i=1:prof
            Matrice_a(:,:,i) = transpose(Matrice(:,:,i));
        end
        
        %%[Hauteur,Largeur,Profondeur] = size(Matrice_a);
        %%C2 = num2str([Hauteur,Largeur,Profondeur]);

    elseif ndim == 4
       %% export - n .mhd files, each one describing a 3D object (Headersize =)
       %%        - a single .raw file
       %%        
       %% supose 4th D is time
       %% [Hauteur,Largeur,Profondeur,Instants] = size(Matrice_a);
       %% C2 = num2str([Hauteur,Largeur,Profondeur,Instants]);
       %% 
        [haut,largeur,prof,instants] = size(Matrice);
       %% /// \TODO : check the *actual* Matlab order of [haut,largeur,prof,instants] !     
        for i=1:instants
           for j=1:prof
               Matrice_a(:,:,i,j) = transpose(Matrice(:,:,i,j));
           end
        end

    else
       %% Here, deal with dim > 4?
       %% (NOT AT ALL a hurry!)
       %%
       %% -> error 'Object dim> 4' already trapped.
    end
    
%% ---
%% --- Write the Pixels in a single '.raw' file---
%% ---     
    fid = fopen([Path_sortie,'.raw'],'w+');    
    fwrite(fid, Matrice_a, tab_matlab_types{choice});

%% ---
%% --- Trick for 4D objects ---
%% ---    
    if ndim < 4
       instants = 1;
       imagesize = numel(Matrice_a)*pixel_size;     
    else
       imagesize =  haut*largeur*prof*pixel_size;    
    end

%% ---
% --- Write the .mhd file(s) ---
% ---
% a .mhd file looks like :
%====================================
% ObjectType = Image
% NDims = 3
% DimSize = 127 147 187
% ElementType = MET_USHORT
% HeaderSize = 0
% ElementNumberOfChannels = 1
% ElementSpacing = 0.7 0.7 0.7
% Position = 0 0 0
% TransformMatrix = 1 0 0 0 1 0 0 0 1
% Offset = 0 0 0
% CenterOfRotation = 0 0 0
% BinaryData = True
% BinaryDataByteOrderMSB = False
% CompressedData = False
% ElementDataFile = hola.raw
%% ====================================

    %% // \TODO
    %% Check if it works   
    C2 = num2str(size(Matrice_a));

    for ins=1:instants
       A = ['ObjectType = Image']; 
       B = ['NDims = ', num2str(ndim)];  
       C = ['DimSize = ', C2];  
       D = ['ElementType =', mhd_type];
       E = ['HeaderSize = ', num2str((ins-1)*imagesize)] ;
       F = ['ElementNumberOfChannels = 1'] ;
       G = ['ElementSpacing = ', num2str(Pixel_spacing)];
       H = ['BinaryData = True'];
       I = ['ElementByteOrderMSB = False'];
       J = ['CompressedData = False'];
       K = ['ElementDataFile = ', Path_sortie,'.raw'];

       Header = strvcat(A,B,C,D,E,F,G,H,I,J,K);
       if instants < 4
          dlmwrite([Path_sortie, '.mhd'], Header, 'delimiter', '', 'newline', 'pc') ;
        else
          dlmwrite([Path_sortie, num2str(ins), '.mhd'], Header, 'delimiter', '','newline', 'pc') ;
        end
    end