redShirt.m

classdef redShirt < dataRecording
    properties
        %{
        Properties of parent class: dataRecording
        recordingName=[]; %The name of the recording
        recordingDir %Full directory containing the recorded session
        dataFileNames %Name of recording data files
        startDate %Start date (time) of Recording (matlab long format)
        endDate %End date (time) of Recording (matlab long format)
        samplingFrequency %Sampling rate
        recordingDuration_ms %the total duration of the recording in [ms]
        channelNames %a cell array with the names of the channels
        channelNumbers %an array with integer channel numbers
        triggerNames %the names of trigger channels
        dspLowCutFrequency %Low-pass cutoff frequency in the Neuralynx DSP (in raw data)
        dspHighCutFrequency %High-pass cutoff frequency in the Neuralynx DSP (in raw data)
        multifileMode %if multi files were selected
        nRecordings %number of recordings
        %}
%         includeDigitalDataInTriggers=1;    
    end
    properties (Constant, Hidden)
        maxNumberOfDigitalChannels=4;
        defaultRawDataStreamName='VSD_optical_data';
        defaultAnalogDataStreamName='Analog '; %dont erase the space
        defaultTriggerStreamName='Trigger';
        defaultDigitalDataStreamName='Digital';
        fileExtension='tsm';
        defaultLocalDir='C:\Users\Tulip\Documents\Academic\Post-Doc\Experiments'; %Default directory from which search starts
    end
    properties
        fileSize % size of the recording file in bytes
        DataSize % size of primary data in bytes
        dynamicRange % [1x2] vector with min and max values produced by the RedShirt camera
        headerOffset % length in bytes of the header with info on vsd recording (note: Each FITS structure shall consist of an integral 
        % number of FITS blocks which are each 2880 bytes (23040 bits) in length       
    end
    properties (Hidden)
        MCinfo
        getDataConfig
        totalChannels
        totalAnalogChannels
        streamNames
        rawDataStreamNumber
        rawDataInfo
        analogDataStreamNumber
        analogDataInfo
        triggerStreamNumber
        triggerInfo
        digitalDataStreamNumber
        digitalDataInfo
        channelID
        number2ID
        ID2number
        ZeroADValue
        MicrovoltsPerAD
        ZeroADValueAnalog
        MicrovoltsPerADAnalog
        SweepStartTime
        SweepStopTime
        fileOpenStruct
        recordingDurationLocal_ms
        startDateLocal
        endDateLocal
        cumEnd
        cumStart
    end
    
    methods
        
        function obj = readHeaderFiles(obj)
            filename = strcat(obj.recordingDir,obj.recordingName{1});
            info = fitsinfo_ale(filename);
            rec_dim = info.PrimaryData.Size;
            obj.startDate = info.FileModDate;
            obj.endDate = [];
            obj.fileSize = info.FileSize;
            obj.DataSize = info.PrimaryData.DataSize;
            obj.channelNames = cellstr(num2str([1:rec_dim(1)*rec_dim(2)]'))';
            obj.channelNumbers = [1:rec_dim(1)*rec_dim(2)];
            obj.signalBits = info.PrimaryData.Keywords{2,2};
            obj.samplingFrequency = 1/info.PrimaryData.Keywords{11,2};
            obj.recordingDuration_ms = rec_dim(3)*info.PrimaryData.Keywords{11,2}*1000;
            obj.nRecordings = length(filename);
            obj.chLayoutNumbers = reshape(1:rec_dim(1)*rec_dim(2),rec_dim(1),rec_dim(2));
            obj.dynamicRange = [info.PrimaryData.Keywords{14,2} info.PrimaryData.Keywords{13,2}];
            obj.headerOffset = info.PrimaryData.Offset;
        end

        function [data,T_ms]=getData(obj,channels,startTime_ms,window_ms,mode,name,varargin)
            %Extract MCRack recording raw data from file
            %Usage: [V_uV,T_ms]=obj.getData(channels,startTime_ms,window_ms,name);
            %Input : channels - [1xN] a vector with channel numbers as appearing in the data folder files
            %        startTime_ms - a vector [1xN] of start times [ms]. If Inf, returns all time stamps in recording (startTime_ms is not considered)
            %        window_ms - a scalar [1x1] with the window duration [ms].
            %        name - the name of the stream (if not entered, default name is used)
            %        varargin - optional parameters for reading data in entered in the varargin cell array of the function call with the names of
            %        the parameters and their values in Matlab style
            %Output: I - A 3D matrix [nChannels x nTrials x nSamples] with values of light intensity across specified channels and trials
            %        T_ms - A time vector relative to recording start (t=0 at start)
%             if nargin~=5 | nargin ~=6
%                 error('method getData was not used correctly: wrong number of inputs');
%             end

            % default flags for reading in data 
            params.saveDataToFileMode = 1;
            params.DarkFrameSubtraction = 1; 
            params.DarkFrame_Last = 0;
            params.recompose_channels = 0;
            params = omex_read_params(params, varargin{:});
            params = params2logical(params);

            conversionFactor = 1/1000*obj.samplingFrequency; % from ms to units in data points
            startTime_ms = round(startTime_ms*conversionFactor)/conversionFactor;
            window_ms = round(window_ms*conversionFactor)/conversionFactor;
            endTime_ms = startTime_ms+window_ms; %no need to conversion factor
            recordingDuration_ms = round(obj.recordingDuration_ms*conversionFactor)/conversionFactor;
            windowSamples = round(window_ms*conversionFactor);
            nTrials = length(startTime_ms);     
            % here put your code
            cd(obj.recordingDir)
            filename = obj.dataFileNames{1}; %here implement cycling through recording repeats from same region and equal conditions (= nTrials)
            
            % index ordering for de-interleaving (first channel pixel
            % ordering is maintaned, second channel row order is flipped
            % about its middle point, column ordering is preserved)
            ind1 = NaN(1,256*256/2);
            ind2 = NaN(1,256*256/2);
            for c = 1:256
                ind1( [1:128]+128*(c-1) ) = [1:128] + 256*(c-1);
                ind2( [1:128]+128*(c-1) ) = fliplr([129:256] + 256*(c-1));
            end
            deinterleaving_matrix = cat(1,reshape(ind1,128,256),reshape(ind2,128,256));
            pixels_linear_ind = reshape(1:256*256,256,256);
            
            pixels_per_frame = max(obj.channelNumbers);
            if startTime_ms == 0
                first_frame = 1;
            else
                first_frame = round(startTime_ms*conversionFactor);
            end
            bit_depth = obj.signalBits;
            byte_conversion = bit_depth/8;
            offset = obj.headerOffset;
            
            switch mode
                case 'fullframe'
                    data = ones(size(obj.chLayoutNumbers,1),size(obj.chLayoutNumbers,2),max(windowSamples),nTrials,'int16');     
                    fid = fopen(filename,'r');
                    % skip to the end of frame_skip and start reading from the following frame   
                    for t = 1: nTrials
                        fseek(fid, offset + (pixels_per_frame * byte_conversion * (first_frame(t)-1)), 'bof'); 
                        datatemp = fread(fid,[256*256,windowSamples(t)],'*uint16'); 
                        if params.saveDataToFileMode
                           datatemp = reshape(datatemp,256,256,windowSamples(t));
                           datatemp = datatemp([1:2:255 256:-2:2],:,:);
                        else
                           datatemp = reshape(datatemp,256,256,windowSamples(t));
                        end

                        if params.DarkFrameSubtraction
                            if params.DarkFrame_Last
                                tail_bytes = obj.fileSize - obj.DataSize - offset;
                                tail_bytes >= numel(obj.channelNumbers)*2
                                data_pixels = numel(obj.channelNumbers)*obj.recordingDuration_ms * obj.samplingFrequency/1000;
                                position = offset + data_pixels;
                                fseek(fid,position,'bof')
                                disp(['start of dark frames: ',mat2str(ftell(fid))]);
                                disp(['bytes after end of recording: ',num2str(tail_bytes)]);
                                dark_frame = fread(fid,pixels_per_frame,'*uint16');
                                dark_noise = int16(nanstd(double(dark_frame)));
                                DF = reshape(dark_frame,256,256); 
                                DF = DF([1:2:255 256:-2:2],:);
                            else
                                % alternatively select dark frame from
                                % first frame
                                DF = datatemp(:,:,1);
                            end
                        end        
                            datatemp = datatemp - repmat(DF,1,1,size(datatemp,3));
                    end
                    
                    if params.recompose_channels
                        datatemp = cat(1,datatemp(1:128,:,:),flipdim(datatemp(129:end,:,:),1));
                    end
                    data(:,:,:,t) = datatemp; clear datatemp
                    fclose(fid);
                case 'single_pixels' % % % TO BE MODIFIED AND UPDATED % % %
                    fid = fopen(filename,'r');
                    data = ones(numel(channels),nTrials,max(windowSamples),'int16'); 
                    % convenient to look for contiguous pixels (channels)
                    channels_transformation =  deinterleaving_matrix(channels);
                    [channels1, ord] = sort(channels_transformation);
                    
                    diff_ind = diff(channels1);
                    one_pixel_step = diff_ind;
                    one_pixel_step(find(one_pixel_step ~=1)) = 0;
                    tempvector = one_pixel_step;
                    c = 0; jumps = []; continuity = [];
                    while length(tempvector) > 0
                        c = c +1;
                        if c == 1
                            jumps(c) = find(~tempvector,1,'first');
                            continuity(c) = length(1:jumps-1);
                            tempvector = one_pixel_step(jumps(1)+1:end);
                        else
                            if isempty(find(~tempvector,1,'first'))
                                continuity(c) = length(tempvector);
                                break
                            end
                            jumps(c) = jumps(c-1) + find(~tempvector,1,'first');
                            continuity(c) = length(jumps(c-1)+1:jumps(c)-1);
                            tempvector = one_pixel_step(jumps(c)+1:end);
                        end  
                    end
                    % check this equality at the end !!!!!
                    %
                    % numel(channels) == sum(continuity) + numel(jumps) + 1
                    % this is also equivalent to the fact that the 
                    % length of jumps can be at most total number of pixels
                    % minus 1 in that case all continuity values are 0
                    %
                    jumps = jumps(find(jumps < length(diff_ind)));
                    n_reads = length(jumps) + 1;
                    for r = 1: n_reads
                        if r == 1
                            pixel_indexing = 1:jumps(1);
                            position = channels1(1);
                        elseif r == n_reads
                            pixel_indexing = jumps(r-1)+1:length(diff_ind)+1;
                            position = channels1(pixel_indexing(1));
                        else
                            pixel_indexing = jumps(r-1)+1:jumps(r);
                            position = channels1(pixel_indexing(1));
                        end
                        for t = 1: nTrials
                            current_frame = round( startTime_ms(t) *conversionFactor );
                            fseek(fid, offset + (position * byte_conversion * current_frame) - 1*byte_conversion, 'bof');
                            data(pixel_indexing,t,1:windowSamples(t)) = reshape(fread(fid,[numel(pixel_indexing)*windowSamples(t) 1],'int16=>int16'),...
                                numel(pixel_indexing),1,windowSamples(t));
                        end   
                    end
                    % order pixels as they will be displayed in the final
                    % pixel_linear_order
                    datatemp = zeros(size(data));
                    datatemp(ord,:,:) = data;
                    data = datatemp; clear Itemp
                    fclose(fid);
            end
            
            if obj.convertData2Double
                data = (double(data)) ;
            end       
            if nargout==2
                T_ms = zeros(max(windowSamples),nTrials);
                for t = 1: nTrials
                    T_ms(1:windowSamples(t),t) = [1:windowSamples(t)]*(1e3/obj.samplingFrequency);
                end
            end
        end
        
        function [V_uV,T_ms]=getAnalogData(obj,channels,startTime_ms,window_ms,name)
            %Extract MCRack recording analog data
            %Usage: [V_uV,T_ms]=obj.getAnalogData(channels,startTime_ms,window_ms,name);
            %Input : channels - [1xN] a vector with channel numbers as appearing in the data folder files
            %        startTime_ms - a vector [1xN] of start times [ms]. If Inf, returns all time stamps in recording (startTime_ms is not considered)
            %        window_ms - a scalar [1x1] with the window duration [ms].
            %        name - the name of the stream (if not entered, default name is used)
            %Output: V_us - A 3D matrix [nChannels x nTrials x nSamples] with voltage waveforms across specified channels and trials
            %        T_ms - A time vector relative to recording start (t=0 at start)
            
            %obj.analogDataStreamNumber=find( cellfun(@(x) all(x(1:numel(obj.defaultAnalogDataStreamName))==obj.defaultAnalogDataStreamName),obj.streamNames) );
                %obj.analogDataInfo
            if nargin==4
                obj.getDataConfig.streamname=obj.streamNames{obj.analogDataStreamNumber};
            elseif nargin==5
                obj.getDataConfig.streamname=name;
                %this option should be revised because currently all parameters are derived from the raw data stream
            else
                error('method getData was not used correctly: wrong number of inputs');
            end
            conversionFactor=1/1000*obj.samplingFrequency;
            startTime_ms=round(startTime_ms*conversionFactor)/conversionFactor;
            window_ms=round(window_ms*conversionFactor)/conversionFactor;
            endTime_ms=startTime_ms+window_ms; %no need to conversion factor
            recordingDuration_ms=round(obj.recordingDuration_ms*conversionFactor)/conversionFactor;
            windowSamples=round(window_ms*conversionFactor);
            nTrials=length(startTime_ms);
            V_uV=ones(numel(channels),nTrials,windowSamples,'uint16')*obj.ZeroADValueAnalog;
            
            cumStart=[-Inf obj.cumStart obj.cumEnd(end)];
            cumEnd=[0 obj.cumEnd Inf];
            obj.getDataConfig.StreamNumber=obj.analogDataStreamNumber-1;
            if obj.multifileMode %this mode currently does not support extraction from edges of the recording
                for i=1:nTrials
                    tmpStartTime=startTime_ms(i);
                    startSample=1;
                    
                    pFileStart=find(startTime_ms(i)>=cumStart,1,'last');
                    pFileEnd=find((startTime_ms(i)+window_ms)<=cumEnd,1,'first');
                    
                    for f=pFileStart:pFileEnd
                        tmpEndTime=min([cumEnd(f) endTime_ms(i)]);
                        endSample=round(startSample+(tmpEndTime-tmpStartTime)/1000*obj.samplingFrequency)-1;
                        
                        if f>1 && f<=(obj.nRecordings+1) % data in inside recording range
                            mcstreammex(obj.fileOpenStruct(f-1));
                            obj.getDataConfig.startend=[tmpStartTime;tmpEndTime]-cumStart(f);
                            data=mcstreammex(obj.getDataConfig);
                            data=reshape(data.data,obj.totalAnalogChannels,length(data.data)/obj.totalAnalogChannels);
                            V_uV(:,i,startSample:endSample)=data(channels,:);
                        else % some of the data is outside the recording range - add zeros
                            V_uV(:,i,startSample:endSample)=obj.ZeroADValueAnalog;
                        end
                        startSample=endSample+1;
                        tmpStartTime=tmpEndTime;
                    end
                end
            else
                for i=1:nTrials
                    obj.getDataConfig.startend=[startTime_ms(i);startTime_ms(i)+window_ms];
                    if startTime_ms(i)>=0 && (startTime_ms(i)+window_ms)<=recordingDuration_ms
                        data=mcstreammex(obj.getDataConfig);
                        data=reshape(data.data,obj.totalAnalogChannels,length(data.data)/obj.totalAnalogChannels);
                        V_uV(:,i,:)=data(channels,:);
                    else
                        startSample=min(0,round(startTime_ms(i)*conversionFactor));
                        endSample=min(windowSamples,round((recordingDuration_ms-startTime_ms(i))*conversionFactor)); %end sample in window (not in recroding)
                        obj.getDataConfig.startend=[max(0,startTime_ms(i));min(startTime_ms(i)+window_ms,recordingDuration_ms)];
                        data=mcstreammex(obj.getDataConfig);
                        data=reshape(data.data,obj.totalAnalogChannels,length(data.data)/obj.totalAnalogChannels);
                        V_uV(:,i,1-startSample:endSample)=data(obj.number2ID(channels),:);
                        disp('Recording at edge');
                    end
                end
            end
            if obj.convertData2Double
                V_uV = (double(V_uV) - obj.ZeroADValueAnalog) * obj.MicrovoltsPerADAnalog;
            end
            
            if nargout==2
                T_ms=(1:windowSamples)*(1e3/obj.samplingFrequency);
            end
        end
        function [D,T_ms]=getDigitalData(obj,startTime_ms,window_ms,name)
            %Extract MCRack digital data from recording
            %Usage: [D,T_ms]=getDigitalData(startTime_ms,window_ms,name)
            %Input : channels - [1xN] a vector with channel numbers as appearing in the data folder files
            %        startTime_ms - a vector [1xN] of start times [ms]. If Inf, returns all time stamps in recording (startTime_ms is not considered)
            %        window_ms - a scalar [1x1] with the window duration [ms].
            %        name - the name of the stream (if not entered, default name is used)
            %Output: D - A 3D matrix [nChannels x nTrials x nSamples] with digitalData waveforms across specified channels and trials
            %        T_ms - A time vector relative to recording start (t=0 at start)
            
            if nargin==3
                obj.getDataConfig.streamname=obj.streamNames{obj.digitalDataStreamNumber};
            elseif nargin==4
                obj.getDataConfig.streamname=name;
                %this option should be revised because currently all parameters are derived from the raw data stream
            elseif nargin==1
                startTime_ms=0;
                window_ms=obj.recordingDuration_ms;
                obj.getDataConfig.streamname=obj.streamNames{obj.digitalDataStreamNumber};
            else
                error('method getData was not used correctly: wrong number of inputs');
            end
            conversionFactor=1/1000*obj.samplingFrequency;
            startTime_ms=round(startTime_ms*conversionFactor)/conversionFactor;
            window_ms=round(window_ms*conversionFactor)/conversionFactor;
            endTime_ms=startTime_ms+window_ms; %no need to conversion factor
            recordingDuration_ms=round(obj.recordingDuration_ms*conversionFactor)/conversionFactor;
            windowSamples=round(window_ms*conversionFactor);
            nTrials=length(startTime_ms);
            
            D=false(obj.maxNumberOfDigitalChannels,nTrials,windowSamples); %up to 4 digital bits are allowed
            
            obj.getDataConfig.StreamNumber=obj.digitalDataStreamNumber-1;
            if obj.multifileMode %this mode currently does not support extraction from edges of the recording
                for i=1:nTrials
                    pFileStart=find(startTime_ms(i)>=obj.cumStart,1,'last');
                    pFileEnd=find((startTime_ms(i)+window_ms)<=obj.cumEnd,1,'first');
                    tmpStartTime=startTime_ms(i);
                    startSample=1;
                    for f=pFileStart:pFileEnd
                        mcstreammex(obj.fileOpenStruct(f));
                        
                        tmpEndTime=min([obj.cumEnd(f) endTime_ms(i)]);
                        obj.getDataConfig.startend=[tmpStartTime;tmpEndTime]-obj.cumStart(f);
                        
                        if tmpStartTime>=0 && tmpEndTime<=recordingDuration_ms
                            data=mcstreammex(obj.getDataConfig);
                            endSample=startSample+numel(data.data)-1;
                            D(:,i,startSample:endSample)=rem(floor(data.data*pow2(0:-1:(1-obj.maxNumberOfDigitalChannels))),2)';
                        else
                            error('Requested data is outside stream limits - this is currently not supported in multi file mode');
                        end
                        startSample=endSample+1;
                        tmpStartTime=tmpEndTime;
                    end
                end
            else %single file mode
                for i=1:nTrials
                    obj.getDataConfig.startend=[startTime_ms(i);startTime_ms(i)+window_ms];
                    if startTime_ms(i)>=0 && (startTime_ms(i)+window_ms)<=recordingDuration_ms
                        data=mcstreammex(obj.getDataConfig);
                        D(:,i,:)=rem(floor(data.data*pow2(0:-1:(1-obj.maxNumberOfDigitalChannels))),2)';
                    else
                        startSample=min(0,round(startTime_ms(i)*conversionFactor));
                        endSample=min(windowSamples,round((recordingDuration_ms-startTime_ms)*conversionFactor));
                        obj.getDataConfig.startend=[max(0,startTime_ms(i));min(startTime_ms(i)+window_ms,recordingDuration_ms)];
                        data=mcstreammex(obj.getDataConfig);
                        D(:,i,1-startSample:endSample)=rem(floor(data.data*pow2(0:-1:(1-obj.maxNumberOfDigitalChannels))),2)';
                        disp('Recording at edge');
                    end
                end
            end
            if nargout==2
                T_ms=(1:windowSamples)*(1e3/obj.samplingFrequency);
            end
        end
        function [T_ms]=getTrigger(obj,startTime_ms,window_ms,name)
            %Extract triggers from recording
            %Usage : [T_ms]=obj.getTrigger(,startTime_ms,window_ms,name)
            %Input : startTime_ms - start time [ms].
            %        window_ms - the window duration [ms]. If Inf, returns all time stamps in recording (startTime_ms is not considered)
            %        name
            %Output: T_ms - trigger times [ms] - different triggers are arranged in a cell array
            if isempty(obj.triggerStreamNumber) && isempty(obj.digitalDataStreamNumber) 
                disp('Warming!!!! No trigger or digital data entities in recording');
                T_ms={};
                return;
            end
            if nargin==4
                obj.getDataConfig.streamname=name;
            elseif nargin==3
                obj.getDataConfig.streamname=obj.streamNames(obj.triggerStreamNumber);
            elseif nargin==1
                startTime_ms=0;
                window_ms=obj.recordingDuration_ms;
                obj.getDataConfig.streamname=obj.streamNames(obj.triggerStreamNumber);
            end
            
            endTime_ms=startTime_ms+window_ms; %no need to conversion factor
            
            pFileStart=find(startTime_ms>=obj.cumStart,1,'last');
            pFileEnd=find((startTime_ms+window_ms)<=obj.cumEnd,1,'first');
            nFiles=pFileEnd-pFileStart+1;
            
            nTriggers=numel(obj.triggerStreamNumber);
            nTriggersDigital=2*obj.maxNumberOfDigitalChannels;
            T=cell(nTriggers+nTriggersDigital,nFiles);
            
            tmpStartTime=startTime_ms;
            for i=pFileStart:pFileEnd
                tmpEndTime=min([obj.cumEnd(i) endTime_ms]);
                        
                obj.MCinfo=mcstreammex(obj.fileOpenStruct(i)); %a file can be read only after opening, and if not other file was opened afterwards
                obj.getDataConfig.startend=[tmpStartTime;tmpEndTime]-obj.cumStart(i);
                
                for j=1:nTriggers
                    obj.getDataConfig.StreamNumber=obj.triggerStreamNumber-1;
                    triggerData=mcstreammex(obj.getDataConfig);
                    T{j,i}=triggerData.times+obj.cumStart(i);
                end
                
                if obj.includeDigitalDataInTriggers & ~isempty(obj.digitalDataStreamNumber)
                    [D,Ttmp]=getDigitalData(obj,tmpStartTime,tmpEndTime-tmpStartTime);
                    validChannels=any(any(D,3),2);
                    D=D(validChannels,:,:);
                    for j=1:size(D,1)
                        T{nTriggers+2*j-1,i}=Ttmp(find(diff(squeeze(D(j,:,:)))>0))+obj.cumStart(i);
                        T{nTriggers+2*j,i}=Ttmp(find(diff(squeeze(D(j,:,:)))<0))+obj.cumStart(i);
                    end
                end
                tmpStartTime=tmpEndTime;
            end
            
            for i=1:(nTriggers+nTriggersDigital)
                T_ms{i}=cell2mat(T(i,:));
            end
            T_ms(cellfun(@isempty,T_ms))=[];
            
        end
    end
    
    methods (Hidden)
        %class constructor
        function obj = redShirt(recordingFile)
            %Usage: obj = redShirt(recordingFile),
            %structure of resulting file: [nChannels x nTrials x nSamples]
            
            obj.signalBits = 16; %the quantization of the sampling card
            
            %get data files
            if nargin==0
                recordingFile=[];
            elseif nargin>1
                disp('Object was not constructed since too many parameters were given at construction');
                return;
            end
            obj=obj.getRecordingFiles(recordingFile,obj.fileExtension);
            obj = obj.readHeaderFiles;
        end
    end
end

%Making a new layout
%{
for i=1:252
    tmp=obj.channelNames{i};
    x(i)=tmp(1)-0;
    y(i)=str2num(tmp(2:end));
end
x=x-min(x)+1;
x(x>10)=x(x>10)-1;
x(x==17)=16;

Ena=cell(16,16);
En=nan(16,16);
for i=1:252
    Ena(x(i),y(i))=obj.channelNames(i);
    En(x(i),y(i))=obj.channelNumbers(i);
end

Ena=Ena';
En=En';
%}