clear all;
addpath(genpath('/home/dch/ismrmrd/matlab'));
%% Generating a agilent ISMRMRD data set

% File Name

importname = 'radialFSE_13C_golden_02';
fidname = strcat(importname,'.fid');
filename = strcat(importname,'.h5');
% Create an empty ismrmrd dataset
if exist(filename,'file')
    delete(filename); % Hmmm, doesn't work
    %error(['File ' filename ' already exists.  Please remove first'])
end
dset = ismrmrd.Dataset(filename);

imagesize = 64;
%% Read in Agilent data
[RE,IM,NP,NB,NT,HDR]=load_fid(importname);


nCoils = length(queryPP(fidname,'rcvrs'))-2

nX = queryPP(fidname,'np')/2;
nY =queryPP(fidname,'nv');
nZ = 1;

pad = length(queryPP(fidname,'pad'));

tr=queryPP(fidname,'tr');
te=queryPP(fidname,'te');
radialAngleRange = queryPP(fidname,'radialAngleRange')
pelist=queryPP(fidname,'pelist');
%pe2list=queryPP(fidname,'pe2list');

fid = complex(RE,IM);
fid = reshape(fid,[nX nY nCoils pad]);


%% Compute trajectory
matrix =linspace(-nX/2,nX/2-3,nX);
matrix = matrix/nX; %scaling required by Gadgetron
matrix = matrix;
%phi=linspace(0,radialAngleRange,nY+1);
%phi = phi(1:end-1);
phi = pelist*111.25;


clear Kx;
clear Ky;
%clear Kz;

% Our reconstruction consists of a series of 2D iNFFTs
% - with corresponding 2d trajectories (identical for each slice).
for Y=1:nY
    Kx(:,Y) = matrix*cos(phi(Y)*pi/180);
    Ky(:,Y) = matrix*sin(phi(Y)*pi/180);
end
k = reshape(complex(Ky,Kx),nCoils,nY*nX);


acqblock = ismrmrd.Acquisition(nY);

% Set the header elements that don't change
acqblock.head.version(:) = 1.0;
acqblock.head.number_of_samples(:) = nX;
acqblock.head.center_sample(:) = floor(nX/2);
acqblock.head.active_channels(:) = nCoils;

% Set the number of trajectory dimension 
acqblock.head.trajectory_dimensions(:) = 3;

% !!! TODO !!! 
% Set the matrix orientation for the current slice.
% We really should set these direction vectors to
% provide the actual rotation matrix of the current slice.
acqblock.head.read_dir  = repmat(single([1 0 0]'),[1 nY]);
acqblock.head.phase_dir = repmat(single([0 1 0]'),[1 nY]);
acqblock.head.slice_dir = repmat(single([0 0 1]'),[1 nY]);

dcw = zeros(1,nX);
%fid = ifftshift(fft(abs(ifft(fid,[],1)),[],1),1);
% Loop over the acquisitions, set the header, set the data and append
for rep = 1:1
    for acqno = 1:nY
        
        % Set the header elements that change from acquisition to the next
        % c-style counting
        acqblock.head.scan_counter(acqno) = (rep-1)*nY + acqno-1;
        acqblock.head.idx.kspace_encode_step_1(acqno) = acqno-1; 
        acqblock.head.idx.repetition(acqno) = rep - 1;
        
        % Set the flags
        acqblock.head.flagClearAll(acqno);
        if acqno == 1
            acqblock.head.flagSet(acqblock.head.FLAGS.ACQ_FIRST_IN_ENCODE_STEP1, acqno);
            acqblock.head.flagSet(acqblock.head.FLAGS.ACQ_FIRST_IN_SLICE, acqno);
            acqblock.head.flagSet(acqblock.head.FLAGS.ACQ_FIRST_IN_REPETITION, acqno);
        elseif acqno==nY
            acqblock.head.flagSet(acqblock.head.FLAGS.ACQ_LAST_IN_ENCODE_STEP1, acqno);
            acqblock.head.flagSet(acqblock.head.FLAGS.ACQ_LAST_IN_SLICE, acqno);
            acqblock.head.flagSet(acqblock.head.FLAGS.ACQ_LAST_IN_REPETITION, acqno);
        end
        
        % fill the data
        
        acqblock.data{acqno} = squeeze(fid(:,acqno,rep));

        % set the trajectories
        K = zeros(3*nX,1);
        K(1:3:end) = Kx(:,acqno); %2D
        K(2:3:end) = Ky(:,acqno); %2D
        %K(3:3:end) = dcf(1,:,acqno);
        K(3:3:end) = 2*abs(matrix);
        
        

        acqblock.traj{acqno} = squeeze(K(:));   
      
     
    end

    % Append the acquisition block
    dset.appendAcquisition(acqblock);
        
end % rep loop


%%%%%%%%%%%%%%%%%%%%%%%%
%% Fill the xml header %
%%%%%%%%%%%%%%%%%%%%%%%%

%% Experimental Conditions
header = struct();
header.experimentalConditions.H1resonanceFrequency_Hz = 400000000; % 9.4T
% 
% expcond = org.ismrm.ismrmrd.ExperimentalConditionsType;
% %expcond.setH1ResonanceFrequencyHz(128000000); % 3T
% expcond.setH1ResonanceFrequencyHz(400000000); % 9.4T
% cd 
% dset.xmlhdr.setExperimentalConditions(expcond);

%% Acquisition System Information
% acqsysinfo = org.ismrm.ismrmrd.AcquisitionSystemInformationType;
% acqsysinfo.setReceiverChannels(java.lang.Integer(nCoils));
header.acquisitionSystemInformation.receiverChannels = nCoils;
% dset.xmlhdr.setAcquisitionSystemInformation(acqsysinfo);

%% The Encoding
% encoding = org.ismrm.ismrmrd.Encoding;
% %encoding.setTrajectory(org.ismrm.ismrmrd.TrajectoryType.CARTESIAN);
% encoding.setTrajectory(org.ismrm.ismrmrd.TrajectoryType.RADIAL); %?
encoding = struct();
encoding.trajectory = 'radial';



% Encoded Space
% fov = org.ismrm.ismrmrd.FieldOfViewMm;
% fov.setX(queryPP(fidname,'lro'));
% fov.setY(queryPP(fidname,'lpe'));
% fov.setZ(queryPP(fidname,'lpe2'));

encodedSpace = struct();
encodedSpace.fieldOfView_mm = struct('x',queryPP(fidname,'lro'), 'y', queryPP(fidname,'lpe'), 'z', queryPP(fidname,'lpe2') );

encodedSpace.matrixSize = struct('x',imagesize, 'y',imagesize, 'z', 1);
% matrix = org.ismrm.ismrmrd.MatrixSize;
% matrix.setX(nX);
% matrix.setY(nY);
% matrix.setZ(nZ);

encoding.encodedSpace = encodedSpace;
% encodedspace = org.ismrm.ismrmrd.EncodingSpaceType;
% encodedspace.setMatrixSize(matrix);
% encodedspace.setFieldOfViewMm(fov);
% encoding.setEncodedSpace(encodedspace);

% Recon Space
% (same as encoding space)
encoding.reconSpace =encodedSpace;
% encoding.setReconSpace(encodedspace);

% Encoding Limits

encodingLimits = struct();
kspace_encoding_step_0 = struct();
kspace_encoding_step_0.minimum = 0;
kspace_encoding_step_0.center = floor(nX/2);
kspace_encoding_step_0.maximum = nX-1;
encodingLimits.kspace_encoding_step_0 = kspace_encoding_step_0;
% encodinglimits = org.ismrm.ismrmrd.EncodingLimitsType;
% 
% limitsa = org.ismrm.ismrmrd.LimitType;
% limitsa.setMinimum(0);
% limitsa.setCenter(floor(nX/2)); 
% limitsa.setMaximum(nX-1); 
% encodinglimits.setKspaceEncodingStep0(limitsa);

kspace_encoding_step_1 = struct();
kspace_encoding_step_1.minimum = 0;
kspace_encoding_step_1.center = floor(nY/2);
kspace_encoding_step_1.maximum = nY-1;
encodingLimits.kspace_encoding_step_1 = kspace_encoding_step_1;

% limitsb = org.ismrm.ismrmrd.LimitType;
% limitsb.setMinimum(0);
% limitsb.setCenter(floor(nY/2)); 
% limitsb.setMaximum(nY-1); 
% encodinglimits.setKspaceEncodingStep1(limitsb);




% limitsc = org.ismrm.ismrmrd.LimitType;
% limitsc.setMinimum(0);
% limitsc.setCenter(floor(nZ/2)); 
% limitsc.setMaximum(nZ-1);
% encodinglimits.setRepetition(limitsc);



% All the rest are zero and we can use the same one for all of them
limits = struct('minimum',0,'center',0,'maximum',0);

% limits = org.ismrm.ismrmrd.LimitType;
% limits.setMinimum(0);
% limits.setCenter(0); 
% limits.setMaximum(0);
encodingLimits.repetition = limits;
encodingLimits.average = limits;
encodingLimits.slice = limits;
encodingLimits.contrast = limits;
encodingLimits.phase = limits;

encodingLimits.set = limits;
encodingLimits.segment = limits;
encodingLimits.kspace_encoding_step_2 = limits;
encoding.encodingLimits = encodingLimits;

% encodinglimits.setAverage(limits);
% encodinglimits.setContrast(limits);
% encodinglimits.setKspaceEncodingStep2(limits);
% encodinglimits.setPhase(limits);
% encodinglimits.setSegment(limits);
% encodinglimits.setSet(limits);
% encodinglimits.setSlice(limits);
% Stuff
% encoding.setEncodingLimits(encodinglimits);
% dset.xmlhdr.getEncoding.add(encoding);

header.encoding = encoding;

%% Subject information
subjectInformation = struct();
% subjectInformation.patientName = 'Rat O Ratson';
% subjectInformation.patientWeight_kg = 0.1;
% subjectInformation.patientID = 42;
% subjectInformation.patientBirthdate = '20000101';
% subjectInformation.patientGender = 'o';
header.subjectInformation  = subjectInformation;
%% Study Information
studyInformation = struct();
studyInformation.studyDate = '20000101';
studyInformation.studyTime = '0000';
studyInformation.studyID = '42';
studyInformation.accessionNumber = '1';
studyInformation.referringPhysicianName = 'Dr. Who';
studyInformation.studyDescription = '';
studyInformation.studyInstanceUID = '1.2.3.4.5.6.7';
header.studyInformation = studyInformation;

%% MeasurementInformation
measurementInformation = struct();
measurementInformation.patientPosition='prone';

measurementDependency = struct();
measurementDependency.dependencyType='x';
measurementDependency.measurementID = '42';
measurementInformation.measurementDependency = measurementDependency;
referencedImageSequence = struct();
referencedImageSequence.referencedSOPInstanceUID = 'x';
measurementInformation.referencedImageSequence = referencedImageSequence;
header.measurementInformation = measurementInformation;

%% AcquisitionSystemInformation
acquisitionSystemInformation = struct();
header.acquisitionSystemInformation = struct();

%% Sequence parameters

sequenceParameters = struct();
sequenceParameters.TR = 1.0;
sequenceParameters.TE = 1.0;
sequenceParameters.TI = 1.0;
sequenceParameters.flipAngle_deg = 90.0;

header.sequenceParameters = sequenceParameters;

%% Write the dataset
% dset.close();

xmlheader = ismrmrd.xml.serialize(header);
dset.writexml(xmlheader);

D = dset.readAcquisition();
dset.close();
