view release on metacpan or  search on metacpan

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

#!/usr/bin/env perl

package App::AFNI::SiemensPhysio;
use strict; 
use warnings;
use Carp;
use List::MoreUtils qw/minmax uniq/;
use File::Basename;
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

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

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',
   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 ...
 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/^((?<puls>$pulsStart)|(?<resp>$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
   $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, ",
       "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;

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