0001 function [eexpi,einsp]= find_frc( imgs, ROI, frate, name, ok, remove_pts)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
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
0028 Fseq= fft(seq);
0029
0030
0031
0032
0033 L = round( 0.5*2*lseq/frate );
0034 Fseq([1+L+1:end-L])=0;
0035 Fseq([1,2,end])=0;
0036 seq1= ifft(Fseq);
0037
0038
0039 if std(imag(seq1))>1e-10; error('FFT code'); end
0040 seq1= real(seq1);
0041
0042
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 );
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])= [];
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
0089
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
0124 function x=plotpoints( seq, eexpi, einsp, name, ok, frate);
0125 seq= -seq;
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
0131
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