view release on metacpan or  search on metacpan

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

	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!";
	}
	use PDL::NiceSlice;
	#my $v=zeroes(3);
	# This needs to be converted to dicom fields for navs ...
	# Center of FOV
	my $dims=$self->hdr->{dicom}->{"Image Orientation (Patient)"}->shape;
	my @li=( 3,2,list ( $dims->slice("1:")));
	my $sh= $self->hdr->{dicom}->{"Image Orientation (Patient)"}->reshape(@li);
	my $w=  $self->hdr->{dicom}->{"Pixel Spacing"}->dummy(0,1)*$sh;
	my $vv=($self->hdr->{dicom}->{Columns}*$w->slice(",0")+$self->hdr->{dicom}->{Rows}*$w->slice(",1"))->clump(0,1);
	my  $v=$vv/2+$self->hdr->{dicom}->{"Image Position (Patient)"}; #->transpose;

	#$v(0,).=hpar($self,'ascconv','sSliceArray_asSlice_0__sPosition_dSag') ||0; #x
	#$v(1,).=hpar($self,'ascconv','sSliceArray_asSlice_0__sPosition_dCor') ||0; #y
	#$v(2,).=hpar($self,'ascconv','sSliceArray_asSlice_0__sPosition_dTra') ||0; #z
	#say $v;
	#say "hpar: pos ",hpar($self,'dicom','Image Position (Patient)'),
		#$self->hdr->{dicom}->{'Image Position (Patient)'};
	#say "init_dims: ",hpar($self,'dicom','Rows');
	my $pos_d;
	$pos_d=(hpar($self,'dicom','Image Position (Patient)'))->flat->(:2); #->uniq;
	#my $ir=hpar($self,'ascconv','sSliceArray_asSlice_0__dInPlaneRot') ||0; #radiant
	my $or=pdl(hpar($self,'dicom','Image Orientation (Patient)'))->flat->(:5)
		->reshape(3,2)->transpose; #
	my $pe_dir=hpar($self,'dicom','In-plane Phase Encoding Direction');
# Scaling
# Rotation
	my $srot=zeroes(3,3);
	$srot(:1,).=$or;
	$srot(2,;-).=crossp($or(1,;-),$or(0,;-));
	$srot(2,;-).=pdl[0,0,1] unless any ($srot(2,;-));
	#say "spatial rotation $srot";
	#$pos_d=$pos_d(,0;-);
# 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
	my $pe=$self->hdr->{dicom}->{"In-plane Phase Encoding Direction"} 
	||$self->hdr->{dicom}->{"InPlanePhaseEncodingDirection"}; 
	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';
		}
	}
	my $s=zeroes(3); # matrix size 
	#my $fov=zeroes(3); # FOV
	#$fov(0).=$self->hdr->{ascconv}->{sSliceArray_asSlice_0__dReadoutFOV};
	#$fov(1).=$self->hdr->{ascconv}->{sSliceArray_asSlice_0__dPhaseFOV};
	$s(0).=$self->dim(0);
	$s(1).=$self->dim(1);
	#if ($pe =~ 'COL') {
	#	$s(0).=$self->hdr->{csa}->{Columns};
	#	$s(1).=$self->hdr->{csa}->{Rows};
		#say "COL! $s";
	#} else {
	#	$s(1).=$self->hdr->{dicom}->{Width}||$self->hdr->{dicom}->{Columns};
	#	$s(0).=$self->hdr->{dicom}->{Height}||$self->hdr->{dicom}->{Rows};
	#$fov(1).=$self->hdr->{ascconv}->{sSliceArray_asSlice_0__dReadoutFOV};
	#$fov(0).=$self->hdr->{ascconv}->{sSliceArray_asSlice_0__dPhaseFOV};
	#}
	#say "PE $pe $s $fov " ;
	$s(2).=1;
	# Slice groups!
	
	
	#if ($$opt{split}) {
		#my ($sg,$size)=map_slicegroup($self);
		#$s(2).=$self->hdr->{dicom}->{"Image Orientation (Patient)"}->dim(1);
	#} else {
		#$s(2)*=$self->hdr->{ascconv}->{'sSliceArray_lSize'};
	$s(2).=$self->hdr->{dicom}->{'Image Orientation (Patient)'}->dim(2);
	#}	
	$self->hdr->{'3d'}=1 if (($self->hdr->{dicom}->{MRAcquisitionType}||
		hpar($self,'dicom','MR Acquisition Type')) eq '3D'); # 3D
	if ($self->hdr->{'3d'}) {
		#$s(2).=$self->hdr->{ascconv}->{'sKSpace_lImagesPerSlab'} ;
		#$fov(2).=$self->hdr->{ascconv}->{sSliceArray_asSlice_0__dThickness};
	} else {
		#$fov(2).=$self->hdr->{dicom}->{"Spacing Between Slices"}*$s(2)
		#	||$self->hdr->{ascconv}->{'sSliceArray_asSlice_0__dThickness'};
	}
	#$s(2).=1 if ($s(2)<1);
	#say "FOV $fov matrix $s";
	my $rot=identity($self->ndims);
	my $inc_d=zeroes(3);