function [PatientData,PRxEventTimes,PTM]=CHARISGUIReportPRxEvents(patientID,thresholdPressure,minEventTime,minPreEventTime,numBars,rawICP, rawTime, PRx,PacketAvg_Time)
% The input patientID must be a string
% Example: '101314101514' will be output as ReportPRxEvent 101314101514.mat
% numBars is in 5 second segments. 1 minute = 12 numBars
%Run PRx functions with inputs
[PatientData] = icpEventFinder(patientID,rawICP,rawTime,thresholdPressure,minEventTime,minPreEventTime);
eventTimes=PatientData.eventTimes';
[ PRxEventTimes , PTM] = eventIcp2eventPrx(eventTimes,rawTime,PacketAvg_Time,PRx);
end
function [PatientData] = icpEventFinder(patientID,rawICP,rawTime,thresholdPressure,minEventTime,minPreEventTime)
%make sure icp data is in mmHg before input
%rawTime must be in seconds
%thresholdPressure must be in mmHg
%minEventTime and minPreEventTime must be in minutes
%eventPreTime and eventTimeLengths are in minutes
%outTime is in seconds
if(length(rawICP)==length(rawTime))
%average every minute, assumes 50 Hz
[meanM,varM,x,~] = statEveryK(rawICP,50*60);
%adjust the mean
lowVarMeanX = find(varM<50);
lowVarMean = meanM(lowVarMeanX);
adjustedMean = interp1(lowVarMeanX,lowVarMean,1:length(meanM));
%find possible event start locations
eventThresholdPressure = thresholdPressure;
highX = find(adjustedMean>=eventThresholdPressure);
DhighX = diff(highX);
highXstarters = [0,find(DhighX > 1)]+1;
highXstartersCount = zeros(1,length(highXstarters));
for starter = 1:length(highXstarters)-1
highXstartersCount(starter) = length(find(DhighX(highXstarters(starter):highXstarters(starter+1)-1)==1));
end
highXstartersCount(end) = length(find(DhighX(highXstarters(end):end)==1));
highXlocations = highX(highXstarters);
%Choose the possible events that last for n minutes or more
nMinutes = minEventTime;
eventLocations = highXlocations(find(highXstartersCount>=nMinutes));
eventTimeLengths = highXstartersCount(find(highXstartersCount>=nMinutes));
%Choose the events that proceed a segment of regular icp
regTimeThresh = minPreEventTime;
regLowThresh = 0;
eventRegTimeLengths = zeros(1,length(eventLocations));
tempSeg = adjustedMean(1:eventLocations(1));
eventRegTimeLengths(1) = backBandCounter(tempSeg(1:end-1),regLowThresh,eventThresholdPressure);
for event = 2:length(eventLocations)
tempSeg = adjustedMean(eventLocations(event-1):eventLocations(event));
eventRegTimeLengths(event) = backBandCounter(tempSeg(1:end-1),regLowThresh,eventThresholdPressure);
end
eventsWithReg = find(eventRegTimeLengths>=regTimeThresh);
eventLocations = eventLocations(eventsWithReg);
eventTimeLengths = eventTimeLengths(eventsWithReg);
outIndex = x(eventLocations);
outTime = rawTime(x(eventLocations));
eventPreTime = eventRegTimeLengths(eventsWithReg);
%Go through each event
goodEvents = ones(length(eventLocations),1);
set(0,'DefaultFigureWindowStyle','docked');
figure;
for event = 1:length(eventLocations)
clc;
domain = x(eventLocations(event)):x(eventLocations(event)+eventTimeLengths(event));
plot(rawTime(domain),rawICP(domain));
hold on;
plot(rawTime(x),adjustedMean,'g');
axis([rawTime(x(eventLocations(event)-eventPreTime(event))) rawTime(x(eventLocations(event)+eventTimeLengths(event))) 0 50])
xlabel('Time in seconds');
ylabel('mmHg');
title(sprintf('Event %d of %d.\n',event,length(eventLocations)));
zoom xon;
pan xon;
hold off;
choice=questdlg('Is this a non-artifact Event?','Find Events','Yes','No','Yes');
switch choice
case 'Yes'
answer=1;
case 'No'
answer=2;
case 'Cancel'
end
if(answer == 2)
goodEvents(event)=0;
end
end
outIndex = outIndex(find(goodEvents==1));
outTime = outTime(find(goodEvents==1));
eventPreTime = eventPreTime(find(goodEvents==1));
eventTimeLengths = eventTimeLengths(find(goodEvents==1));
close
%create eventResult PatientDataect
crThreshold = sprintf('%d mmHg Threshold',thresholdPressure);
crPreTime = sprintf('%d minutes pre-time',minPreEventTime);
crPostTime = sprintf('%d minutes post-time',minEventTime);
criteria = struct('threshold',crThreshold,'preTime',crPreTime,'postTime',crPostTime);
atPreTimes = eventPreTime;
atPostTimes = eventTimeLengths;
eventAttributes = struct('postTimes',atPostTimes,'preTimes',atPreTimes);
PatientData = CHARISicpEventResult(patientID,criteria,outTime,eventAttributes);
else
fprintf('The ICP and time data are of different lengths.\n');
end
end
function [meanM,varM,indexStatStart,indexStatMid] = statEveryK(matrix,k)
[R,C] = size(matrix);
numSeg = floor(R./k);
if(numSeg<1)
numSeg = 1;
end
meanM = zeros(numSeg+1,C);
varM = zeros(numSeg+1,C);
x = 1:numSeg+1;
for c = 1:C
for s = 1:numSeg
fprintf(' Segment %d of %d \n',s,numSeg);
x(s) = (s.*k)-k+1;
tempData = matrix((s.*k)-k+1:s.*k,c);
meanM(s,c) = mean(tempData(find(~isnan(tempData))));
varM(s,c) = var(tempData(find(~isnan(tempData))));
end
if((numSeg)+1<=R)
x(numSeg+1) = numSeg.*k;
tempData = matrix((numSeg.*k)+1:end,c);
meanM(numSeg+1,c) = mean(tempData(find(~isnan(tempData))));
varM(numSeg+1,c) = var(tempData(find(~isnan(tempData))));
end
end
indexStatStart = x;
indexStatMid = x+(round(k/2));
end
function [ count ] = backBandCounter(vector,low,high)
%Counts the number of values that fall within the band from last to first.
%Count stops when a value falls outside the band.
flippedVector = flip(vector);
count = 0;
for k = 1:length(vector);
if(flippedVector(k)>= low && flippedVector(k)<=high)
count = count + 1;
else
break;
end
end
end
function [location]=closestPoint(Vector,GivenPoint)
% Check unique Vector
if length(unique(Vector))==length(Vector);
distance=nan(length(Vector),1);
% Get distances
for i=1:length(Vector)
distance(i,1)=abs(Vector(i)-GivenPoint);
end
[~,location]=min(distance);
else
fprintf('ERROR: Values are not Unique')
location=[];
end
end
function [ PRxEventTimes , PTM] = eventIcp2eventPrx(eventTimes,rawTime,PRxTime,PRx)
PTM = getPRxIcpTimeMatrix(PRx,PRxTime);
PRxEventTimes = nan(length(eventTimes),3);
for c = 1:3;
vector = PTM(:,c);
for event = 1:length(eventTimes);
time = eventTimes(event);
location = closestPoint(vector,time);
if(location == length(find(~isnan(vector))))
PRxEventTimes(event,c) = nan;
else
PRxEventTimes(event,c) = location;
end
end
end
end
function [PTM] = getPRxIcpTimeMatrix(PRx,time)
[R,C] = size(PRx);
PTM = nan(R,C);
time = correctPRxTime(time);
min5Start = length(time)-length(PRx)+1;
min10Start = min5Start+60;
min20Start = min5Start+180;
PTM(:,1) = time(min5Start:end);
PTM(1:end-60,2) = time(min10Start:end);
PTM(1:end-180,3) = time(min20Start:end);
end
function [outTime] = correctPRxTime(inTime)
dTime = diff(inTime);
delta = mode(dTime);
X = find(dTime == delta);
V = inTime(X);
Xq = 1:length(inTime);
outTime = interp1(X,V,Xq)';
if(isnan(outTime(end)))
outTime(end) = outTime(end-1)+5;
end
end