dcmtab_bids

#!/usr/bin/env perl
package BIDS;
use strict;
use feature 'signatures';
no warnings 'experimental::signatures';
use Carp;
use Data::Dumper;

# these are hard coded translations of modality label suffices to BIDS directory
# can also pass e.g. anat/T1w
my %mode_folder = (dwi=>'dwi',
   T1w=>'anat', T2w=>'anat', MTR=>'anat',
   bold=>'func', sbref=>'func',
   epi=>'fmap',
   'magnitude'=>'fmap', 'magnitude1'=>'fmap', 'magnitude2'=>'fmap',
   'phase'=>'fmap', 'phase1'=>'fmap', 'phase2'=>'fmap', 'phasediff'=>'fmap');

# https://bids-specification.readthedocs.io/en/stable/appendices/entities.html
my @attributes_ordered = qw/task acq ce dir rec run echo mode/;

sub combine_attr($attr, $val){
   # mode (T1w, bold, dwi, ...) is by itself at the end of the file name
   # all others are e.g. run-1, acq-label
   return $attr eq "mode" ? $val : "$attr-$val" if $val;
   return "";
}

# BIDS->new() is only called once we have a protocol to match against a criteria
# %opts is criteria->{bids} -- regex match group of mode, subj, maybe ses and other @attributes_ordered
# NB most "entities" in attributes_ordered will not be provided here.
# only a few are able to be expressed. see $cfgstring in criteria->dsl
sub new($class, %opts){
   # to build out a bids, the very minimum info is subject and mode
   # and if mode == "task" we also need a task name
   my @req = qw/mode subj/;
   push @req, "task" if $opts{mode} =~ m/bold/;
   # only sbref for func needs a task
   push @req, "task" if ($opts{mode} =~ /sbref/ and $opts{folder} =~ /func/);

   for my $req (@req){
      croak "BIDS: must have '$req' in '".Dumper(%opts)."'" if not $opts{$req};
   }

   # given e.g. anat/T1w or perf/asl instead of hard coded suffix
   # add to folder and update mode
   if($opts{mode} =~ m:^(anat|func|fmap|dwi|perf)/([^/]+)$:){
      $opts{folder}=$1;
      $opts{mode}  =$2;
   } else {
      croak "BIDS: bad mode in parsed criteria: ",
             Dumper(\%opts),
        "\nmode='$opts{mode}'! not in @{[keys(%mode_folder)]}"
        if not $mode_folder{$opts{mode}};

      $opts{folder}=$mode_folder{$opts{mode}};
   }

   # check that we have attributes known to be in the BIDS standard
   # e.g. don't try to make up 'mynewkey-value' to put in the file name
   for my $attr (keys(%opts)){
      next if $attr =~ /^(mode|ses|subj|folder)$/;
      croak "BIDS: bad attribute '$attr' not in @attributes_ordered"
        if not grep(/^$attr$/, @attributes_ordered);
   }

   # perl object
   my $self = \%opts;
   bless $self, $class;
   return $self;
}

sub path($self, $ext=".nii.gz") {
   # use BIDS object's subj and ses if they are provided
   my $subj=$self->{subj};
   my $ses =$self->{ses};
   my $prefix="sub-$subj/$self->{folder}/sub-${subj}_";
   $prefix = "sub-$subj/ses-$ses/$self->{folder}/sub-${subj}_ses-${ses}_" if $ses;

   my $file = join "_",
      grep {!/^-?$/}
      map {combine_attr($_, $self->{$_})} @attributes_ordered;
   return "${prefix}$file$ext";
}

package criteria;
use strict;
use feature 'signatures';
use Carp;
use Data::Dumper;
no warnings 'experimental::signatures';

# criteria->new
sub new($class, %opts){
#      (n_expect => 1,
#       bids => {mode=>"T1w"},
#       matches => { pname => qr/mprage/, ndcm => 176,}),
   $opts{n_expect} ||= 1;
   for my $req (qw/bids matches/){
       croak ("criteria: must have '$req' given ". join(", ", keys(%opts))) unless $opts{$req};
   }
   my %subset = %opts{qw/n_expect bids matches/};
   my $self = \%subset;
   bless $self, $class;
}

sub dsl($class, $cfgstring) {
   # mode;pattern;runs
   # bold=mytask;pname=rest,ndcm=180;3
   # bold=rest;pname=rest,ndcm=180
   # T1w;pname=mprage,ndcm=176
   # epi;pname=spinecho,ndcm=4;dir=AP

   # excpliclty check number of seps b/c ; and , is confusing and hard to debug
   # expect at most 2 ';' seps between 3 pieces: mode ; criteria ; #runs
   # but also can have acq= and dir=
   my $help_msg="Want eg. bold=mytask;pname=MyTask,ndcm=180;3;dir=AP;acq=mylabel;run=2";
   my $max_seps = 2;
   my $extra_colons = scalar(()=$cfgstring =~ m/(acq|dir)=/g);
   my $n_seps_seen  = scalar(()=$cfgstring =~ m/;/g);
   croak "bad separator: too many ';' in '$cfgstring'. $help_msg" if $n_seps_seen > ($max_seps + $extra_colons);

   # TODO: support more entities. do not force order (currently dir then acq)
   # maybe write with a grammer
   # consider (runs=)?(?<n_expect>\d+)
   $cfgstring =~ m/^(?<mode>[^;=]+)=?(?<task>[^;]+)?;(?<matches>[^;]+);?(?<n_expect>\d+)?(;dir=(?<dir>[^;=]+))?(;acq=(?<acq>[^;=]+))?(;fixrun=(?<run>\d+))?;?$/ or
     croak "criteria: failed to parse '$cfgstring'. $help_msg";
   my %bids = %+{qw/mode task dir acq run/};
   my $n_expect = $+{n_expect};
   my %matches = map {split /=/} (split /,/, $+{matches});
   return new($class,
              bids=>\%bids,
              matches=>\%matches,
              n_expect=>$n_expect);

}

sub match_seq($self, $seq) {
   my $tests = $self->{matches};
   my $does_match = 0;
   #print "test: ", Dumper($tests), "\n";
   #print "seq ", Dumper($seq), "\n";
   for my $col (keys(%$tests)){
      # TODO: ref($tests->{$col}) eq CODE, Regexp, etc
      #print "# comparing $col: $seq->{$col} ?=? $tests->{$col}\n";
      if($seq->{$col} =~ /$tests->{$col}/){
         $does_match = 1;
      } else {
         return 0;
      }
   }
   return $does_match;
}

package heuristic;
use strict;
use feature 'signatures';
no warnings 'experimental::signatures';

use File::Basename;
use Data::Dumper;
$ENV{DEBUG}="NONE" unless $ENV{DEGUB};

sub bids_out($sequences, @criteria){
   # filter sequences to those with properties matching heuristic criteria
   # check criteria for every session so we can enumerate matches (into run#s)

   #print "ses: ", Dumper($sequences), "\n";
   my %matches = ();
   for my $crt (@criteria) {
      my $prevsesid="";
      my $run;
      for my $seq (@$sequences) {
        my $sesid = join "_", map {$_||""} @{$seq}{qw/subj ses/};

        # 20240401 - if given 'fixrun' will start with run at non-0
        $run = ($crt->{bids}->{run}||0) if $prevsesid ne $sesid;
        $prevsesid = $sesid;

        #print "seq: ", Dumper($seq), "\n";
        if (not $crt->match_seq($seq)){
           print "NO MATCH: ".Dumper($seq).Dumper($crt->{matches}) if $ENV{DEBUG}=~/MISS/;
           next;
        }

        ++$run if $crt->{n_expect} != 1;

        my $nii_out;
        if(ref($crt->{bids}) eq "CODE") {
          $nii_out = $crt->{bids}->($seq, $run);
        } else {
          # will pass run=>$run along with seq info
          # ses will only be passed if it exists
          my @runs=$run?("run" => $run):();
          $nii_out = BIDS->new(%{$crt->{bids}}, @runs, %{$seq}{qw/subj ses/})->path();
        }


        print "# MATCH: $seq->{fullpath} => $nii_out\n" if $ENV{DEBUG}=~/HIT/;

        # using hash so we only write out a file once. later we use values(%results)
        # this then always takes the last sequence if n_expected == 1
        #
        # maybe we want to keep a list of what we are replacing?
        # push @{$matches{$nii_out}}, $cmd;
        $matches{$nii_out}=$seq->{fullpath};
      }
   }

   return %matches;
}

sub parse($sequences, $criteria){
   my %bids = bids_out($sequences, @$criteria);
   #print "cmds: ",Dumper(%cmds), "\n";
   return map {[$_, $bids{$_}]} keys(%bids);
}

package main;
use v5.26;
use strict; use warnings; use autodie;
use feature qw/signatures/; no warnings qw(experimental::signatures);
use File::Slurp; # read_file
use Pod::Usage;
use Data::Dumper;

sub get_criteria(){
   # help for anything that looks like -h --help
   # or nothing at all
   pod2usage(2) if $#ARGV==-1;
   pod2usage(1) if "@ARGV" =~ /-h\>|--help/;

   # everything else should be a criteria
   return map {criteria->dsl($_)} @ARGV;
}

sub read_tsv(){
   my @header;
   my @all;
   while(<STDIN>){
      chomp;
      my @line = split /\t/, $_;
      if(not @header) {
         @header=@line;
	 die "missing subj and/or fullpath in first line (header): '$_'"
	   unless(/subj/ and /fullpath/);
         next;
      }
      my %l;
      @l{@header} = @line;
      push @all, \%l;
   }

   return sort {
      #($a->{subj} cmp $b->{sub} || $a->{subj} <=> $b->{subj}) ||
      $a->{subj} cmp $b->{subj} ||
      ($a->{ses}||0)   <=> ($b->{ses}||0) ||
      ($b->{seqno}||0) <=> ($b->{seqno}||0) } @all
}

# if we aren't testing/using as a module,
# run as command line tool
if(not caller){
   my @criteria = get_criteria;
   my @sequences = read_tsv;
   my @in_out = heuristic::parse(\@sequences, \@criteria);

   print join "\n", map {join "\t", @{$_}} @in_out;
   print "\n";
}

1;

__END__
=head1 NAME

dcmtab_bids

=head1 SYNOPSIS

   dcmdirtab -s 'E\d+' -d 'E07/scans/*/' |
    dcmtab_bids 'bold=rest;ndcm=574,pname=EYESFIXED;acq=slow' 'T1w;pname=mprage,ndcm=176'


Generate commands from dicom info tsv and bids criteria.

Input likely from C<dcmdirtab>. See C<--help> for more.

output ready for C<mknii>. eg:
   dcmtab_bids ... | parallel --colsep '\t' mknii bids/WPC7225_7tclz/"{1}" "{2}"

=head1 OUTPUT

one line per matching pattern with tab separated columns (1) BIDS name (2) matching dicom folder

 sub-1/anat/sub-1_T1w.nii.gz /path/to/dicom/MP2RAGEPTX_TR6000_1MMISO_UNI-DEN_0025



=head1 INPUT

Tab separated fields are read from stdin. Input is expected to be from C<dcmdirtab>. The first row is assumed to be the column labels. Patterns match against these line by line. Matches are translated to BIDS output file names.

=head1 OPTIONS

=over 4

=item -h --help:

this message

=item ENV

Environmental variable 'DEBUG' can be used to see criteria comparisons

=over 6

=item DEBUG=HIT

=item DEBUG=MISS

=item DEBUG=HITMISS

=back

for example:

  dcmdirtab -s 'E\d+' -d 'E07/scans/*/' |
  DEBUG=HIT dcmtab_bids 'bold=rest;ndcm=574;pname=EYESFIXED' 'T1w;pname=mprage,ndcm=176'


=item criteria

'mode;pattern,pattern[;runs][;acq=label][;dir=label]'

Each arguments is interpreted as a disticnt critera specification for parsing dicom metadata to BIDS file structure.

Provide as many as criteria as needed.  Quote each as a single argument.

Criteria components are separated by ';'

Example criteria specifications

=over 4

   bold=mytask;pname=rest,ndcm=180;3
   bold=rest;pname=rest,ndcm=180
   sbref=rest;dname=rest.*sbref,ndcm=1
   T1w;pname=mprage,ndcm=176
   magnitude1;pname=FieldMap,ndcm=144
   phasediff;pname=FieldMap,ndcm=72
   epi;pname=SpinEchoAP,ndcm=4;dir=AP
   perf/asl;pname=asl,ndcm=100
   MTR;pname=MTR,ndcm=38;acq=NM
   anat/UNIT1;ndcm=172,dname=UNI-DEN

=back


=over 6

=item mode

mode can be like 'perf/asl' or just the modality ('T1w') for well known modality labels (T1w, bold)

use 'bold=mytask' to add task-mytask to output filesnames

=item pattern,pattern

Patterns match against the dicom metadata streamin from the input.
You can have as many patterns as needed. separate with a comma

Despite the '=' syntax, the compareison is not strictly equal but "containing": the right side is used as a regular expression to match the field on the left side.

some examples:

=over 8

=item pname=mprage

matches all protocols that contain 'mprage'

=item ndcm=180|246

matches ndcm field containing 180 or 246

=item dname=256x256

matches dname (dicom directory name) containing FOV 256x256


=item tr=1.5

match revolution time of 1.5

=back

=back

run, acq, and dir are optional criteria

=over 6

=item runs

number indicating number of repeats
digit to count the exptect number of matches. each will be enumerated like run-1, run-2, ... in the output filenames

=item acq=label

string label. freeform string to label files with information not inherent provided by the BIDS structure

=item dir

likely dir=AP or dir=PA. useful for 'epi' fmap

=back

When more than one runs, acq, or dir is specified, they must be provided in order: always 'acq=lowres;dir=AP' never 'dir=AP;acq=lowres'


=back


=head1 Testing

pathalogical commitment to a single file
use Perl::RunEND's C<perl-run-end>

=cut

require 'dcmtab_bids';
# in iperl, reply, or pdl2: Module::Refresh->refresh
use Data::Dumper;
use Test2::V0;
use Test::Exception;
is(BIDS::combine_attr("anything",""), "", "combine_attr empty");
is(BIDS::combine_attr("run","1"), "run-1", "combine_attr normal");
is(BIDS::combine_attr("mode","bold"), "bold", "combine_attr mode");
my $sequences = [{subj=>1, pname=>"rest",  ndcm=>180,fullpath=>'xxx/yyyy'  },
                {subj=>2, pname=>"rest",ndcm=>76, fullpath=>'short1'},
                {subj=>2, pname=>"rest",ndcm=>76, fullpath=>'short2'},
                {subj=>3,ses=>1,pname=>"rest",ndcm=>76, fullpath=>'short3'},
                {subj=>3,ses=>2,pname=>"rest",ndcm=>76, fullpath=>'short4'},
                {subj=>3,ses=>2,pname=>"asl",ndcm=>180, fullpath=>'short4'},
                {subj=>3, pname=>"mprage",ndcm=>176,fullpath=>'y/x'       }];

my %c=(n_expect=>1, matches=>{pname => qr/rest/, ndcm=>180}, bids=>{mode=>"bold",task=>"rest"});
my $crt = criteria->new(%c);
ok criteria->new(%c), "make criteria";
ok criteria->dsl('bold=mytask;pname=rest,ndcm=180;3'), "criteria from string";
dies_ok {criteria->dsl('bold=mytask;pname=rest;ndcm=180;3')} "bad separator";
dies_ok {criteria->dsl('bold=mytask,pname=rest;ndcm=180;3;')} "bad separator";
dies_ok {criteria->dsl('bold=mytask;pname=rest;ndcm=180;3;dir=AP')} "bad separator w/dir";

my $cdsl = criteria->dsl('bold=mytask;pname=rest,ndcm=180;3');
is $cdsl->{n_expect}, 3, "criteria dsl nexpect";
is $cdsl->{bids}->{task}, "mytask", "criteria dsl bids:task";
is $cdsl->{matches}->{ndcm}, 180, "criteria dsl matches:ndcm";
$cdsl = criteria->dsl('T1w;pname=mprage');
is $cdsl->{matches}->{pname}, "mprage", "criteria dsl minimal matches:pname";
is $cdsl->{bids}->{mode}, "T1w", "criteria dsl minimal bids:mode";

is $crt->match_seq($sequences->[0]), 1, "heuristic match (rest)";
is $crt->match_seq($sequences->[1]), 0, "heuristic no match (bad rest)";
is $crt->match_seq($sequences->[2]), 0, "heuristic no match (mprage)";
ok heuristic::bids_out($sequences, $crt), "generic sequences cmds";
ok heuristic::parse($sequences,[$crt]), "generic parse";
my @res = heuristic::parse($sequences,[$crt]);
is @res, 1, "parse returns only one";
ok $res[0][0] =~ m:sub-1/.*func:, "output is sub-1 rest xxx/yyyy";
ok $res[0][1] =~ m:xxx/yyyy:, "input is xxx/yyyy";

# criteral->{bids} function instead of {mode=>xxx, ...}
$c{bids} = sub {return "./nonsense.nii.gz"};
my $crtf = criteria->new(%c);
@res = heuristic::parse($sequences,[$crtf]);
is $res[0][0], "./nonsense.nii.gz", "use criteria->{bids} is function instead of template";

# run numbers were previously off
@res = heuristic::parse($sequences,[criteria->dsl('bold=rest;pname=rest,ndcm=76;2')]);
is @res, 4, "found all short rests ";
is scalar(grep {$_->[0] =~ /run-2/} @res), 1, "only 1 run2";
is scalar(grep {$_->[0] =~ /run-1/} @res), 3, "2 run 1s";

$cdsl = criteria->dsl('epi;pname=rest,ndcm=180;dir=AP');
is $cdsl->{n_expect}, 1, "criteria dsl nexpect omited but have dir";
is $cdsl->{bids}->{dir}, "AP", "criteria dsl: direction AP";
is $cdsl->{bids}->{mode}, "epi", "criteria dsl: mode epi";

$cdsl = criteria->dsl('sbref=mytask;pname=rest.*SBRef,ndcm=1');
is $cdsl->{bids}->{mode}, "sbref", "criteria dsl: mode sbref";
is $cdsl->{bids}->{task}, "mytask", "criteria dsl: single band task name";

$cdsl = criteria->dsl('bold=mytask;pname=rest,ndcm=180;dir=PA;acq=xyz;');
is $cdsl->{n_expect}, 1, "criteria dsl nexpect omited but have dir";
is $cdsl->{bids}->{acq}, "xyz", "criteria dsl: acquisition";
is $cdsl->{bids}->{dir}, "PA", "criteria dsl: direction PA";
is $cdsl->{bids}->{mode}, "bold", "criteria dsl: mode bold";

@res = heuristic::parse([$sequences->[0]], [$cdsl]);
#print("$res[0][0]");
is $res[0][0], "sub-1/func/sub-1_task-mytask_acq-xyz_dir-PA_bold.nii.gz",
   "output is sub-1 rest xxx/yyyy with acq and dir";

$cdsl = criteria->dsl('bold=mytask;pname=rest,ndcm=180;dir=PA;fixrun=3');
is $cdsl->{bids}->{run}, "3", "criteria dsl: hardcode run";

$cdsl = criteria->dsl('perf/asl;pname=asl,ndcm=180');
is $cdsl->{bids}->{mode}, "perf/asl", "criteria dsl: mode from folder/mode";
@res = heuristic::parse($sequences,[$cdsl]);
is $res[0][0],  "sub-3/ses-2/perf/sub-3_ses-2_asl.nii.gz", "output is sub-3/perf/.*asl";
done_testing;