view release on metacpan or  search on metacpan

lib/PDL/IO/Dcm/Plugins/MRSSiemens.pm  view on Meta::CPAN

	use PDL::NiceSlice;
	my $self=shift;
	my $opt=shift;
	# we need these modules, return undef otherwise.
	require PDL::Dims || return;
	require PDL::Transform || return;
	PDL::Dims->import(qw/is_equidistant dmin dmax vals hpar dinc initdim dimsize drot idx diminfo /);
	#say "init_dims: ",$self->hdr->{dcm_key};
	#say "init_dims: ",$self->hdr->{dicom}->{Rows} ;
	#say "init_dims: hpar ",hpar($self,'dicom','Rows');
	PDL::Transform->import(qw/t_linear/);
	#say diminfo ($self),$self->hdr->{ascconv}->{sGroupArray_asGroup_1__nSize};
# center=inner(rot*scale,dim/2)+Image Position (Patient)
# xf=t_linear(matrix=>Pixel Scale*$rot->transpose,post=>Image Position (Patient))
	if (hpar($self,'ascconv','sGroupArray_asGroup_1__nSize')){
		warn "Multiple slice groups, support dubious!";
	}
	say "Dims for scan ",$self->hdr->{dicom}->{"Series Number"};
	my $pos_d=zeroes(3);
	my $fov=zeroes(3);
	$fov(0).=hpar($self,'csa','VoiReadoutFoV');
	$fov(1).=hpar($self,'csa','VoiPhaseFoV');
	$fov(2).=hpar($self,'csa','VoiThickness');
	my $t=hpar($self,'csa','VoiInPlaneRotation');
	#$v.=hpar($self,'csa','sSliceArray_asSlice_0__sPosition_dSag') ||0; #x
	$pos_d.=pdl (hpar($self,'csa','VoiPosition')); # ||0; #y
	my $nv=pdl (hpar($self,'csa','VoiOrientation')); # ||0; #y
	my $av=pdl([[0,-$nv(2),$nv(1)],[$nv(2),0,-$nv(0)],[-$nv(1),$nv(0),0]])->squeeze;
	my $srot=identity(3)*cos($t)+sin($t)*$av+(1-cos($t))*$nv->transpose*$nv;
	my $rot=identity($self->ndims);
	say $srot;
	my $s=zeroes(3);
	# This is apparently not so easy to get right ...
	# May be it has to do with phase encoding direction ...
	# or software version ..
	#$s(0).=(hpar($self,'csa','SpectroscopyAcquisitionPhaseColumns')); # ||0; #y
	#$s(1).=(hpar($self,'csa','SpectroscopyAcquisitionPhaseRows')); # ||0; #y
	#$s(0).=(hpar($self,'dicom','Columns')); # ||0; #y
	#$s(1).=(hpar($self,'dicom','Rows')); # ||0; #y
	$s(0).=hpar($self,'ascconv','sSpecPara_lFinalMatrixSizePhase')||1;
	$s(1).=hpar($self,'ascconv','sSpecPara_lFinalMatrixSizeRead')||1;
	$s(2).=hpar($self,'ascconv','sSpecPara_lFinalMatrixSizeSlice')||1;
	#$s(2).=(hpar($self,'csa','SpectroscopyAcquisitionOut-of-planePhaseSteps')); # ||0; #y
	my $pe=hpar($self,'csa','PhaseEncodingDirection');
# Calculate and initialise the transformation
	my @ors=qw/Sag Cor Tra/;
	say $srot;
	$self->hdr->{orientation}=$ors[maximum_ind(abs($srot(2;-)))];  # normal vector lookup

	if ($self->hdr->{orientation} eq 'Tra'){ # transversal slice
		#say $self->hdr->{orientation}. " Orientation";
		$self->hdr->{sl}='z';
		if ($pe eq 'COL'){
			$self->hdr->{ro}='x'; 
			$self->hdr->{pe}='y';
		} else { 
			$self->hdr->{ro}='y'; 
			$self->hdr->{pe}='x';
		}
	}
	if ($self->hdr->{orientation} eq 'Cor'){ # coronal slice
		$self->hdr->{sl}='y';
		if ($pe eq 'COL') {
			$self->hdr->{ro}='x'; 
			$self->hdr->{pe}='z';
		} else { 
			$self->hdr->{ro}='z'; 
			$self->hdr->{pe}='x';
		}
	}
	if ($self->hdr->{orientation} eq 'Sag'){ # sagittal slice
		$self->hdr->{sl}='x';
		if ($pe eq 'COL') {
			$self->hdr->{ro}='y'; 
			$self->hdr->{pe}='z';
		} else { 
			$self->hdr->{ro}='z'; 
			$self->hdr->{pe}='y';
		}
	}
	# Slice groups! Does that make sense?
	if ($$opt{split}) {
		my ($sg,$size)=PDL::IO::Dcm::Plugins::MRISiemens::map_slicegroup($self);
		$s(2)*=$size;
	} else {
		#$s(2)*=$self->hdr->{csa}->{'NumberOfFrames'};
		$s(2)*=$self->hdr->{ascconv}->{'sSliceArray_lSize'};
	}	
	#$s(2).=1 if ($s(2)<1);
	say "FOV $fov matrix $s";
	my $inc_d=zeroes(3);
	$inc_d=$fov/$s;
	#say $srot;
	$rot(2:4,2:4).=$srot;
	#say "FOV $fov matrix $s, pixels $inc_d";
	barf $self->hdr->{dicom}->{'Series Number'},": dims don't fit! $s vs. ",$self->shape->(:4) if any($self->shape->(2:4)-$s);
	#say "Rot: $rot";
	#if ($pe eq 'COL') {
	initdim($self,'c');
	initdim($self,'fid'	,min=>0,inc=>$self->hdr->{ascconv}->{sRXSPEC_alDwellTime_0_},unit=>'us');
	initdim($self,getx($self->hdr),size=>sclr($s(0)),min=>sclr($pos_d(0)),inc=>sclr($inc_d(0)),unit=>'mm');
	initdim($self,gety($self->hdr),size=>sclr($s(1)),min=>sclr($pos_d(1)),inc=>sclr($inc_d(1)),unit=>'mm');
	#} else {
	#initdim($self,getx($self->hdr),size=>$s(0),min=>sclr($pos_d(0)),inc=>sclr($inc_d(0)),unit=>'mm');
	#initdim($self,gety($self->hdr),size=>$s(1),min=>sclr($pos_d(1)),inc=>sclr($inc_d(1)),unit=>'mm');
	#}
	initdim($self,getz($self->hdr),size=>sclr($s(2)),min=>sclr($pos_d(2)),inc=>sclr($inc_d(2)),unit=>'mm',);
	#say "initdim for x,y,z done.";
	#say "after init dim ",(diminfo ($self));
	#say "size $s min $pos_d inc $inc_d rot $rot";
	idx($self,'x',sclr(dimsize($self,'x')/2));
	idx($self,'y',sclr(dimsize($self,'y')/2));
	idx($self,'z',sclr(dimsize($self,'z')/2));
	say "Index x ",idx($self,'x');
	say "Dimsize ",dimsize($self);
	# other dimensions
	for my $n (5..$#{$$opt{MRSSiemens}->{Dimensions}}) { # x,y,z are handled above
		my $dim=$$opt{MRSSiemens}->{Dimensions}->[$n];
		print "Init Dim $dim - $n\n";
		my $str=('(0),' x ($n-2)).','.('(0),' x ($#{$$opt{Dimensions}}-$n));
		#say "$str ";