function numRecords = UWListBin(TextFileName,OutPutImageFileName,ImageMatrix,E1low, E1high,E2low,E2high,force,applyZoom, rotate,dbgflag)
%USAGE:  numRecords = UWListBin('TextFileName','OutPutImageFileName',ImageMatrix,E1low, E1high,GEFlag,dbgflag);
% 
% TextFileName = is the name befoe the .txt extension for the input file
%
% OutPutImageFileName = name for the binned stack of projection images, it
% will have a .proj extension added.  The code also generates energy
% spectra - both uncorrected and corrected for all events.  The file
% extension is .espec and for each view angle it bins three spectra as
% one array of 512 pixels for uncorrected followed by an array of 512
% pixels for corrected followed by those events inside the energy windows
%.  The binning is 1 KeV/array element, uint16.
%
% ImageMatrix = output image matrix size, N,  to bin to (NxN)
% E1low = lower pha threshold for energy window 1
% E1high = upper PHA threshold for energy window 1
% force = true/flase - force weights of 0 to 1
% applyZoom = true/false = apply zoom factor in .txt file (true) or set to
%             
% rotate = rotation angle for final storage of image stack = rotate*90
% degrees.  Negative and positive numbers accepted.
% dbgflag = true turns on run time debugginb messages to the display
%
% numRecords return = number of records processed by the code.
%
%Author:
%Tom Lewellen, University of Washington, Imaging Research Lab
%2017-02-08 wrote the inital binning code
%
% Revision history after 2017-02-09
% version 0.1 - NOTE this version assumes only a max of two Energy window defined by
% E1low and E1high the lower and upper PHA settings to use.
%
% 13 Feb 2017:    Still need to check why images come out rotated by 90
% degrees from the GE files.   For now we rotate them in the viewing
% programs, but it is an issue that needs to be understood and fixed.
%
% 14 Feb 2017   Added a new output file OutPutImageFileName.param that is a
% small file with two entries, the output image matrix size and the number
% of view angles to be used by the quick look viw tool or similar tool.
%
% 15 Feb 2017:  Added code to get Zoom facotrs and offsets and apply them
% to the final positioning of events in the output file.   Also added code
% to bin up more than one set of energy windows.  The resulting images
% stacks are written in order - that is, energy window1 is the first stack
% immediately followed by a second stack for energyt window2, etc.
%
% 6 March 7 2017 - fixed issue in UWList bin that messed up espectra
% added correct code for zoom apply true/false

numRecords = 0;
countMove = 0;
count = 0;

% first we get some data from the .txt file to allow us to setup needed
% arrays and constants as well as the name of the .DAT file to use.
uwTXT = strcat(TextFileName,'.txt');
imageFile = strcat(OutPutImageFileName,'.proj');
espectFile = strcat(OutPutImageFileName,'.espec');

fidUwText = fopen(uwTXT,'r');
fidImageFile =fopen(imageFile,'w');   % we will store the data in memory and 
                                      % write out the image stack whendone.

 fidEspecFile = fopen(espectFile,'w'); %also stored in memory then written out.



                                      
% for reading in the infomration in our .txt file we know it is a series of
% strings with a new line (\n) as the deliminator.  Matlab will readin the
% entire text file as a string, ignoring the \n characters.

studyInfo = fscanf(fidUwText,'%s');
fclose(fidUwText);

% now we extract what we want by looking for /key/ and the following / 
% for the excersize of testing the extraction code we will read out all the
% parameters and assign them to local varicables even though we will not
% need it all for the output file this program will create.
studyType = extractBetween(studyInfo,'/StudyType/','/');
studyName = extractBetween(studyInfo,'/StudyName/','/');
CT = extractBetween(studyInfo,'/CT/','/');
spectFile = extractBetween(studyInfo,'/SpectFile/','/');
vendor = extractBetween(studyInfo,'/Vendor/','/');
model = extractBetween(studyInfo,'/Model/','/');
collimator = extractBetween(studyInfo,'/Collimator/','/');
mode = extractBetween(studyInfo,'/Mode/','/');
isotope1 = extractBetween(studyInfo,'/Isotope1/','/');
isotope2 = extractBetween(studyInfo,'/Isotope2/','/');
energyUnits = extractBetween(studyInfo,'/EnergyUnits/','/');
numEsets = extractBetween(studyInfo,'/NumEsets/','/');

numberEsets = int16(str2double(numEsets));
for index=1:numberEsets
    stringTemp = strcat('/Energy',num2str(index),'/');
    energyStringTemp = extractBetween(studyInfo,stringTemp,'/');
    
end

correctionsInWindow = extractBetween(studyInfo,'/CorrectionsInWindow/','/');
correctionsOutOfWindow = extractBetween(studyInfo,'/CorrectionsOutOfWindow/','/');
gantryPositionsPerHead = extractBetween(studyInfo,'/GantryPositionsPerHead/','/');
angleRangePerHead = extractBetween(studyInfo,'/AngleRangePerHead/','/');
startAngle = extractBetween(studyInfo,'/StartAngle/','/');
timePerStop = extractBetween(studyInfo,'/TimePerStopInSeconds/','/');
pixelScale = extractBetween(studyInfo,'/PixelScale/','/');
pixelMatrixSize = extractBetween(studyInfo,'/MatrixSize/','/');
detectorMatrixSize = extractBetween(studyInfo,'/DetectorMaxBins/','/');
detectorOrigin = extractBetween(studyInfo,'/DetectorOrigin/','/');
% zoom is the last entry and does not have a following /
zoom = extractBetween(studyInfo,'/ZoomFactor/','/');
xshift = extractBetween(studyInfo,'/Xshift/','/');
yshfit = extractBetween(studyInfo,'/Yshift/','/');
bodyContour = extractAfter(studyInfo,'/BodyContour/');

%now convert some of the strings to numbers needed for the binning
zoomFactor = str2double(zoom);
detScale = str2double(pixelScale);
scaleMatrix = str2double(pixelMatrixSize);
detectorMax = str2double(detectorMatrixSize);
angSep = str2double(mode);
angles = int16(str2double(gantryPositionsPerHead));
detMatSize = str2double(detectorMatrixSize);

% we now save a couple of parameters for use by other software tools
paramFile = strcat(OutPutImageFileName,'.param');
fidParam = fopen(paramFile,'w');
ans = fwrite(fidParam,ImageMatrix,'uint16');
ans = fwrite(fidParam,(angles*2),'uint16');
ans = fclose(fidParam);

% now we define the in memory array of floats (due to event weights) for
% storing data as we process the list mode file.  We will enter it as
% (angle, X, Y) for each energy set.  For this code we assume only 2 energy
% sets.  Will need to find a structure that allows us to make it more
% generic when we use multi primary windows with multipe energies.

    imageArray1 = zeros((angles*2),ImageMatrix,ImageMatrix,'double');
    imageArray2 = zeros((angles*2),ImageMatrix,ImageMatrix,'double');
    especUcArray = zeros((angles*2),512,'uint16');
    especCrArray = zeros((angles*2),512,'uint16');
    especWnArray = zeros((angles*2),512,'uint16');
    
%open source file
fidSpect = fopen(spectFile,'r');

% read first  byte.  We do this here so that later on when we read a record
% byte that is actually past the end of file the while loop will stop
% correctly
recordType=fread(fidSpect,1,'uint8=>uint8');
    
angleIndex = int16(0);   % this will be the index for the angles field in the imageArray


energyScale = 1.0/(str2double(energyUnits));

if applyZoom 
    offset = (detMatSize*(zoomFactor-1.0))/2.0;
else
    offset = 0;
    zoomFactor = 1.0;
end
% for now assume no x or y shfit, just zooming on center


while ~feof(fidSpect)
    if dbgflag
        disp(['record = ',num2str(count)]);  %heartbeat
        count = count + 1;
    end
    
    switch recordType
        case  241 %time record
             
            if dbgflag, disp('time record');end
           
           % for the 0.1 version of this code we ignore the time records,
           % but need to read them in anyway both to check the code and to
           % set the file pointer correctly
           gate =fread(fidSpect,1,'uint8=>uint8');
           timeStamp = fread(fidSpect,1,'uint32=>uint32');
            
        case 240 %event record
            
            if dbgflag, disp('event record');end
           % here we process the event that is in the E window, using its
           % weight and storing it in the output array.
           UncEnergy = fread(fidSpect,1,'uint16=>uint16');
           CorEnergy = fread(fidSpect,1,'uint16=>uint16');
           Detector = int16(fread(fidSpect,1,'uint8=>uint8'));
           Wieght = fread(fidSpect,1,'uint16=>uint16');
           Xdet = double(fread(fidSpect,1,'uint16=>uint32'));
           Ydet = double(fread(fidSpect,1,'uint16=>uint32'));
           
           % convert E into KeV
           CorEnergy = energyScale * double(CorEnergy);
           UncEnergy = energyScale * double(UncEnergy);
           %write out energy spectra data
                UncEnergyIndex = round(UncEnergy);
                CorEnergyIndex = round(CorEnergy);
                if CorEnergyIndex <= 0 
                    CorEnergyIndex = 1;
                end
                if UncEnergyIndex <= 0
                    UncEnergyIndex = 1;
                end
                if CorEnergyIndex > 512
                    CorEnergyIndex = 512;
                end
                if UncEnergyIndex > 512
                    UncEnergyIndex = 512;
                end
                A = angleIndex + (Detector*angles);
                especUcArray(A,UncEnergyIndex) = especUcArray(A,UncEnergyIndex) + 1;
                especCrArray(A,CorEnergyIndex) = especCrArray(A,CorEnergyIndex) + 1;
                 
           
           E1InWindow = ((CorEnergy >= E1low)&& (CorEnergy <= E1high));
           E2InWindow = ((CorEnergy >= E2low)&& (CorEnergy <= E2high));
           
           % is it in the E window?
           if (E1InWindow || E2InWindow)
               % now calculate the weighted event to place in the output array
               % including any shfits and zoom factors.  
               % we apply the zoom and shfit and if the modifed XY events
               % are outside of the specified image matrix they are not
               % saved in the image arrays.
               
               % save windowed energy spectra using energy correcrted data
               especWnArray(A,CorEnergyIndex) = especWnArray(A,CorEnergyIndex) + 1;
               if ((force) && (Wieght == 0))
                    Wieght = 1000;
               end
                
                
                % if no shifts specified, then we use offset, add code for
                % shifts after we debug zooming.
                
                X = round((double(ImageMatrix)*(((Xdet-offset)*zoomFactor)/detMatSize)))+1;
                Y = round((double(ImageMatrix)*(((Ydet-offset)*zoomFactor)/detMatSize)))+1;
                if X >= 129
                    X = 128;
                    Wieght = 0.0;   % if out of bounds, do not store!
                end
                if X < 1
                    X = 1;
                    Wieght = 0.0;
                end
                if Y >= 129 
                    Y = 128;
                    Wieght = 0.0;   % if out of bounds, do not store!
                end
                if Y < 1 
                    Y = 1;
                    Wieght = 0.0;   % if out of bounds, do not store!
                end 
                
                event = double(Wieght)/1000.0;
                %disp(['event weight = ', num2str(event)]);  % put in for
                
           
           % now place in output array based on current gantry angle 
                 
                  % temp debug info
                %if ((X>128)|| (Y>128))
                   %disp([num2str(A),'  ',num2str(X),'   ',num2str(Y),'   ',num2str(event)]);
                    %pause;
                %end
                
                if E1InWindow 
                    imageArray1(A,X,Y) = imageArray1(A,X,Y) + event; 
                end
                
                if E2InWindow 
                    imageArray2(A,X,Y) = imageArray2(A,X,Y) + event; 
                end
                
                
               
           end
           
           
            
        case 242 %motion record
           
            disp(['motion record number ', num2str(countMove)]);   % a sort of hearbeat to mark passage of time
            disp(['file pointer: ',num2str(ftell(fidSpect))]);
            
            countMove = countMove + 1;
            angleIndex = angleIndex + 1;
            frameStart = fread(fidSpect,1,'uint8=>uint8');
            rotationPosition = fread(fidSpect,1,'int32=>int32');
            head1RadialPosition = fread(fidSpect,1,'uint32=>uint32');
            head2RadialPosition = fread(fidSpect,1,'uint32=>uint32');
            TableLinearPosition = fread(fidSpect,1,'uint32=>uint32');
            
            
        otherwise
            disp(['Unknown record type: ',num2str(recordType)]);
            disp(['last file pointer: ',num2str(ftell(fidSpect))]);
            break;
        
     
    end
    
    numRecords = numRecords + 1;
    
    %read next byte to get record type
      recordType=fread(fidSpect,1,'uint8=>uint8');
    
      
end
fclose(fidSpect);
imageArray1 = round(imageArray1);
endArray = angles*2;
for index = 1:endArray
    outArray = imageArray1(index,1:128,1:128);
    if (rotate ~= 0) 
        outArray = rot90(outArray,rotate);
    end
    fwrite(fidImageFile, outArray,'uint16');
end

imageArray2 = round(imageArray2);
endArray = angles*2;
for index = 1:endArray
    outArray = imageArray2(index,1:128,1:128);
    if (rotate ~= 0) 
        outArray = rot90(outArray,rotate);
    end
    fwrite(fidImageFile, outArray,'uint16');
    outArray = especUcArray(index,1:512);
    fwrite(fidEspecFile,outArray,'uint16');
    outArray = especCrArray(index,1:512);
    fwrite(fidEspecFile,outArray,'uint16');
    outArray = especWnArray(index,1:512);
    fwrite(fidEspecFile,outArray,'uint16');
    
end

fclose(fidImageFile);
fclose(fidEspecFile);



   
end





