view release on metacpan or  search on metacpan

lib/App/AFNI/SiemensPhysio.pm  view on Meta::CPAN

=back

if using matlab+retroTS, the path to retroTS.m should be in your MATLABPATH
  export MATLABPATH="$HOME/afni_matlab/matlab/:$MATLABPATH"

=cut


sub retroTS {
 my $self=shift;
 my $runtype=shift;
 return if $runtype and $runtype =~ /none/i;

 # default to using matlab
 my $matlabbin="matlab";

 # find where matlab script points to
 my $acutalmatlab=`perl -ne 'print \$& if /^.*matlab /' \`which matlab\``;
 $matlabbin=$acutalmatlab if $acutalmatlab;

 # or use env MATLABBIN
 $matlabbin= $ENV{MATLABBIN} if $ENV{MATLABBIN};


 # we need to have both data types
 croak "need writeMRPhys for both puls and resp" 
   if ! $self->{dat} || ! -e $self->{dat}->{resp} || ! -e $self->{dat}->{resp};

 my %params = (
   "Opts.Respfile"   => "'".$self->{dat}->{resp}."'", # Respiration data file
   "Opts.Cardfile"   => "'".$self->{dat}->{puls}."'", # Cardiac data file
   "Opts.PhysFS"     => 1/$self->{PhRate},    # Physioliogical signal sampling frequency in Hz.
   "Opts.Nslices"    => $self->{nslice},      # Number of slices
   "Opts.VolTR"      => $self->{TR},          # Volume TR in seconds
   "Opts.SliceOrder" => "'".$self->{sliceOrder}."'"  # ['alt+z']/'alt-z'/'seq+z'/'seq-z'/'Custom'/filename.1D
 );

 # McRetroTS Respdatafile ECGdatafile VolTR Nslices SamplingFreq(PhysFS) ShowGraphs
 my @mcrts = qw/Opts.Respfile Opts.Cardfile Opts.VolTR Opts.Nslices Opts.PhysFS/;
 my $mccmd =  "McRetroTs @params{@mcrts}";
 say $mccmd if $runtype !~ /matlab|McRetroTs/g ;;


 # if have matlab and singal toolbox, can use this
 my $cmd = join("; ", map { join("=",$_,$params{$_}) } keys %params);
 $cmd .= "; Opts.ShowGraphs=0;Opts.Quiet=0;"; # turn off graphs, turn on verbose
 $cmd .= " rts = RetroTS(Opts)";

 # we should wrap matlab up in a try+quit so we dont hang in ML command window on a failure
 my $matlabwrap= qq/$matlabbin -nodisplay -r "try; $cmd; catch err; err, exit(1); end; rts, quit;"/;

 say $matlabwrap if $runtype !~ /matlab|McRetroTs/i;
 # eg
 # matlab -nodisplay -r "try; Opts.Cardfile='rest_164627.359000.puls.dat'; Opts.VolTR=1.5; Opts.Nslices=29; Opts.SliceOrder='alt+z'; Opts.PhysFS=50.0074711455304; Opts.Respfile='rest_164627.359000.resp.dat'; rts = RetroTS(Opts); catch; exit(666); end...

 # with either command, the original output name will be "oba.slibase.1D"
 # change that to our basename (assume resp and puls have same basename, use one from resp)
 my $outputname = $self->{dat}->{resp};
 $outputname =~ s/.resp.dat$/.slibase.1D/;
 
 # or rename to specified input
 #my $outputname=shift if $#_;
 
 
 my $runcmd="";
 if($runtype =~ /matlab/i){
  $runcmd=$matlabwrap
 }elsif($runtype =~ /McRetroTs/i){
  $runcmd=$mccmd;
 }

 system($runcmd) if $runcmd;

 if( $runcmd && ! -e "oba.slibase.1D" ){
   croak "failed to run $runcmd";
 } else {
   # move file to output name if told to
   rename "oba.slibase.1D", $outputname if $outputname;
 }
}



#################
#### helpers ####
#################


# do start end and tr times make sense?
sub timeCheck {
 my ($start,$end,$n,$tau) = @_;

 #my $maxDiffSec=$tau*2;
 my $maxDiffSec=$tau;
 my $dur=$end-$start;
 my $ideal=$n *$tau;

 # start and end time are sane
 croak "time starts ($start) before or on end time ($end)!" 
   if($end<=$start);

 # samples * sample rate == actual duration
 my $offby  = sprintf('%.3f',$ideal-$dur);
 my $offbyN = sprintf('%.0f',$offby/$tau);
 croak "off by $offby s ($offbyN samples) >$maxDiffSec diff: $n samples at $tau should be $ideal secs not $dur ($end - $start)" 
   if(abs($offby) > $maxDiffSec);

 return 1;
}

# DICOM Acq Time is fmt like HHMMSS.SS
sub getMRAcqSecs {
  $_=shift;
  m/^(?<HH>\d{2})(?<MM>\d{2})(?<SS>\d{2}\.\d+)$/ 
    or croak "$_ from MR time does not look like HHMMSS.sssss";

  my $secs = ($+{HH}*60*60) + ($+{MM}*60) + $+{SS} ;
  return $secs;
}



# returns vector of phys for the timing of an MR file
sub getMRPhys {
 my $self=shift;
 my ($s,$e) = sandwichIdx( 
                    [@{$self}{qw/physStart physEnd/}], 
                    [@{$self}{qw/MRstart MRend/}], 
                    $#{$self->{measures}},
                    $self->{PhRate} );
 $self->{MRstartIdx} = $s;
 $self->{MRendIdx}   = $e;

 ## print out where data is coming from/how sandwiching worked
 $self->sayIndex if $self->{VERB};

 my @pval = @{$self->{measures}}[$s..$e];