function [x,y,z,M]=GetVoxelPositions_OUH(Image)
%This code calculates coordinates in the image cube coordinate system
%which is not nesesarily aligned with the XYZ coordinate system fixed in
%the patient. In the Dicom protocol there is no standard for absolute
%coordinates in the datacube system only for the patient fixed system.
%The follwoing code calculates the coordinates such that the center voxel
%gets the same coordinates as in the patient fixed coordinate system and
%then uses pixel widht and hight and slice distances to calculate the
%other coordinates.

%The output of this funtion is thus not the Dicom coordinates fixed in
%the patient but fixed in the datacube. For data aliged with the patient
%coordinate system (e.g. as typpical for CT data) the two system will be
%identical.
%The returned x y and z are thus the coordiantes that in the Cerr image
%format can be used for entry 3, 4 and 5 

%Relation between imageplane (the system of the imagecube) and the Dicom
%patient coordinates which is fixed in the patient is given in the Dicom protocol
%section C.7.6.2.1.1 "Image Position And Image Orientation"

%The Dicom patient coordinate are defines as follows:
%X axis: From the patients right to left
%Y axis: From front to back of patient
%Z axis: from feet to head of patient

%Let [P_X;P_Y;P_Z] be the patien coordinates and define then the vector
%P=[P_X;P_Y;P_Z; 1]
%Let a translation/rotation matrix be given as
%M=[X_x*Di, Y_x*Dj, 0, Sx; X_y*Di, Y_y*Dj, 0, Sy; X_z*Di, Y_z*Dj, 0, Sz; 0, 0, 0, 1]
%Where:
%X_xyz (e.g. X_x or X_y or X_z) the values of the row (X) direction cosines
%of the Image Orientation (Patient) (0020,0037) attribute
%Y_xyz (e.g. Y_x or Y_y or Y_z) the values of the collumn (Y) direction cosines
%of the Image Orientation (Patient) (0020,0037) attribute
%Sxyz (Sx or Sy or Sz) is the three values of the Image Position
%(Patient) (0020,0032) attributes. It is the location in mm from the
%origin of the RCS (refernce coordinate system which the Dicom coordinate
%system fixed in the patient)
%Di Collumn pixel resolution of the Pixel spacing (0028,0030) attribute
%in units of mm.
%Dj row pixel resolution of the pixel spacing (0028,0030) attribute in
%units of mm

%The coordinates can then be linked by the eqation
%P=M*[i; j; 0; 1]
%i collumn index of the image plane. The first collumn is index zero
%j row index of the image plane. The first row indes is zero

%8th of July 2016 /CaB

%Have added at matrix M which links the returned x y z coordinates with the
%Dicom coorinates. The Dicom coordinates can be calculated by
% DicomCoor = M*(x y z 1)'

%18th of August 2016 /CaB

%Small correction in order to support imagesets of just one image (e.g.
%Overview image)
%31-07-2017 /CaB

%The update was made to support other patient positions than Head first
%supine. The main reason for the need for an update was that the voxel
%coordinates are not necessarily increased as a function of the voxel index
%for other patient positions. This required a bit of rewrite of the code,
%which somewhat slimmed the code.
%13-03-2023 /CaB

%A check was added to ensure that the distance between images is constant.
%If not, a warning will be shown, and the returned values will be NaN
%18-03-2024 /CaB

%Bug related to the offset of the coordinates for tiltel images corrected
%19-01-2025 /CaB

%Start with some check of data
DataProblem=false;
%Dose data do have slightly different datarepresentation.
%Correction such that they can be handled by the same code 
if size(Image{1},2)>1
  if strcmp(Image{1}{2}.Modality,'RTDOSE')
    temp{1}{1}=Image{1}{2};
    Image(1)=temp(1);
  end
end
  
if isfield(Image{1}{1}, 'PixelSpacing')
  temp=cellfun(@(x) x.PixelSpacing(1),Image{1}(:));
  if (max(temp)-min(temp))>0.01*temp(1)
    DataProblem=true;
    warning('PixelSpacing(1) is not constant in data');
  end
  temp=cellfun(@(x) x.PixelSpacing(2),Image{1}(:));
  if (max(temp)-min(temp))>0.01*temp(1)
    DataProblem=true;
    warning('PixelSpacing(2) is not constant in data');
  end
else
  DataProblem=true;
  warning('PixelSpacing is not defined in data');
end
if isfield(Image{1}{1}, 'ImagePositionPatient')
  temp=sort(cellfun(@(x) x.ImagePositionPatient(3),Image{1, 1}));
  %Avoid cheking dose files which does only have one temp value
  if length(temp)>1
  temp=temp(2:end)-temp(1:end-1);
    if (max(temp)-min(temp))>0.01*temp(1)
      DataProblem=true;
      warning('Image position between slices is not constant in the data');
    end
  end
else
  DataProblem=true;
  warning('ImagePositionPatient is not defined in data');
end


%ImageOrientationPatient(1:3)
if isfield(Image{1}{1}, 'ImageOrientationPatient')
  temp=zeros(size(Image{2},3),6);
  for k=1:6
    temp(:,k)=cellfun(@(x) x.ImageOrientationPatient(k),Image{1}(:));
  end
  if max(max(temp,[],1)-min(temp,[],1))>1e-4
    DataProblem=true;
    warning('ImageOrientationPatient is not constant in data');
  end
else
  DataProblem=true;
  warning('ImageOrientationPatient is not defined in data');
end
if isfield(Image{1}{1}, 'Rows') && isfield(Image{1}{1}, 'Columns')
  temp=cellfun(@(x) x.Rows,Image{1}(:));
  if (max(temp)-min(temp))>0.5
    DataProblem=true;
    warning('Number of rows not constant in data');
  end
  temp=cellfun(@(x) x.Columns,Image{1}(:));
  if (max(temp)-min(temp))>0.5
    DataProblem=true;
    warning('Number of collumns not constant in data');
  end
else
  DataProblem=true;
  warning('Rows and Columns are not defined in data');
end

%If okay start caluclation of coordinates
if ~DataProblem
  %Calculate the center index which will get coordinates similar to the
  %Dicom coordinates of the patient fixed coordinates
  zeroindex=round(size(Image{2})/2);
  %Check whether image is just one slice since then we need explicitly to
  %add a 1 for the third dimension
  %31-07-2017 CaB
  if length(zeroindex)==2
      zeroindex=[zeroindex 1];
  end
  %Remember that the first coordinate of matlab is the second coordinate in
  %the dicom standard (and second matlab first in Dicom). Also remember
  %that the index in matlab start at one but in Dicom at zero thus one unit
  %needs to be subtracted
  P=M_matrix(Image,zeroindex(3))*[zeroindex(2)-1;zeroindex(1)-1;0;1];
  P=P(1:3);
  signdirection=1;
  if dot(Image{1}{1}.ImageOrientationPatient(1:3),[1;0;0])<0
    signdirection=-1;
  end
  x=Image{1}{1}.PixelSpacing(2)*signdirection*((1:size(Image{2},2))-zeroindex(2))+P(1);
  signdirection=1;
  if dot(Image{1}{1}.ImageOrientationPatient(4:6),[0;1;0])<0
    signdirection=-1;
  end
  y=Image{1}{1}.PixelSpacing(1)*signdirection*((1:size(Image{2},1))-zeroindex(1))+P(2);
  %Calc z direction of the datacube (cross product of the image vectors in
  %the image plane)
  %zdir=cross(Image{1}{1}.ImageOrientationPatient(1:3),Image{1}{1}.ImageOrientationPatient(4:6));
  %z=zeros(1,size(Image{2},3));
  PixelVector_z=[0;0;1;0];
  if size(Image{2},3)>1
    PixelVector_z=(M_matrix(Image,2)*[1;1;0;1]-M_matrix(Image,1)*[1;1;0;1]);
  end
  PixelSpacing_z=sqrt(dot(PixelVector_z,PixelVector_z));
  signdirection=1;
  if PixelVector_z(3)<0
    signdirection=-1;
  end
  z=PixelSpacing_z*signdirection*((1:size(Image{2},3))-zeroindex(3))+P(3);
  %for k=1:size(Image{2},3)
  %  temp=M_matrix(Image,k)*[zeroindex(2)-1;zeroindex(1)-1;0;1];
  %  z(k)=dot(temp(1:3)-P,zdir)+P(3);
  %end
  M=GetMatrixRotationDicom(Image,x,y,z);
else
  x=NaN(1,size(Image{2},1));
  y=NaN(1,size(Image{2},2));
  z=NaN(1,size(Image{2},3));
  M=NaN(4,4);
end

end
function I=M_matrix(Image,z)
  if strcmpi(Image{1}{1}.Modality,'RTDOSE')
    %For RTDOSE for which only one dicom header exist we need to use the
    %existing header
    X_xyz=Image{1}{1}.ImageOrientationPatient(1:3);
    Y_xyz=Image{1}{1}.ImageOrientationPatient(4:6);
  else
    X_xyz=Image{1}{z}.ImageOrientationPatient(1:3);
    Y_xyz=Image{1}{z}.ImageOrientationPatient(4:6);
  end
  %In the protocol it is stated that the pixel spacing (0028,0030)is
  %defined as:
  %Physical distance in the patient between the center of each pixel,
  %specified by a numeric pair - adjacent row spacing (delimiter) adjacent
  %column spacing in mm.
  %Thus the first value is related dicom Y direction (not X) and the second
  %to the X direction. Di is thus Di=Image{1}{z}.PixelSpacing(2) although
  %it seems strange.
  if strcmpi(Image{1}{1}.Modality,'RTDOSE')
    %For RTDOSE for which only one dicom header exist we need to use the
    %existing header and then offset using the GridOffsetVector
    Di=Image{1}{1}.PixelSpacing(2);
    Dj=Image{1}{1}.PixelSpacing(1);
    S_xyz=Image{1}{1}.ImagePositionPatient;
    S_xyz(3)=S_xyz(3)+Image{1}{1}.GridFrameOffsetVector(z);
   
  else
    Di=Image{1}{z}.PixelSpacing(2);
    Dj=Image{1}{z}.PixelSpacing(1);
    S_xyz=Image{1}{z}.ImagePositionPatient;
  end
  I=[[X_xyz*Di;0], [Y_xyz*Dj;0], zeros(4,1), [S_xyz;1]];
end
function I=GetMatrixRotationDicom(Image,x,y,z)
  %MeanIndex=round(size(Image{2})/2);
  %Check whether image is just one slice since then we need explicitly to
  %add a 1 for the third dimension
  %31-07-2017 CaB
  %if length(MeanIndex)==2
  %    MeanIndex=[MeanIndex 1];
  %end
  %DicomMeanPos=[x(MeanIndex(2));y(MeanIndex(1));z(MeanIndex(3))];
  %I=M_matrix(Image,MeanIndex(3));
  %I(1:3,1)=I(1:3,1)/sqrt(dot(I(1:3,1),I(1:3,1)));
  %I(1:3,2)=I(1:3,2)/sqrt(dot(I(1:3,2),I(1:3,2)));
  %I(1:3,3)=cross(I(1:3,1),I(1:3,2));
  %I(1:4,4)=I(1:4,4)+[DicomMeanPos;1]-I*[DicomMeanPos;1];
  
  M=zeros(4,4);
  M(:,1)=(M_matrix(Image,1)*([1;0;0;1])-M_matrix(Image,1)*([0;0;0;1]))/(x(2)-x(1));
  M(:,2)=(M_matrix(Image,1)*([0;1;0;1])-M_matrix(Image,1)*([0;0;0;1]))/(y(2)-y(1));
  M(:,3)=(M_matrix(Image,2)*([0;0;0;1])-M_matrix(Image,1)*([0;0;0;1]))/(z(2)-z(1));
  %M(:,4)=M_matrix(Image,1)*([0;0;0;1])-[x(1);y(1);z(1);1];
  M(:,4)=M_matrix(Image,1)*([0;0;0;1])-M*[x(1);y(1);z(1);1];
  %M(4,4)=1;
  I=M;
end  