I have the following equation from the Dicom standard to calculate the voxel position in space from a point in an image plane where
Pxyz The coordinates of the voxel (i,j) in the frame's image plane in units of mm.
Sxyz The three values of Image Position (Patient) (0020,0032). It is the location in mm from the origin of the RCS.
Xxyz The values from the row (X) direction cosine of Image Orientation (Patient) (0020,0037).
Yxyz The values from the column (Y) direction cosine of Image Orientation (Patient) (0020,0037).
i Column index to the image plane. The first column is index zero.
Δi Column pixel resolution of Pixel Spacing (0028,0030) in units of mm.
j Row index to the image plane. The first row index is zero.
Δj Row pixel resolution of Pixel Spacing (0028,0030) in units of mm.
What I want to do is solve the whole thing for the image plane coordinate (i,j) parts. I think what I need to do is multiply by the inverse of the XxΔiYxΔj0Sx (etc) part.
1) Is that correct? Answer from comment:
matrix is not invertible because its third column is a zero column
1.A) Maybe this is stupid, but it looks to me like that third column is unnecessary. Isn't this 3x3 equivalent?
2) How do I get the inverse of that 4x4 3x3 matrix?
1) Yes, the $3 \times 3$ system is equivalent to the previous one (I should have noted it...)
2) Let $M$ denote this $3 \times 3$ matrix.
Inverting the equation means that
$U:=\begin{pmatrix}i\\j\\1\end{pmatrix}$ should be equal to $V:=M^{-1} \begin{pmatrix}P_x\\P_y\\P_z\end{pmatrix}.$
The apparently big problem is that the third coordinate of $V$ has no reason to be equal to $1$.
But in fact, as this $1$ coordinate has to be considered in the context of homogeneous coordinates (have you been introduced to them ?), the real relationship betwen $U$ and $V$ is that their coordinates are not equal but are proportional.
Therefore, take
$$i=V(1)/V(3), \ \ j=V(2)/V(3).$$
Remark: I advise you to do the computation of $M^{-1}$ with the usual formulas:
if $M=\begin{pmatrix}a&d&g\\b&e&h\\c&f&i\end{pmatrix}$ (I take the simplest notations ; beware that the $i$ here is not the coordinate $i$ seen before), then:
$$M^{-1}=\dfrac{1}{\det(M)}\begin{pmatrix} ei - fh& fg - di& dh - eg\\ ch - bi& ai - cg& bg - ah\\ bf - ce& cd - af& ae - bd \end{pmatrix},$$
where $\det(M)=aei + bfg + cdh - afh - bdi -ceg.$
with a warning if this determinant is too close from zero.