%% This function creates Pinfo struct for the selected .spar file path on the MAIN GUI Page.
% Pinfo includes main path, sparname, coordpath, spectrapath, ncol, nrow,
% nsli, voxoffap, voxofflr, voxoffcc, voxangap, voxanglr, voxangcc, apVOI,
% lrVOI, ccVOI, FOV, slidistance, slithickness, littlevoxels: [3×14×14
% struct] which includes each voxels' location in anatomical position.

% Inputs: filename and pathname of the .spar file
% Output: Pinfo struct

% Author: Sevim Cengiz, Bogazici University, Turkey (2020-11-27)
% Contact: sevim_cengiz@icloud.com

function [Pinfo]= loadfile(filename,pathname)

[fpath,fname,extfile] = fileparts([pathname,filename]);

if strcmp(extfile,'.spar')

    Pinfo.ext='spar';
    sparfilepath=[fpath,filesep,fname];
    coordfilespath= [fpath,filesep,'lcmodel_outputs',filesep,fname,filesep];
    Pinfo.sparname=fname;
    Pinfo.coordpath=coordfilespath;
    Pinfo.spectrapath=fpath;
    Pinfo.sparfilepath=sparfilepath;
    exam1=strsplit(Pinfo.spectrapath,'exam_1');
    Pinfo.mainpath=exam1{1,1};
    patientno=strsplit(Pinfo.mainpath,filesep);
    Pinfo.name=patientno{1,(end-1)};

    % Get slice, row, column info from .spar file.
    sparheadinfo = textread([pathname,filename], '%s');

    finderinfo=find(ismember(sparheadinfo, 'dim2_pnts')==1); % col  --mm
    Pinfo.ncol= str2double(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'dim3_pnts')==1); % row  --mm
    Pinfo.nrow= str2double(sparheadinfo{finderinfo+2});

    finderinfo=find(ismember(sparheadinfo, 'dim4_pnts')==1); % slice --mm dim4_pnts
    try
        Pinfo.nsli= str2double(sparheadinfo{finderinfo+2});
    catch
        Pinfo.nsli=1;
    end

    finderinfo=find(ismember(sparheadinfo, 'ap_off_center')==1);
    Pinfo.voxoffap = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'lr_off_center')==1);
    Pinfo.voxofflr = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'cc_off_center')==1);
    Pinfo.voxoffcc = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'ap_angulation')==1);
    Pinfo.voxangap= str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'lr_angulation')==1);
    Pinfo.voxanglr = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'cc_angulation')==1);
    Pinfo.voxangcc = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'ap_size')==1); % ap voi size
    Pinfo.apVOI = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'lr_size')==1); % lr voi size
    Pinfo.lrVOI = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'cc_size')==1); % cc voi size
    Pinfo.ccVOI = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'phase_encoding_fov')==1); % FOV --mm
    Pinfo.FOV = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'slice_distance')==1); %slice distance --mm
    Pinfo.slidistance = str2num(sparheadinfo{finderinfo+2});
    finderinfo=find(ismember(sparheadinfo, 'slice_thickness')==1); %slice thickness --mm
    Pinfo.slithickness = str2num(sparheadinfo{finderinfo+2});

    little_vox_dim_lr= (Pinfo.FOV) /  (Pinfo.ncol); % --mm
    little_vox_dim_ap= (Pinfo.FOV) /  (Pinfo.nrow); % --mm
    little_vox_dim_cc= (Pinfo.slithickness) / (Pinfo.nsli); %--mm

    % Calculate each little voxel's location in anatomical position
    clear ncol
    count=1;
    for npoint=0:((Pinfo.ncol/2)-1)
        ncol(1,count)=(-little_vox_dim_lr*npoint)+(Pinfo.voxofflr)-(little_vox_dim_lr/2);
        count=count+1;
    end
    clear pcol
    count=1;
    for npoint=1:((Pinfo.ncol/2))
        pcol(1,count)=(little_vox_dim_lr*npoint)+(Pinfo.voxofflr)-(little_vox_dim_lr/2);
        count=count+1;
    end
    clear loc_col
    loc_col=[pcol ncol];
    loc_col= sort( loc_col);

    clear prow
    count=1;
    for npoint=0:((Pinfo.nrow/2)-1)
        prow(1,count)=(little_vox_dim_ap*npoint)+(Pinfo.voxoffap)+(little_vox_dim_ap/2);
        count=count+1;
    end
    clear nrow
    count=1;
    for npoint=1:((Pinfo.nrow/2))
        nrow(1,count)=(-little_vox_dim_ap*npoint)+(Pinfo.voxoffap)+(little_vox_dim_ap/2);
        count=count+1;
    end
    clear loc_row
    loc_row=[prow nrow];
    loc_row = sort( loc_row);

    clear psli
    count=1;
    for npoint=0:(Pinfo.nsli/2)
        psli(1,count)=(little_vox_dim_cc*npoint)+(Pinfo.voxoffcc);
        count=count+1;
    end
    clear nsli
    count=1;
    for npoint=1:((Pinfo.nsli/2))
        nsli(1,count)=(-little_vox_dim_cc*npoint)+(Pinfo.voxoffcc);
        count=count+1;
    end
    clear loc_sli

    if Pinfo.nsli == 1
        loc_sli=[Pinfo.nsli];
        for isli=1: Pinfo.nsli
            for irow=1:Pinfo.nrow
                for icol=1:Pinfo.ncol
                    Pinfo.littlevoxels(isli,irow,icol).name= sprintf('litvox_sl%d_r%d_c%d',isli,irow, icol);
                    Pinfo.littlevoxels(isli,irow,icol).loc= [loc_sli(1,isli) loc_row(1,irow) loc_col(1,icol)] ;
                end
            end
        end
    else
        loc_sli=[psli nsli];
        loc_sli= sort( loc_sli);

        for isli=1: Pinfo.nsli
            for irow=1:Pinfo.nrow
                for icol=1:Pinfo.ncol
                    Pinfo.littlevoxels(isli,irow,icol).name= sprintf('litvox_sl%d_r%d_c%d',isli,irow, icol);
                    Pinfo.littlevoxels(isli,irow,icol).loc= [loc_sli(1,isli) loc_row(1,irow) loc_col(1,icol)] ;
                end
            end
        end
    end

elseif strcmp(extfile,'.rda')

    [fpath,fname,extfile] = fileparts([pathname,filename]);
    
    rda_filename = fullfile(fpath, [fname, extfile]);
    fid = fopen(rda_filename);
       
    out.sparname = fname;
    coordfilespath = [fpath,filesep,'lcmodel_outputs',filesep,fname,filesep];
    out.coordpath = coordfilespath;
    out.spectrapath = fpath;
    rdafilepath = [fpath,filesep,fname];
    out.sparfilepath = rdafilepath;
    exam1 = strsplit(out.spectrapath,'exam_1');
    out.mainpath = exam1{1,1};
    patientno = strsplit(out.mainpath,filesep);
    out.name = patientno{1,(end-1)};
    rdaheadinfo = textread([pathname,filename], '%s');
     
    head_start_text = '>>> Begin of header <<<';
    head_end_text   = '>>> End of header <<<';
    tline = fgets(fid);

    while (isempty(strfind(tline , head_end_text)))
        tline = fgets(fid);

        if ( isempty(strfind (tline , head_start_text)) + isempty(strfind (tline , head_end_text )) == 2)
            
            occurence_of_colon = findstr(':',tline);
            variable = tline(1:occurence_of_colon-1) ;
            value    = tline(occurence_of_colon+1 : length(tline)) ;

            switch variable
                case { 'PatientID' , 'PatientName' , 'StudyDescription' , 'PatientBirthDate' , 'StudyDate' , 'StudyTime' , 'PatientAge' , 'SeriesDate' , ...
                        'SeriesTime' , 'SeriesDescription' , 'ProtocolName' , 'PatientPosition' , 'ModelName' , 'StationName' , 'InstitutionName' , ...
                        'DeviceSerialNumber', 'InstanceDate' , 'InstanceTime' , 'InstanceComments' , 'SequenceName' , 'SequenceDescription' , 'Nucleus' ,...
                        'TransmitCoil' }
                    eval(['rda.' , variable , ' = value; ']);
                case { 'PatientSex' }
                    switch value
                        case 0
                            rda.sex = 'Unknown';
                        case 1
                            rda.sex = 'Male';
                        case 2
                            rda.sex = 'Female';
                    end

                case {  'SeriesNumber' , 'InstanceNumber' , 'AcquisitionNumber' , 'NumOfPhaseEncodingSteps' , 'NumberOfRows' , 'NumberOfColumns' , 'VectorSize' }
                    eval(['rda.' , variable , ' = str2num(value); ']);
                case { 'PatientWeight' , 'TR' , 'TE' , 'TM' , 'DwellTime' , 'NumberOfAverages' , 'MRFrequency' , 'MagneticFieldStrength' , 'FlipAngle' , ...
                        'SliceThickness' ,  'FoVHeight' , 'FoVWidth' , 'PercentOfRectFoV' , 'PixelSpacingRow' , 'PixelSpacingCol'}
                    value = strrep(value, ',', '.');
                    eval(['rda.' , variable , ' = str2num(value); ']);
                    
                % Extract orientation vectors properly
                case {'VOINormalSag'}
                    value = strrep(value, ',', '.');
                    rda.VOINormalSag = str2num(value);
                case {'VOINormalCor'}
                    value = strrep(value, ',', '.');
                    rda.VOINormalCor = str2num(value);
                case {'VOINormalTra'}
                    value = strrep(value, ',', '.');
                    rda.VOINormalTra = str2num(value);
                case {'VOIRotationInPlane'}
                    value = strrep(value, ',', '.');
                    rda.VOIRotationInPlane = str2num(value);
                    
                case {'SoftwareVersion[0]' }
                    rda.software_version = value;
                case {'CSIMatrixSize[0]' }
                    rda.CSIMatrix_Size(1) = str2num(value);
                case {'CSIMatrixSize[1]' }
                    rda.CSIMatrix_Size(2) = str2num(value);
                case {'CSIMatrixSize[2]' }
                    rda.CSIMatrix_Size(3) = str2num(value);
                case {'PositionVector[0]' }
                    value = strrep(value, ',', '.');
                    out.PositionVector(1) = str2num(value);
                case {'PositionVector[1]' }
                    value = strrep(value, ',', '.');
                    out.PositionVector(2) = str2num(value);
                case {'PositionVector[2]' }
                    value = strrep(value, ',', '.');
                    out.PositionVector(3) = str2num(value);
                case {'RowVector[0]' }
                    value = strrep(value, ',', '.');
                    out.RowVector(1) = str2num(value);
                case {'RowVector[1]' }
                    value = strrep(value, ',', '.');
                    out.RowVector(2) = str2num(value);
                case {'RowVector[2]' }
                    value = strrep(value, ',', '.');
                    out.RowVector(3) = str2num(value);
                case {'ColumnVector[0]' }
                    value = strrep(value, ',', '.');
                    out.ColumnVector(1) = str2num(value);
                case {'ColumnVector[1]' }
                    value = strrep(value, ',', '.');
                    out.ColumnVector(2) = str2num(value);
                case {'ColumnVector[2]' }
                    value = strrep(value, ',', '.');
                    out.ColumnVector(3) = str2num(value);

                otherwise
                    % Unrecognised variable
            end
        else
            % Header start/end markers
        end
    end

    % Read spectral data
    bytes_per_point = 16;
    complex_data = fread(fid , rda.CSIMatrix_Size(1) * rda.CSIMatrix_Size(1) *rda.CSIMatrix_Size(1) *rda.VectorSize * 2 , 'double');
    fclose(fid);

    % Process spectral data
    hmm = reshape(complex_data,  2 , rda.VectorSize , rda.CSIMatrix_Size(1) ,  rda.CSIMatrix_Size(2) ,  rda.CSIMatrix_Size(3) );
    fids = complex(hmm(1,:,:,:,:),hmm(2,:,:,:,:));
    for rows=1:rda.CSIMatrix_Size(1)
        for cols=1:rda.CSIMatrix_Size(2)
            fids(:,:,rows,cols)=fids(:,:,rows,cols)';
        end
    end
    specs = fftshift(fft(fids,[],2),2);

    % Calculate acquisition parameters
    sz = size(specs);
    dwelltime = rda.DwellTime/1000000;
    spectralwidth = 1/dwelltime;
    txfrq = rda.MRFrequency*1000000;
    dims.t = 1;
    dims.subSpecs = 0;
    Bo = rda.MagneticFieldStrength;
    rawAverages = rda.NumberOfAverages;
    averages = 1;
    subspecs = 1;
    rawSubspecs = 'na';
    date = rda.StudyDate;
    seq = rda.SequenceDescription;
    TE = rda.TE;
    TR = rda.TR;
    pointsToLeftShift = 'N/A';

    % Calculate frequency axis
    f = [(-spectralwidth/2)+(spectralwidth/(2*sz(1))):spectralwidth/(sz(1)):(spectralwidth/2)-(spectralwidth/(2*sz(1)))];
    ppm = -f/(Bo*42.577);
    ppm = ppm+4.6082;
    t = [0:dwelltime:(sz(1)-1)*dwelltime];

    % Fill output structure
    out.fids = fids;
    out.specs = specs;
    out.sz = sz;
    out.ppm = ppm;
    out.t = t;
    out.spectralwidth = spectralwidth;
    out.dwelltime = dwelltime;
    out.txfrq = txfrq;
    out.date = date;
    out.dims = dims;
    out.Bo = Bo;
    out.averages = averages;
    out.rawAverages = rawAverages;
    out.subspecs = subspecs;
    out.rawSubspecs = rawSubspecs;
    out.seq = seq;
    out.te = TE;
    out.tr = TR;
    out.pointsToLeftshift = pointsToLeftShift;
    out.Header = rda;
    out.slithickness = rda.SliceThickness;
    out.nrow = rda.CSIMatrix_Size(1);
    out.ncol = rda.CSIMatrix_Size(2);
    out.nsli = rda.CSIMatrix_Size(3);
    out.FOV = rda.CSIMatrix_Size(1)*rda.PixelSpacingRow;

    % Extract voxel offsets from header
    finderinfo = find(ismember(rdaheadinfo, 'MidSlabPosition[1,1]:')==1);
    out.voxoffap = str2num(rdaheadinfo{finderinfo+1});
    finderinfo = find(ismember(rdaheadinfo, 'MidSlabPosition[0,0]:')==1);
    out.voxofflr = str2num(rdaheadinfo{finderinfo+1});
    finderinfo = find(ismember(rdaheadinfo, 'MidSlabPosition[2,2]:')==1);
    out.voxoffcc = str2num(rdaheadinfo{finderinfo+1});

    % Extract voxel dimensions
    finderinfo = find(ismember(rdaheadinfo, 'SlabThickness[2]:')==1);
    out.apVOI = str2num(rdaheadinfo{finderinfo+1});
    finderinfo = find(ismember(rdaheadinfo, 'SlabThickness[1]:')==1);
    out.lrVOI = str2num(rdaheadinfo{finderinfo+1});
    finderinfo = find(ismember(rdaheadinfo, 'SlabThickness[0]:')==1);
    out.ccVOI = str2num(rdaheadinfo{finderinfo+1});

    % ============== FIXED: EXTRACT NORMAL VECTORS ==============
  
    if isfield(rda, 'VOINormalSag') && isfield(rda, 'VOINormalCor') && isfield(rda, 'VOINormalTra')
        
        out.NormSag = rda.VOINormalSag;
        out.NormCor = rda.VOINormalCor;
        out.NormTra = rda.VOINormalTra;
        fprintf('Using VOINormal fields from RDA header\n');
    else
        % Calculate from RowVector and ColumnVector if not available
        if isfield(out, 'RowVector') && isfield(out, 'ColumnVector')
            fprintf('VOINormal not in header, calculating from ColumnVector x RowVector\n');
            % FIXED: REVERSED CROSS PRODUCT ORDER for CSI data
            normalVector = cross(out.ColumnVector, out.RowVector);
            normalVector = normalVector / norm(normalVector);
            % FIXED: No sign negation needed with reversed cross product
            out.NormSag = normalVector(1);
            out.NormCor = normalVector(2);
            out.NormTra = normalVector(3);
        else
            warning('Cannot determine voxel orientation! Using default (transversal)');
            out.NormSag = 0;
            out.NormCor = 0;
            out.NormTra = 1;
        end
    end
    
    % Extract or default in-plane rotation
    if isfield(rda, 'VOIRotationInPlane')
        out.VoI_InPlaneRot = rda.VOIRotationInPlane;
    else
        out.VoI_InPlaneRot = 0;
    end
    
    % Copy to Pinfo
    Pinfo = out;
    Pinfo.ext = 'rda';
    
    % Store orientation parameters in Pinfo
    Pinfo.NormSag = out.NormSag;
    Pinfo.NormCor = out.NormCor;
    Pinfo.NormTra = out.NormTra;
    Pinfo.VoI_InPlaneRot = out.VoI_InPlaneRot;

    % Set flags
    out.flags.writtentostruct = 1;
    out.flags.gotparams = 1;
    out.flags.leftshifted = 0;
    out.flags.filtered = 0;
    out.flags.zeropadded = 0;
    out.flags.freqcorrected = 0;
    out.flags.phasecorrected = 0;
    out.flags.averaged = 1;
    out.flags.addedrcvrs = 1;
    out.flags.subtracted = 0;
    out.flags.writtentotext = 0;
    out.flags.downsampled = 0;
    if out.dims.subSpecs==0
        out.flags.isFourSteps = 0;
    else
        out.flags.isFourSteps = (out.sz(out.dims.subSpecs)==4);
    end

    % Calculate voxel grid coordinates
    little_vox_dim_lr = (out.FOV) / (rda.CSIMatrix_Size(1));
    little_vox_dim_ap = (out.FOV) / (rda.CSIMatrix_Size(2));
    little_vox_dim_cc = (rda.SliceThickness) / rda.CSIMatrix_Size(3);

    % LR direction
    clear ncol_rda;
    count = 1;
    for npoint = 0:((rda.CSIMatrix_Size(1)/2)-1)
        ncol_rda(1, count) = (-little_vox_dim_lr * npoint) + (Pinfo.voxofflr) - (little_vox_dim_lr / 2);
        count = count + 1;
    end

    clear pcol_rda;
    count = 1;
    for npoint = 1:((rda.CSIMatrix_Size(1)/2))
        pcol_rda(1, count) = (little_vox_dim_lr * npoint) + (Pinfo.voxofflr) - (little_vox_dim_lr / 2);
        count = count + 1;
    end

    clear loc_col_rda;
    loc_col_rda = [pcol_rda, ncol_rda];
    loc_col_rda = sort(loc_col_rda);

    % AP direction
    clear prow_rda;
    count = 1;
    for npoint = 0:((rda.CSIMatrix_Size(2)/2)-1)
        prow_rda(1, count) = (little_vox_dim_ap * npoint) + (Pinfo.voxoffap) + (little_vox_dim_ap / 2);
        count = count + 1;
    end

    clear nrow_rda;
    count = 1;
    for npoint = 1:((rda.CSIMatrix_Size(2)/2))
        nrow_rda(1, count) = (-little_vox_dim_ap * npoint) + (Pinfo.voxoffap) + (little_vox_dim_ap / 2);
        count = count + 1;
    end

    clear loc_row_rda;
    loc_row_rda = [prow_rda, nrow_rda];
    loc_row_rda = sort(loc_row_rda);

    % FH/CC direction
    clear psli_rda;
    count = 1;
    for npoint = 0:(rda.CSIMatrix_Size(3)/2)
        psli_rda(1, count) = (little_vox_dim_cc * npoint) + (Pinfo.voxoffcc);
        count = count + 1;
    end

    clear nsli_rda;
    count = 1;
    for npoint = 1:((rda.CSIMatrix_Size(3)/2))
        nsli_rda(1, count) = (-little_vox_dim_cc * npoint) + (Pinfo.voxoffcc);
        count = count + 1;
    end

    clear loc_sli_rda;
    if rda.CSIMatrix_Size(3) == 1
        loc_sli_rda = Pinfo.voxoffcc;
        for isli = 1:rda.CSIMatrix_Size(3)
            for irow = 1:rda.CSIMatrix_Size(2)
                for icol = 1:rda.CSIMatrix_Size(1)
                    Pinfo.littlevoxels(isli, irow, icol).name = sprintf('litvox_sl%d_r%d_c%d', isli, irow, icol);
                    Pinfo.littlevoxels(isli, irow, icol).loc = [loc_sli_rda(1, isli), loc_row_rda(1, irow), loc_col_rda(1, icol)];
                end
            end
        end
    else
        loc_sli_rda = [psli_rda, nsli_rda];
        loc_sli_rda = sort(loc_sli_rda);
        
        for isli = 1:rda.CSIMatrix_Size(3)
            for irow = 1:rda.CSIMatrix_Size(2)
                for icol = 1:rda.CSIMatrix_Size(1)
                    Pinfo.littlevoxels(isli, irow, icol).name = sprintf('litvox_sl%d_r%d_c%d', isli, irow, icol);
                    Pinfo.littlevoxels(isli, irow, icol).loc = [loc_sli_rda(1, isli), loc_row_rda(1, irow), loc_col_rda(1, icol)];
                end
            end
        end
    end

else
    % Other file types
end

end