File Coverage

blib/lib/App/AFNI/SiemensPhysio.pm

Criterion	Covered	Total	%
statement	55	196	28.0
branch	7	88	7.9
condition	2	37	5.4
subroutine	13	23	56.5
pod	7	14	50.0
total	84	358	23.4

#!/usr/bin/env perl

package App::AFNI::SiemensPhysio;
use strict;
use warnings;
use Carp;
use List::MoreUtils qw/minmax uniq/;
use File::Basename;
use File::Copy 'mv';
use feature 'say';


=pod

=head1 NAME

App::AFNI:SiemensPhysio - Physio from Siemens into format suitable for AFNI's RetroTS retroicor routine

=head1 SYNOPSIS

Get slice based respiration volume per time (RVT) regressors from physio collected on Siemens scanner

  my $p = SiemensPhysio->new({VERB=>1});
  # read MR data (get times, TR, nslices)
  #  looks at all files in this directory with "dicom_hinfo"
  $p->readMRdir('MRRaw/10824_20111108/rest_384x384.21/');

  # read pulse
  $p->readPhysio('10824/20111108/wpc4951_10824_20111108_110811.puls');
  # write card: $protocol_$sessionTime.puls.dat
  $p->writeMRPhys;

  # read card
  $p->readPhysio('10824/20111108/wpc4951_10824_20111108_110811.resp');
  # write resp: $protocol_$sessionTime.resp.dat
  $p->writeMRPhys;

  #
  $p->retroTS('matlab')

  # we could get the raw data
  #   in this case, card was resp was loaded last
  #   thats what will be returned)
  my @pval = $p->getMRPhys();

=head1 DESCRIPTION



=head2 Pipeline

=over

=item read Siemens physio files

=item read timing from MR DICOM files

=item snip physio files relative to MR

=item prepare/run AFNI's RetroTS

=back


=head2 prior art


=over

=item  https://cfn.upenn.edu/aguirre/public/exvolt/

=item https://cfn.upenn.edu/aguirre/wiki/public:pulse-oximetry_during_fmri_scanning

=back

=cut



=head2 new

initialize object

=head3 OPTIONS

=over

=item timetype

MDH (default) or MPCU

=item PhRate

Physio sample rate
defaults to .2

=item pulsStart and respStart
numeric sequence to be remove, also identifies stream type

defaults:

  pulsStart =>'1 2 40 280'
  respStart =>'1 2 20 2'


=item sliceOrder

alt+z (default)
other options: alt-z,seq+z,seq-z,filename # slice order

=item VERB

set to true to be verbose, defaults false

=item trustIdx

don't check sample rate and index count against end-start time
  none=> check both
  MR  => trust MR (TR)
  phys=> trust physio (PhRate as set by init)
  all => trust both

Note: just have to match reg exp, so MRphys is same as all

=back

=cut

sub new {
  my $class = shift;
  my $self  = shift;
  # default to MDH becaues MPCPU doesn't align to MR time
  my %defaults = (
   #run       => 'matlab', # matlab|McRetroTs|none
   timetyp   => 'MDH',    # MPCU|MDH

   PhRate =>'.02',
   # parameters of the acquisition.
   # third is ticktime, but unlcear what the transfor is to get to freq/samplerate
   pulsStart =>'1 2 40 280',
   respStart =>'1 2 20 2',
   sliceOrder => 'alt+z',
   MRDiscardNum => 0, # how many MR volues to discard
   VERB=>0,
   trustIdx=>'none',
  );
  # use defaults when we didn't set anything
  for my $k (keys %defaults) {
    $self->{$k} = $defaults{$k} unless $self and $self->{$k};
  }

  return bless $self, $class;
}

=head2 readPhysio

after intializing p, provide a file name

 $p->readPhysio('10824/20111108/wpc4951_10824_20111108_110811.puls');

=head3 input file format

 1 2 40 280 ... [long space delimn. list of measurements, maybe 5000 for trigger?]
 ECG  Freq Per: 0 0
 PULS Freq Per: 74 807
 RESP Freq Per: 20 2860
 EXT  Freq Per: 0 0
 ECG  Min Max Avg StdDiff: 0 0 0 0
 PULS Min Max Avg StdDiff: 527 1586 828 4
 RESP Min Max Avg StdDiff: 2380 6700 3477 86
 EXT  Min Max Avg StdDiff: 0 0 0 0
 NrTrig NrMP NrArr AcqWin: 0 0 0 0
 LogStartMDHTime:  66439690
 LogStopMDHTime:   71116595
 LogStartMPCUTime: 66439512
 LogStopMPCUTime:  71114802
 6003

=cut

sub readPhysio {
   my $self=shift;
   my $fname=shift;

   croak "cannot open $fname"
     unless $fname and open my $fh, '<', $fname;

   my $pulsStart = $self->{pulsStart};
   my $respStart = $self->{respStart};
   my $timetyp   = $self->{timetyp};

   # first line is physio measures
   # puls and resp have unique start sequences
   my $values = <$fh>;

   croak "$fname does not start with expected resp or pulse prefix sequence values"
     unless $values=~s/^((?$pulsStart)|(?$respStart))\W*//;

   # we can get the type by which regex we matched
   $self->{ptype}  = join('',map { $+{$_}?$_:"" } qw/puls resp/);

   # break values by whitespace, remove 5000 and above
   # 5000 is a scanner trigger, 5003 is end
   $self->{measures} = [ grep { $_ < 5000 } split(/\W+/,$values) ];

   # get settings matching
   my %settings;
   while($_=<$fh>){
    # remove dos chars
    s/
//g;

    # parse settings, primarily for start and end time
    $settings{$1} = $2 if m/(.*):\W+(.*)/;

    # check file integrity; assume last line is 6003
    croak "corrupt file $fname. does not end with 6003"
      if eof and ! m/^6003$/;

   }



   # get start and end from settings
   if(!$settings{"LogStart${timetyp}Time"} or !$settings{"LogStop${timetyp}Time"}) {
     croak "Cannot find Log(Start|Stop)*Time in $fname";
   }

   $self->{physStart} = $settings{"LogStart${timetyp}Time"}/1000;
   $self->{physEnd}   = $settings{"LogStop${timetyp}Time"}/1000;




   # reset rate if its only off by a very small amount
   # and we don't trust the sample rate we provided
   my $newrate = abs($self->{physStart}- $self->{physEnd})/$#{$self->{measures}};
   $self->{PhRate} = $newrate
     if abs($newrate-$self->{PhRate}) < .00001  and
        $self->{trustIdx}!~/All|Phys/i;


   say "file is $self->{ptype} with $#{$self->{measures}} samples, " ,
       "$self->{physStart}s - $self->{physEnd}s (",
       sprintf("%.2f",($self->{physEnd}-$self->{physStart})/60),
       "min), sample rate adjusted to $self->{PhRate}"
     if $self->{VERB};

   # does the time match the sample rate and number of samples
   timeCheck($self->{physStart},
             $self->{physEnd},
             $#{$self->{measures}},
             $self->{PhRate} ) unless $self->{trustIdx}=~/All|Phys/i;
}


=head2 readMRdir

after intializing p, read in MR info from raw DICOM directory

  $p->readMRdir('MRRaw/10824_20111108/rest_384x384.21/');

sets

=over

=item timing (MRstart and MRend)

=item protcol info (protocol,TR,ET,nslices,Series)

=back

=head3 Example Info

dicom header info

  dicom_hdr MRRaw/10824_20111108/rest_384x384.21/MR* |egrep 'protocol|acquisition Time|Echo Time|Repetition Time' -i
    0008 0031       14 [620     ] //                 ID Series Time//164627.359000
    0008 0032       14 [642     ] //            ID Acquisition Time//164932.315000
    0018 0080        4 [1418    ] //            ACQ Repetition Time//1500
    0018 0081        2 [1430    ] //                  ACQ Echo Time//29
    0018 1030        4 [1612    ] //              ACQ Protocol Name//rest
    0019 100a        2 [1788    ] //                               // 29

shortend to

  dicom_hinfo -tag 0008,0032 0008,0031 0018,0080 0018,0081 0018,1030 MR*

=cut

# sets MRstart MRend protocol TR ET protocol nslices Series
sub readMRdir {
 my $self=shift;
 my $dicomdir=shift;
 croak "$dicomdir is not a directory!" if ! -d $dicomdir;
 my $dcmcmd = "dicom_hinfo -tag 0008,0031 0008,0032 0018,0080 0018,0081 0018,1030  0019,100a $dicomdir/*";
 our @returns =     qw/Filename   Series    AcqTime       TR        ET     protocol nslice/;
 # N.B.  nslices/"Number Of Images In Mosaic" (0019,100a) is Siemens specific
 # which is okay, because thats the kind of physio we have


 # the index at which we can find the item we want
 sub getidx {
  my $name=shift;
  return (grep { $returns[$_] eq $name } (0..$#returns))[0];
 }

 # @v is an element for each dcm (line of output)
 my @v=`$dcmcmd` or croak "could not run $dcmcmd";


 # make each line an array
 # so we have an array of arrays
 # v[0] is all info on first dicom
 # v[0][0] is the first dicom's file name
 @v= map { [split / /] } @v;

 # Sort by acq time (just in case file names are out of order)
 # only needed to ensure trunction takes out the right dicoms
 my $ati=getidx('AcqTime');
 @v = sort {$a->[$ati] <=> $b->[$ati]} @v;

 # Truncate dicom list if we have volumes to discard
 # for some multi-band protocols, the first dicom is junk
 # for older protocols, we have to manually discard volumes
 #  -- eitherway we dont want them
 @v=@v[($self->{MRDiscardNum})..$#v] if($self->{MRDiscardNum});

 # record some constant settings/values
 for my $vals (qw/protocol TR ET Series nslice/) {
    my $vidx=getidx($vals);
    # make sure it's constant
    my @allvals = uniq(map {$_->[$vidx]} @v);
    croak "$vals is not constant: ".join(",",@allvals) if $#allvals>0;

    chomp($allvals[0]);
    $self->{$vals} = $allvals[0];
 }
 $self->{nDcms} = $#v+1;
 $self->{TR} /=1000;



 # Acquistion index
 my $ATidx= getidx('AcqTime');

 # find max and min acq time from all MR*s
 my ($starttime, $endtime) = minmax( map {$_->[$ATidx] } @v);


 ## set start and end
 $self->{MRstart} = getMRAcqSecs($starttime);
 $self->{MRend}   = getMRAcqSecs($endtime);

 say "MR starts $self->{MRstart} (DICOM $starttime) and ends $self->{MRend} (DICOM $endtime)"
   if $self->{VERB};

 timeCheck($self->{MRstart},
           $self->{MRend},
           $self->{nDcms},
           $self->{TR}) unless $self->{trustIdx}=~/All|MR/i;
}

=head2 readBIDSJson
read TR MRstart and crate MRend from BIDS style json output
(e.g. created by dcm2niix).

must already have nDcms (number of volumes in 4d) set

=cut
sub readBIDSJson() {
 use JSON qw/decode_json/;
 my $self=shift;
 my $jsonfile=shift;
 open my $JS, '<', $jsonfile or
   croak "could not open BIDS JSON '$jsonfile'";
 my $data = decode_json(do{local $/; <$JS>})  or
   croak "could not parse JSON in '$jsonfile'";

 $self->{'MRstart'} = getMRAcqSecs($data->{'AcquisitionDateTime'});
 # ($+{HH}*60*60) + ($+{MM}*60) + $+{SS} ..
 $self->{TR} = $data->{RepetitionTime}; # in seconds. eg 1.5
 $self->{MRend} = $self->{MRstart} + $self->{nDcms}*$self->{TR};
 # timeCheck is meaningless -- MRend is derived the same way as the check

 # for output, set protocol and series
 $self->{protocol} = $data->{ProtocolName};
 $self->{Series} = $data->{SeriesDescription};
 # Also have meaningful attributes
 # PhaseEncodingDirection SliceTiming
}

=head2 writeMRPhys

write phys during MR to file
works on most recently loaded physio file

  $p->writeMRPhys

=over

=item use getMRphys to get values

=item use writeDat to write values

=back

=cut

sub writeMRPhys {
 my $self=shift;
 my $outfile = shift;
 $outfile="$self->{protocol}_$self->{Series}.$self->{ptype}.dat" if ! $outfile;

 # if we provided a prefix
 if($self->{prefix}){
   # get dir
   my $bn=dirname($self->{prefix});
   croak "prefix directory ($bn) does not exist!" if ! -d $bn;
   $outfile=$self->{prefix}.$outfile;
 }

 $self->{dat}->{$self->{ptype}}=$outfile;

 my @pvals = $self->getMRPhys;

 say "saving to $outfile" if $self->{VERB};
 writeDat($outfile,@pvals)
}




=head2 retroTS

This is kludgy code hacked together and untested :)
=over

=item use a bunch of hacks to find the matlab binary

=item construct a matlab call using options in self object

=item execute matlab or McRetroTS

=item move outputs to fit local naming convention

=back



get/run command to get Resp. Vol./Time (RVT) via AFNI's retroTS a la http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2715870/
MUST have already read MR and written card and resp dat files

  $p->retroTS('matlab')

how this step is handled is defined by the first argument



=over

=item matlab: use RetroTS.m

=item McRetroTs: use compiled matlab verson of RetroTS.m

=item show: print commands for both matlab and McRetroTS

=item none: do nothing (why'd you call me then!?)

=back


=head3 External Commands

see

=over

=item http://afni.nimh.nih.gov/afni/matlab/

=item http://afni.nimh.nih.gov/sscc/dglen/McRetroTS

=item http://fieldtrip.googlecode.com/svn/trunk/external/afni/RetroTS.m

=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;

 # 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};


 # 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};


 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; quit;"

 # 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;
 }


 if($runcmd) {
   # retroTS just dumpes out oba.slibase.1D whereever it is run from
   # ..so run from where we want to save the file
   chdir dirname($outputname);
   system($runcmd);

   # check if we have the expected output
   if(! -e "oba.slibase.1D" ){
     croak "failed to run\n\n: $runcmd\n\n";
   } else {
     # move file to output name
     mv "oba.slibase.1D", $outputname or
        croak "could not move oba.slibase.1D to $outputname";
   }
   print "$outputname # saved as final output " , `date +%F\\ %H:%M` , "\n";
 }


}



#################
#### 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 (172745.487500)
# JSON BIDs AcqTime like 2018-06-15T17:27:45.487500
sub getMRAcqSecs {
  $_=shift;
  m/^(?\d{2})(?\d{2})(?\d{2}\.\d+)$/ or
   m/^\d{4}-\d{2}-\d{2}T(?\d{2}):(?\d{2}):(?\d{2}\.\d+)$/ or
   croak "$_ from MR time does not look like HHMMSS.sssss or like YYYY-MM-DDTHH:MM:SS.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];

 croak "no pvals!? bad indexes?" if $#pval < 0;
 return @pval;
}



# write values to file
sub writeDat {
 my ($outname,@pvals) = @_;
 croak "do not have a file name to save data" if(!$outname);
 croak "have 1 or fewer values to write"      if($#pvals < 1);

 open my $OH, '>', $outname or croak "could not open $outname to save physio";
 print $OH join("\n",@pvals);
 close $OH;
}

#
# get the start and end index of the meat given start,end time of meat and bread
# need: start,end pairs in array ref + size and sample rate
#  sandwichIdx([bread start,end],[meat start, end], size, rate)
#
sub sandwichIdx {
  my ($bread, $meat, $n, $r) = @_;

  # make sure the meat is inside the bread
  my ($MRstart,$MRend)     = @{$meat};
  my ($physStart,$physEnd) = @{$bread};

  croak sprintf("ERROR: MR starts ($MRstart) %.2f min before physio ($physStart)",
		($physStart-$MRstart)/60)
    if( $MRstart  < $physStart );

  croak sprintf("ERROR: MR ends ($MRend) %.02f min after physio ($physEnd)",
		($MRend-$physEnd)/60)
    if( $MRend  > $physEnd);

  # calc start and end by adding sandwitch top (start as ref)
  my $sIdxS = timeToSamples($bread->[0],$meat->[0],$r);
  my $eIdxS = timeToSamples($meat->[0] ,$meat->[1],$r ) + $sIdxS;

  # calc start and end by subt sandwitch bottom (use end as ref)
  my $eIdxE = $n     - timeToSamples($bread->[1],$meat->[1],$r );
  my $sIdxE = $eIdxE - timeToSamples($meat->[0], $meat->[1],$r );

  # are the two the same?
  carp "Inconsistant index calculation. ".
       "Using start time vs end time as ref grabs different sample of measurements\n".
       "(ref=start) $sIdxS to $eIdxS VS $sIdxE to $eIdxE (ref=end) @ $r secs/sample\n"
       if $sIdxS != $sIdxE || $eIdxS != $eIdxE;

  return ($sIdxS,$eIdxS);
}

# the diff of two times over the sampling rate gives the number of samples
sub timeToSamples {
  my ($start, $end, $rate) = @_;
  return  sprintf( '%.0f', abs($end - $start)/$rate ) ;
}

## print whats up
sub sayIndex {
 my $self=shift;

 # print out start and stop index for ps and pe
 my $ps = timeToSamples($self->{physStart},$self->{physStart},$self->{PhRate});
 my $pe = timeToSamples($self->{physStart},$self->{physEnd},$self->{PhRate});
 my $s  = $self->{MRstartIdx} || undef;
 my $e  = $self->{MRendIdx}   || undef;
 # lets talk about what we have
 say "(ps  ) $self->{physStart}  | MR $self->{MRstart} $self->{MRend} | $self->{physEnd} (pe)  ";
 say "(sidx) $ps          | MR $s  $e     | $pe   (eidx)" if $s and $e;

}



1;