# Bruker::Utils package; parsing for Bruker data format
#
# $Id: Utils.pm,v 1.4 2004/11/10 19:39:58 matthewbrett Exp $
package Bruker::Utils;
use Carp;
use File::Basename;
use Math::Matrix 0.4;
require Exporter;
@ISA = qw( Exporter );
@EXPORT = qw( );
@EXPORT_OK = qw(bruker_series bruker_find_dir parse_bruker_params
bruker2generic bruker_transposition
bruker_text_headers best_of);
use strict;
use vars qw($VERSION );
$VERSION = 0.18;
# Bruker text files to load, with error messages
my($subjfile) = [["subject"],
"no subject info - minor problem"];
# Format; file names to try, message if missing, flag=1 if fatal
my(@paramfiles) = (
[["%d/imnd", "%d/method"],
"input params - transformation matrix will fail", 0],
[["%d/acqp"],
"output params - expect problems with size, timing etc", 0],
[["%d/pdata/%d/reco"],
"reco params - cannot do conversion", 1],
[["%d/pdata/%d/d3proc"],
"image display params - used if other files are missing", 0],
);
# option defaults
my(%text_opt_defs) = (
'anon', 1,
);
sub bruker_series{
# gets series nos from Bruker directory
my $pv_dir = shift;
my(@flist, $f, @seriesnos);
@flist = <${pv_dir}/[0-9]*>;
foreach $f(@flist) {
if (-d $f) {
$f =~ /.*?(\d+)$/;
@seriesnos = (@seriesnos, $1);
}
}
return @seriesnos;
}
sub bruker_find_dir{
# does check / find for Bruker directory
my ($pv_start, $one_multf, $options, @rest) = @_;
my @searchpath = @rest ? @rest :
defined($options->{dbpath}) ? @{$options->{dbpath}} :
();
my ($pv_name, $pv_path, $pv_ext, $spath,
@flist, $s_ext, $s_targ, @dirlist);
$pv_start =~ s/\/$//; # strip trailing path sep
($pv_name, $pv_path, $pv_ext) = fileparse($pv_start,'\..*');
if ($pv_path eq './' && $pv_start !~ "^./") { # path not given
# check each search directory, last first
@searchpath = ('.', @searchpath);
} else {
@searchpath = $pv_path;
}
@searchpath = reverse(@searchpath);
# add asterisk for search
$s_ext = $pv_ext eq "" ? "*" : $pv_ext;
SEARCHLOOP: foreach $spath(@searchpath) {
$s_targ = "$spath/$pv_name$s_ext";
print "Searching for $s_targ...\n" if $options->{verbose};
if ($pv_ext ne "") { # we know what we are looking for
if (-d $s_targ) {
push(@dirlist, $s_targ);
last SEARCHLOOP if $one_multf;
}
} else { # use input as root file name
@flist = ();
foreach (glob($s_targ)) {
push(@flist, $_) if (-d);
}
if (@flist) {
@dirlist = (@dirlist, @flist);
last SEARCHLOOP if $one_multf;
}
}
}
return (\@dirlist, \@searchpath);
}
# get header files info
sub bruker_text_headers{
my($pv_dir, $seriesno, $recono, $in_opts) = @_;
my(%options) = %text_opt_defs;
foreach (keys(%$in_opts)) {
$options{$_} = $in_opts->{$_};
}
my($pfile, $file_loc, $bhstr, $warning, $hdrdata, $ptry, $foundf);
local $/;
# add subject info if required
@paramfiles = ($subjfile, @paramfiles) if (!$options{anon});
# initialize file string
$hdrdata = "";
$warning = "";
# cycle through the files passed
PFILE: foreach $pfile(@paramfiles) {
$foundf = 0;
PFILETRY: foreach $ptry(@{$pfile->[0]}) {
# do variable substitution
$file_loc = sprintf($ptry, $seriesno, $recono);
$bhstr = $pv_dir . "/" . $file_loc;
if (open(TD,"<$bhstr") ) {
$foundf = 1;
last PFILETRY;
}
}
unless ($foundf) {
if ($pfile->[2]) { # abort
$warning .= "No $file_loc file; $pfile->[1]; aborting...\n";
$hdrdata = undef;
last PFILE;
} else {
$warning .= "No $file_loc file; $pfile->[1]; continuing...\n";
next PFILE;
}
}
# Read input files as one long record.
{
local $/;
$hdrdata=$hdrdata . <TD>;
}
close TD;
}
# dos2unix strip
$hdrdata =~ s/\r//sg if $hdrdata;
return($hdrdata, $warning);
}
# Parse Bruker parameter file (passed as string) and return as hash
sub parse_bruker_params{
my($hdrdata) = shift;
my($i, $j, $k, $l, $record, $fname, $fval, @arr, @arels, %fields);
# strip out vis sections
$hdrdata =~ s/(^\$\$ \@vis=(.|\n)*?)(^\#\#)/$3/mg;
# to records
foreach $record (split(/\#\#/, $hdrdata)){
($fname, $fval) = split(/=/,$record);
if (defined($fname)) {
if ($fname =~ s/^\$// ) {
chomp($fval);
if ($fval ne "") {
if ($fval =~ s/\((.*)\)\n//) {
@arels = split(/,/, $1);
if ($fval =~ /^<(.*)>$/) {
$fields{$fname} = $1;
}
else {
# deals with 1D or 2D or 3D arrays
@arr = split(/\s+/, $fval);
if (@arels==3){ # 3D
# fills rowwise - is this correct?
$l=0;
for ($i=0;$i<$arels[0];$i++){
for ($j=0;$j<$arels[1];$j++){
for ($k=0;$k<$arels[2];$k++){
$fields{$fname}[$i][$j][$k]=$arr[$l];
$l++;
}
}
}
} elsif (@arels==2) { # 2D
# fills rowwise - is this correct?
$k=0;
for ($i=0;$i<$arels[0];$i++){
for ($j=0;$j<$arels[1];$j++){
$fields{$fname}[$i][$j]=$arr[$k];
$k++;
}
}
} else {
$fields{$fname} = [@arr];
}
}
} else { # not an array
$fields{$fname}=$fval;
}
}
}
}
}
return %fields;
}
# parses Bruker header into a generic format. When adding a format,
# why not add it here, in case your calculations are useful to the
# rest of us? The Bruker header is a bit confusing. Here's some
# terms to start. Each time a subject gets into the scanner and gets
# scanned, is a 'session'. Each individual scan, or time series, is a
# 'series'. A series may be 3d (volume acquisition) or 2d (slice
# acquisition). For each volume (3d) or slice (2d), there may be more
# than one scan acquired. For example, the T2/PD scans are two images
# (for each slice) acquired during the same scan (=series) I will term
# this as a 'layer', so that a normal strucutural has one 'layer',
# whereas the T2/PD scan has two layers. These layers are combined in
# the reconstructed data file (i.e interleaved - volume by volume in
# the 3D case, and slice by slice in the 2D case). The total number
# of images for one volume (slices * layers for 2D, layers for 3D), I
# will term the number of 'images' for that scan. The reconstruction
# can output images with the read/phase or read/slice, phase/slices
# dimensions transposed (this transposition is _not_ reflected in the
# image dimensions from the RECO file, but it _is_ reflected in the
# dimensions in the d3proc file). Each image can be differentially
# transposed. See the paravision documents for more details, and some
# confusing parameter names.
# conversion between Bruker and other datatypes
my(%dataconv) = (
'_8BIT_UNSGN_INT', => {'bytes' => 1,
'd3no' => 2, #?
'min' => 0,
'max', => 255,
'anatype' => 2,
'mnctype' => '-unsigned -byte'},
'_16BIT_SGN_INT' , => {'bytes' => 2,
'd3no' => 3,
'min' => -(2**15),
'max', => (2**15-1),
'anatype' => 4,
'mnctype' => '-signed -short'},
'_32BIT_SGN_INT' , => {'bytes' => 4,
'd3no' => 5,
'min' => -(2**31),
'max', => (2**31-1),
'anatype' => 8,
'mnctype' => '-signed -long'}
);
# RECO transpose coding - see Image_Re.pdf
#1. If (RECO_transposition[i] = 0), no transposition is to be done.
#
#2. If (0 < RECO_transposition[i] < ACQ_dim), the two consecutive data
#directions given by RECO_transposition[i] and RECO_transposition[i] +
#1 will be transposed.
#
#3. If (RECO_transposition[i] = ACQ_dim), the first and final data
#directions will be transposed.
#
#Please note that for ACQ_dim = 2 the two values 1 and 2 are both legal
#for RECO_transposition[i] and will have the same effect; the first and
#second dimensions will be transposed.
sub bruker_transpose {
my($bhdr, $layer) = @_;
my $code = $bhdr->{'RECO_transposition'}[$layer];
return [0, 0, 'none'] if ($code == 0);
my $is2d = ($bhdr->{ACQ_dim} == 2);
return [0, 1, 'read/phase'] if ($is2d || $code == 1);
return [1, 2, 'phase/slice'] if ($code == 2);
return [0, 2, 'read/slice'];
}
# option defaults
my(%opt_defs) = (
'quiet', 1,
'verbose', 0,
'frecofix',1,
'layers', 1
);
sub bruker2generic{
my ($bhdr, $inopts) = @_;
my %options = %opt_defs;
foreach (keys(%$inopts)) {
$options{$_} = $inopts->{$_};
}
my $i4 = Math::Matrix->eye(4);
my $swaps = $i4->new([-1, 0, 0, 0],
[0, -1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]);
my ($key, $code, @t, $tmp, $bruk_trans, @acq_dim);
my ($start, $ft2mm, $trans, $recoff, $so_far, $rotn, $vx2mm, $radio);
my (@ftdim, $is2d, @fov, @vox, $vox, $off, $i, $j, $gradmat, @vector);
my %ghdr = ();
my $i1 = 0;
# get the dimensionality of the acquisition
$ghdr{dimensions} = @{$bhdr->{'RECO_size'}};
if ($ghdr{dimensions} < 2) {
unless ($options{quiet}) {
warn "Do not know how to handle data of less than two dimensions";
}
}
# fill 2D or 3D worth of dimensions etc
$ghdr{dim} = [@{$bhdr->{'RECO_size'}}];
$ghdr{fov} = [@{$bhdr->{'RECO_fov'}}];
# voxel dimensions in mm
foreach $i (0..($ghdr{dimensions}-1)) {
$ghdr{vox}[$i] = $ghdr{fov}[$i] * 10 / $ghdr{dim}[$i];
}
# slice separation mode, see below
$ghdr{sl_mode} = $bhdr->{'ACQ_slice_sepn_mode'};
# test to see if we have an ACQP file
if (defined($ghdr{sl_mode})) {
# We have an ACQ file
if ($ghdr{dimensions} < 3) {
# This is a 2D image
# Therefore we still need slice information
$ghdr{dim}[2] = $bhdr->{'NSLICES'};
$ghdr{sl_thick} = $bhdr->{'ACQ_slice_thick'};
$ghdr{vox}[2]=$bhdr->{'ACQ_slice_thick'};
# slices can be equidistant, contiguous, var parallel
# or var angle. The slice separation as a single no
# only makes sense in the non Var case
if ((($ghdr{sl_mode} eq "Equidistant") ||
($ghdr{sl_mode} eq "Contiguous"))
&& defined($bhdr->{'ACQ_slice_sepn'})) {
$ghdr{sl_sepn} = $bhdr->{'ACQ_slice_sepn'}[0];
$ghdr{vox}[2] = $ghdr{sl_sepn} if ($ghdr{sl_sepn} > 0);
}
} else { # 3D
$ghdr{sl_thick}= $ghdr{vox}[2];
$ghdr{sl_sepn} = $ghdr{sl_thick};
}
if ($ghdr{sl_sepn}) {
$ghdr{sl_gap} = $ghdr{sl_sepn} > 0 ?
$ghdr{sl_sepn} - $ghdr{sl_thick} : 0;
} else {
$ghdr{sl_gap} = $ghdr{sl_thick};
}
# NI is the total number of images (=images) acquired in a
# repetition so, for a 2D volume, and one layer (e.g. EPI),
# this will be equal to the number of slices. If there are
# two layers, NI will be equal to the NSLICES *2. For a
# 3D volume, NSLICES should I think be 1.
$ghdr{image_nr}= $bhdr->{'NI'};
$ghdr{layer_nr} = $ghdr{image_nr} / $bhdr->{'NSLICES'};
# no of volumes
$ghdr{dim}[3]=best_of( $bhdr->{'ACQ_nr_completed'}, 1);
# the TR. This seems to be an array, with one element
# for each multiplex step. We'll take the first I guess.
$ghdr{vox}[3]=$bhdr->{'ACQ_repetition_time'}[0];
# the slice acquisition order. Well in fact the order of
# acquisition of each 2D or 3D image, ie for a 2D scan, with two
# layers and 10 slices, this array will have 20 elements
$ghdr{objorder} = [@{$bhdr->{'ACQ_obj_order'}}];
} else {
# no ACQ, let's try our best, maybe using the d3proc file
# guess no of layers
$ghdr{layer_nr} = 1;
# no of images (as for NI from acqp). These could be different layers,
# could be slices, could be both, we don't know without the ACQP file
$ghdr{image_nr}= @{$bhdr->{'RECO_transposition'}};
if ($ghdr{dimensions} < 3) {
# we don't know how many slices, let's guess that the NI-like parameter
# above actually refers to no of slices
$ghdr{dim}[2] = $ghdr{image_nr};
}
# maybe the total number of volumes can be calculated from the d3proc file
if (!defined($bhdr->{IM_SIZ})) {
# we do seem to have a d3proc
# the following apparently overcomplex expression is to
# get round the situation where the dimensions have been
# transposed, which will be reflected in the d3proc file
my($nslices)=$bhdr->{IM_SIX} * $bhdr->{IM_SIY}* $bhdr->{IM_SIZ}
/ $ghdr{dim}[0] / $ghdr{dim}[1]
* $bhdr->{IM_SIT};
$ghdr{dim}[3]=$nslices / $ghdr{dim}[2];
} else {
# we have to guess at no of volumes
$ghdr{dim}[3]=1;
}
}
# set slice FOV if not already set
if (!defined($ghdr{fov}[2])) {
if (!defined($ghdr{vox}[2])) {
$ghdr{fov}[2] = 1; # no slice volume info
} else {
$ghdr{fov}[2] = $ghdr{vox}[2] * $ghdr{dim}[2] / 10;
}
}
# set use of layers from command line option
if (!$options{layers}) {
if ($ghdr{layer_nr} > 1) {
unless ($options{quiet}) {
warn "Reset layer number to 1 from $ghdr{layer_nr}\n";
}
$ghdr{dim}[$ghdr{dimensions}] *= $ghdr{layer_nr};
$ghdr{layer_nr} = 1;
}
$ghdr{dim}[3] = $options{timelength};
}
# set no of time points from command line option
if ($options{timelength}) {
$ghdr{dim}[3] = $options{timelength};
}
# transposition. For each image, and therefore for each slice
# for a 2D experiment, the image can have its dimensions
# transposed in the output reconstructed image. This is dictated
# by the REC0_transposition parameter, which has one value for
# each layer for a 3D series, and one value for each layer for
# each slice for a 2D series. This transposition is not yet
# reflected in the FOV and image dimensions we have recorded
# above. Calculate the transposed dimensions for each layer and store.
# The paravision parameter is called RECO_transpose_dim in the docs, but not
# in the actual reco files ? version change
my($l, $r, $oimgno);
foreach $l (0..$ghdr{layer_nr}-1) {
$ghdr{transpose}[$l]{code}=$bhdr->{'RECO_transposition'}[$l];
for ($r=$l+$ghdr{layer_nr};$r<$ghdr{image_nr}-1;$r+=$ghdr{layer_nr}){
if ($bhdr->{'RECO_transposition'}[$r] && # maybe we have not enough tranposes
($ghdr{transpose}[$l]{code}!=$bhdr->{'RECO_transposition'}[$r])) {
unless ($options{quiet}) {
warn ("Different transpositions for the same image in ".
"the same volume: confusing, expect ugliness");
}
}
}
# work out the transposed dimensions
$ghdr{transpose}[$l]{dim} = [@{$ghdr{dim}}];
$ghdr{transpose}[$l]{vox} = [@{$ghdr{vox}}];
$ghdr{transpose}[$l]{fov} = [@{$ghdr{fov}}];
$bruk_trans = bruker_transpose($bhdr, $l);
# may need to transpose dimensions
swapels(@{$bruk_trans}[0,1],$ghdr{transpose}[$l]{dim});
swapels(@{$bruk_trans}[0,1],$ghdr{transpose}[$l]{vox});
swapels(@{$bruk_trans}[0,1],$ghdr{transpose}[$l]{fov});
$ghdr{transpose}[$l]{descrip} = $bruk_trans->[2];
# now set the 4x4 transformation matrix
# check that position parameters are defined
# this code needs generalizing
if (best_of(
$bhdr->{IMND_read_offset},
$bhdr->{PVM_ReadOffset},
$bhdr->{PVM_SPackArrReadOffset}) &&
$bhdr->{ACQ_fov}) {
# Compile matrix
$is2d = (@{$bhdr->{ACQ_fov}} == 2);
$start = $i4->clone();
# transposition
$tmp = $start->[$bruk_trans->[0]];
$start->[$bruk_trans->[0]] = $start->[$bruk_trans->[1]];
$start->[$bruk_trans->[1]] = $tmp;
# reco offset
$recoff = $i4->clone();
$recoff->[0][3] = $bhdr->{RECO_offset}[0][$l];
$recoff->[1][3] = $bhdr->{RECO_offset}[1][$l];
if ($is2d) {
$recoff->[2][3] = 0;
} else {
$recoff->[2][3] = $bhdr->{RECO_offset}[2][$l];
}
# ft voxels -> mm
$ft2mm = $i4->clone();
@ftdim = @{$bhdr->{RECO_ft_size}};
@acq_dim = @{$bhdr->{ACQ_fov}};
if ($is2d) {
push(@ftdim, $bhdr->{NSLICES});
push(@acq_dim, $ghdr{fov}[2]);
}
foreach $i(0..2) {
$vox = $acq_dim[$i] / $ftdim[$i] * 10;
$off = ($ftdim[$i]+1)/2 * -$vox;
$ft2mm->[$i][$i] = $vox;
$ft2mm->[$i][3] = $off;
}
# translation offset in read/phase/slice
$trans = $i4->clone();
if (defined($bhdr->{IMND_read_offset})) {
$trans->[0][3] = $bhdr->{IMND_read_offset}[0];
$trans->[1][3] = best_of(
$bhdr->{IMND_phase1_offset}[0],
$bhdr->{ACQ_phase1_offset}[0],
0);
$trans->[2][3] = $bhdr->{IMND_slice_offset};
} else {
$trans->[0][3] = mean(@{
best_of($bhdr->{PVM_ReadOffset},
$bhdr->{PVM_SPackArrReadOffset})
});
$trans->[1][3] = mean(@{
best_of($bhdr->{PVM_Phase1Offset},
$bhdr->{PVM_SPackArrPhase1Offset})
});
$trans->[2][3] = mean(@{
best_of($bhdr->{PVM_SliceOffset},
$bhdr->{PVM_SPackArrSliceOffset})
});
}
# gradient matrix - is transposed
# There is one (3x3) gradient matrix for each slice - see Para_Cla.pdf in the Paravision documentation
# Quoting:
# ACQ_grad_matrix - is a three dimensional array of doubles to describe the
# orientation of the images to be measured.
# � The fast dimension has three items to describe gradient vectors for the Gx/
# Gy/Gz coordinate system. The gradient vectors should also describe a
# unit sphere with a radius of 1.0.
# � The medium dimension has also three items to describe the read-, phase-
# and slice gradient.
# � The slow dimension is used to describe all slices to be measured.
$rotn = $i4->clone();
$gradmat = $bhdr->{ACQ_grad_matrix}[0];
for $i(0..2){
for $j(0..2){
$rotn->[$i][$j] = $gradmat->[$j][$i];
}
}
# whole thing
$vx2mm = $rotn * $trans * $swaps * $ft2mm * $recoff * $start;
# do report
if ($options{showmat}) { # very verbose
$start->print("Tranposition\n");
$so_far = mat_so_far("Reco offset",$recoff, $start);
$so_far = mat_so_far("fft to mm",$ft2mm, $so_far);
$so_far = mat_so_far("Swap xy",$swaps, $so_far);
$so_far = mat_so_far("Translations",$trans, $so_far);
$so_far = mat_so_far("Rotations",$rotn, $so_far);
}
# do radiological conversion if nec
if ($options{radio}) {
$radio = $vx2mm->eye(4);
$radio->[0][0] *= -1;
$vx2mm = $radio * $vx2mm;
if ($options{verbose}) {
print "Converting matrix to radiological orientation\n";
}
}
# and into generic header
$ghdr{mat}[$l] = $vx2mm;
} else {
if ($options{verbose}) {
print "Can't find the data for orientation matrix\n"
}
}
}
# data type
my($datatype);
if ($bhdr->{'RECO_wordtype'}){
$datatype = $bhdr->{'RECO_wordtype'};
} else {
# try and get datatype from d3proc file
if ($bhdr->{DATTYPE}) {
$datatype = $dataconv{$bhdr->{DATTYPE}}{d3no};
}
}
$ghdr{datatype} = {%{$dataconv{$datatype}}};
$ghdr{datatype}{name} = $datatype;
# the size in bytes of the images (potentially) interleaved
$ghdr{layersize} = 1;
foreach $i (0..$ghdr{dimensions}-1) {
$ghdr{layersize} *= $ghdr{dim}[$i];
}
$ghdr{layersize} *= $ghdr{datatype}{bytes};
# endianness of reconstructed data
$ghdr{"endian"} = ($bhdr->{RECO_byte_order} eq 'littleEndian') ?
'le':'be';
# there's an endian problem with FRECO reconstructions
# They report littleEndian, when in fact they are bigEndian
# So, if the raw data is bigEndian, and the reco data is littleEndian,
# then assume this is a FRECO reconstruction, and fix
if ($options{frecofix} &&
($bhdr->{'BYTORDA'} =~ /big/) &&
($ghdr{"endian"} eq "le")) {
if (!$options{quiet}) {
print "Fixing FRECO endian -> bigEndian\n";
}
$ghdr{"endian"} = "be";
}
# some details
if ($bhdr->{'ACQ_time'}){
$bhdr->{'ACQ_time'} =~ /(..:..:..)\s*(\w*\s*\w*\s*\w*)/;
$ghdr{time_of_scan} = $1;
$ghdr{date_of_scan} = $2;
}
# image type
$ghdr{iscomplex}=0;
if ($bhdr->{RECO_image_type}) {
$ghdr{image_type}=$bhdr->{RECO_image_type};
$ghdr{image_type}=~s/_IMAGE//;
if ($ghdr{image_type}=~ /COMPLEX/) {
$ghdr{iscomplex}=1;
}
}
# entire data block max/min values
if ($bhdr->{RECO_maxima}) {
$ghdr{blockmax} = $ghdr{datatype}{min};
foreach $i (@{$bhdr->{RECO_maxima}}) {
$ghdr{blockmax} = $i if ($i > $ghdr{blockmax});
}
}
if ($bhdr->{RECO_minima}) {
$ghdr{blockmin} = $ghdr{datatype}{max};
foreach $i (@{$bhdr->{RECO_minima}}) {
$ghdr{blockmin} = $i if ($i < $ghdr{blockmin});
}
}
return %ghdr;
}
sub swapels {
my($elno1, $elno2, $rarr) = @_;
my($tmp) = $rarr->[$elno1];
$rarr->[$elno1] = $rarr->[$elno2];
$rarr->[$elno2] = $tmp;
}
sub mean {
my($i, $s);
$s = 0;
foreach $i(@_) {
$s += $i;
}
return ($s / @_);
}
sub best_of {
my $a;
A: while (@_) {
$a = shift;
last A if (defined($a));
};
return $a;
}
sub mat_so_far {
my($label, $in, $sofar) = @_;
$in->print("$label\n");
$sofar = $in * $sofar;
$sofar->print("Result so far\n");
return $sofar;
}
1;