find_frc

 FIND_FRC: find candidates for FRC
   [einsp,eexpi]= find_frc( imgs, ROI, frate, name, ok, remove_pts)
 INPUTS:
      points= find_frc( imgs, ROI, frate, name)
      ok ==0 (show graph - default)
           1 (show graph - no click) 
           2 no graph
      name = name for graph ('unspecified' otherwise)
      ROI  = Region of interest ( [] means all the image)
 OUTPUTS: 
      einsp - points corresponding to end inspiration
      eexpi - points corresponding to end expiration

 Find candidates for FRC from a time seq of EIT data
 frate is framerate (in fps units)

0001 function [eexpi,einsp]= find_frc( imgs, ROI, frate, name, ok, remove_pts)
0002 % FIND_FRC: find candidates for FRC
0003 %   [einsp,eexpi]= find_frc( imgs, ROI, frate, name, ok, remove_pts)
0004 % INPUTS:
0005 %      points= find_frc( imgs, ROI, frate, name)
0006 %      ok ==0 (show graph - default)
0007 %           1 (show graph - no click)
0008 %           2 no graph
0009 %      name = name for graph ('unspecified' otherwise)
0010 %      ROI  = Region of interest ( [] means all the image)
0011 % OUTPUTS:
0012 %      einsp - points corresponding to end inspiration
0013 %      eexpi - points corresponding to end expiration
0014 %
0015 % Find candidates for FRC from a time seq of EIT data
0016 % frate is framerate (in fps units)
0017 
0018 % $Id: find_frc.m 1932 2009-10-21 19:41:26Z aadler $
0019 
0020 if isempty(ROI); ROI = ones(1,size(imgs.elem_data,1)); end
0021 if nargin <4; name='unspecified'; end
0022 if nargin <5; ok=0; end
0023 if nargin <6; remove_pts=[]; end
0024 
0025 seq= ROI*imgs.elem_data;
0026 lseq = length(seq);
0027 % LPF the signal
0028 Fseq= fft(seq);
0029 
0030 % Cut off freq
0031 % each point is frate/2/len Hz
0032 % want to cut a 0.25Hz = L *frate/2/len; L=CUTOFF *2*len/frate
0033 L = round( 0.5*2*lseq/frate );
0034 Fseq([1+L+1:end-L])=0; %HPF
0035 Fseq([1,2,end])=0;     %LPF
0036 seq1= ifft(Fseq);
0037 
0038 % Test and take real part
0039 if std(imag(seq1))>1e-10; error('FFT code'); end
0040 seq1= real(seq1);
0041 
0042 % Flow calc
0043 flow = diff(seq1);
0044 thresh = median( abs(flow));
0045 
0046 inout= zeros(lseq,1);
0047 i = 1;
0048 inout(i) = sign( flow(i)+eps ); % eps to force not zero
0049 for i = 2:lseq-1
0050   cval= inout(i-1);
0051   if cval*flow(i) > -thresh
0052      inout(i) = cval;
0053   else
0054      inout(i) = -cval;
0055   end
0056 end
0057 inout(lseq) = inout(lseq-1);
0058 
0059 dinout= diff( inout );
0060 fdiff = find( diff(inout) );
0061 fdiff([1,end])= []; % first and last are unreliable
0062 fdiff = setdiff( fdiff, remove_pts);
0063 
0064 eexpi= fdiff( (dinout(fdiff)>0) );
0065 einsp= fdiff( (dinout(fdiff)<0) );
0066 
0067 window= 3; ww= -window:window;
0068 for i=1:length(eexpi);
0069   wind= seq( eexpi(i)+ ww );
0070   ff = find( wind== min(wind) );
0071   ff= ff(1)-window-1;
0072   eexpi(i)= eexpi(i) + ff;
0073 end
0074 for i=1:length(einsp);
0075   wind= seq( einsp(i)+ ww );
0076   ff = find( wind== max(wind) );
0077   ff= ff(1)-window -1;
0078   einsp(i)= einsp(i) + ff;
0079 end
0080 
0081 [einsp,eexpi] = remove_some_points( einsp, eexpi, remove_pts );
0082 
0083 if ok==2; return; end
0084   
0085 clf;
0086 axes('position',[0,.00,.3,0.9]);
0087 imgss= imgs;
0088 %imgss.elem_data= imgs.elem_data(:,[eexpi(1),einsp(1)]);
0089 %show_slices(imgss);
0090  imgss.elem_data= ...
0091      mean(imgs.elem_data(:,eexpi),2) - ...
0092      mean(imgs.elem_data(:,einsp),2);
0093  show_slices(imgss);
0094  axis normal
0095  text(-2,-2,name,'FontSize',16,'FontWeight','Bold', 'Interpreter','none');
0096 
0097 
0098 axes('position',[.3,.08,.7,.9]);
0099 
0100 if ok==1;
0101    plotpoints( seq, eexpi, einsp, name, ok, frate);
0102    ff= find( name=='/'); if any(ff); nname= name(ff(end)+1:end); end
0103    ff= find(nname=='.'); if any(ff); nname=nname(1:ff(end)-1); end
0104    gg= get(gcf,'paperposition');
0105    set(gcf,'paperposition',[gg(1:3),gg(3)*.3]);
0106    if exist('OCTAVE_VERSION'); 
0107       print([nname,'_sig.png'], '-dpng','-S200,75'); 
0108    else;
0109       print('-dpng','-r100',[nname,'_sig.png']); 
0110    end
0111    set(gcf,'paperposition',gg);
0112    return;
0113 end
0114 
0115 while 1;
0116    xpts=plotpoints( seq, eexpi, einsp, name, ok,frate);
0117      
0118    if isempty(xpts); return; end
0119    fprintf('Removing points: '); fprintf('%d, ',xpts); disp(' ');
0120    [einsp,eexpi] = remove_some_points( einsp, eexpi, xpts);
0121 end
0122 
0123 %plot((0:lseq-1)/frate, seq, 'k' );
0124 function x=plotpoints( seq, eexpi, einsp, name, ok, frate);
0125    seq= -seq; % Air (non-conductivity is now the +ve quantity)
0126    plot( (0:length(seq)-1)/frate, seq, 'k' );
0127    hold on;
0128    maxs= max(seq)*1.1;
0129    mins= min(seq)*1.1;
0130 %  plot( [1;1]*eexpi(:)', [mins;maxs]*ones(1,length(eexpi)), 'r')
0131 %  plot( [1;1]*einsp(:)', [mins;maxs]*ones(1,length(einsp)), 'b')
0132    hh= plot( (eexpi-1)/frate, seq(eexpi), 'bo'); set(hh,'LineWidth',4);
0133    hh= plot( (einsp-1)/frate, seq(einsp), 'ro'); set(hh,'LineWidth',4);
0134    hold off;
0135 
0136 if ok==0
0137    title(['(',name,') CLICK ON LINES TO REMOVE: RETURN TO QUIT']);
0138    [x,jnk] = ginput;
0139     x = round(x*frate) + 1;
0140 else 
0141    title(name);
0142 end
0143 
0144 function [einsp,eexpi] = remove_some_points( einsp, eexpi, xpts);
0145    for xx= xpts(:)'
0146       d_insp= abs(xx - einsp); 
0147       d_expi= abs(xx - eexpi); 
0148       if min(d_insp) < min(d_expi)
0149          ff= find(d_insp == min(d_insp));
0150          einsp(ff)= [];
0151       else
0152          ff= find(d_expi == min(d_expi));
0153          eexpi(ff)= [];
0154       end
0155    end