TPCCLIB
Loading...
Searching...
No Matches
Functions
dcmimage.c File Reference

Extract information from DICOM image. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpcdcm.h"

Go to the source code of this file.

Functions

int dcmImgIsotope (DCMFILE *d, isotope *isot, decaycorrection *dc, const int verbose)
 
int dcmImgPos (DCMFILE *d, double *imgpos, const int verbose)
 
int dcmImgDim (DCMFILE *d, unsigned short int *imgdim, const int verbose)
 
int dcmImgPxlsize (DCMFILE *d, double *pxlsize, const int verbose)
 
int dcmImgOrient (DCMFILE *d, double *iop, const int verbose)
 
int dcmImgXform (double *iop, double *xyzMM, double *imgPos, double *xform, const int verbose)
 
int dcmXformToQuatern (double *xform, double *quatern, double *qoffset, const int verbose)
 

Detailed Description

Extract information from DICOM image.

Definition in file dcmimage.c.

Function Documentation

◆ dcmImgDim()

int dcmImgDim ( DCMFILE * d,
unsigned short int * imgdim,
const int verbose )

Get DICOM image dimensions.

See also
dcmImgOrient, dcmImgPxlsize
Returns
0 if successful, otherwise >0.
Parameters
dPointer to DCMFILE.
imgdimPointer to imgdim array of x, y, z, and frame numbers, to be filled here.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 232 of file dcmimage.c.

239 {
240 if(verbose>0) {printf("%s(dcmfile, imgdim)\n", __func__); fflush(stdout);}
241 if(d==NULL || d->item==NULL || imgdim==NULL) return(1);
242 for(int i=0; i<4; i++) imgdim[i]=0;
243
244 DCMTAG tag; DCMITEM *iptr;
245
246 /* Get image row number (dimy) */
247 dcmTagSet(&tag, 0x0028, 0x0010);
248 iptr=dcmFindTag(d->item, 0, &tag, 0);
249 if(iptr==NULL) {
250 if(verbose>0) {fprintf(stderr, "Error: cannot find matrix row number.\n"); fflush(stderr);}
251 return(11);
252 }
253 imgdim[1]=(unsigned short int)dcmitemGetInt(iptr);
254
255 /* Get image column number (dimx) */
256 dcmTagSet(&tag, 0x0028, 0x0011);
257 iptr=dcmFindTag(d->item, 0, &tag, 0);
258 if(iptr==NULL) {
259 if(verbose>0) {fprintf(stderr, "Error: cannot find matrix column number.\n"); fflush(stderr);}
260 return(21);
261 }
262 imgdim[0]=(unsigned short int)dcmitemGetInt(iptr);
263
264 /* Get image plane/slice number (dimz) */
265 int a, b;
266 if(verbose>1) {printf("trying to read plane number as 'Number of Slices'\n"); fflush(stdout);}
267 dcmTagSet(&tag, 0x0054, 0x0081);
268 iptr=dcmFindTag(d->item, 0, &tag, 0);
269 if(iptr!=NULL) {
270 imgdim[2]=(unsigned short int)dcmitemGetInt(iptr);
271 if(verbose>2) printf("planeNr := %u\n", imgdim[2]);
272 } else {
273 if(verbose>1) {printf("trying to read plane number from 'In Stack Position Number'\n"); fflush(stdout);}
274 dcmTagSet(&tag, 0x0020, 0x9057);
275 if(dcmTagIntRange(d->item, &tag, &a, &b, verbose-100)==0) {
276 imgdim[2]=(unsigned short int)b;
277 if(verbose>2) printf("planeNr := %u\n", imgdim[2]);
278 } else {
279 if(verbose>1) {printf("trying to read plane number from 'Dimension Index Values'\n"); fflush(stdout);}
280 dcmTagSet(&tag, 0x0020, 0x9157);
281 if(dcmTagIntRange(d->item, &tag, &a, &b, verbose-1)==0) {
282 imgdim[2]=(unsigned short int)b;
283 if(verbose>2) printf("planeNr := %u\n", imgdim[2]);
284 } else {
285 if(verbose>0) {fprintf(stderr, "Warning: cannot find plane number.\n"); fflush(stderr);}
286 imgdim[2]=1;
287 }
288 }
289 }
290
291
292 /* Get frame number */
293 if(verbose>1) {printf("trying to get frame number from 'Temporal position index'\n"); fflush(stdout);}
294 dcmTagSet(&tag, 0x0020, 0x9128);
295 if(dcmTagIntRange(d->item, &tag, &a, &b, verbose-100)==0) {
296 imgdim[3]=(unsigned short int)b;
297 if(verbose>2) printf("frameNr := %u\n", imgdim[3]);
298 } else {
299 if(verbose>1) {printf("trying to get frame number from 'Number of Time Slices'\n"); fflush(stdout);}
300 dcmTagSet(&tag, 0x0054, 0x0101);
301 iptr=dcmFindTag(d->item, 0, &tag, 0);
302 if(iptr!=NULL) {
303 imgdim[3]=(unsigned short int)dcmitemGetInt(iptr);
304 if(verbose>2) printf("frameNr := %u\n", imgdim[3]);
305 } else {
306 if(verbose>0) {fprintf(stderr, "Warning: cannot find frame number.\n"); fflush(stderr);}
307 imgdim[3]=1;
308 }
309 }
310
311 return(0);
312}
dcmitemGetInt
long int dcmitemGetInt(DCMITEM *d)
Definition dcmdata.c:295
dcmTagSet
void dcmTagSet(DCMTAG *tag, unsigned short int group, unsigned short int element)
Definition dcmdata.c:483
dcmFindTag
DCMITEM * dcmFindTag(DCMITEM *d, const short int omit, DCMTAG *tag, const int verbose)
Definition dcmdata.c:375
dcmTagIntRange
int dcmTagIntRange(DCMITEM *d, DCMTAG *tag, int *mi, int *ma, const int verbose)
Definition dcmdata.c:631
DCMFILE::item
DCMITEM * item
Definition tpcdcm.h:164

Referenced by imgReadDICOM().

◆ dcmImgIsotope()

int dcmImgIsotope ( DCMFILE * d,
isotope * isot,
decaycorrection * dc,
const int verbose )

Get isotope and status of decaycorrection from DICOM.

See also
isotopeName, lambdaFromIsotope, decayCorrectionFactorFromisotope,
Returns
0 if successful, otherwise >0.
Parameters
dPointer to DCMFILE.
isotPointer to isotope code, to be set here; enter NULL if not needed.
dcPointer to status of decaycorrection.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 24 of file dcmimage.c.

33 {
34 if(verbose>0) {printf("%s(dcmfile, *isotope, *decaycorrection)\n", __func__); fflush(stdout);}
35 if(isot!=NULL) *isot=ISOTOPE_UNKNOWN;
36 if(dc!=NULL) *dc=DECAY_UNKNOWN;
37 if(d==NULL || d->item==NULL) return(1);
38
39
40 DCMTAG tag; DCMITEM *iptr, *jptr;
41
42 /* Get the status of decay correction */
43 if(dc!=NULL) {
44 if(verbose>1) {printf("trying to find 'Decay Correction' tag\n"); fflush(stdout);}
45 dcmTagSet(&tag, 0x0054, 0x1102); iptr=dcmFindTag(d->item, 0, &tag, 0);
46 if(iptr!=NULL) {
47 char *buf=dcmValueString(iptr);
48 if(verbose>12) printf("Decay Correction := '%s'\n", buf);
49 if(strcasestr(buf, "NONE")!=NULL) *dc=DECAY_NOTCORRECTED;
50 else if(strcasestr(buf, "START")!=NULL) *dc=DECAY_CORRECTED_START;
51 else if(strcasestr(buf, "ADMIN")!=NULL) *dc=DECAY_CORRECTED_ADMIN;
52 free(buf);
53 } else {
54 if(verbose>1) {printf(" not found\ntrying to find 'Decay Corrected' tag\n"); fflush(stdout);}
55 dcmTagSet(&tag, 0x0018, 0x9758); iptr=dcmFindTag(d->item, 0, &tag, 0);
56 if(iptr!=NULL) {
57 char *buf=dcmValueString(iptr);
58 if(verbose>12) printf("Decay Correction := '%s'\n", buf);
59 if(strcasestr(buf, "YES")!=NULL) *dc=DECAY_CORRECTED;
60 else if(strcasestr(buf, "NO")!=NULL) *dc=DECAY_NOTCORRECTED;
61 free(buf);
62 } else if(verbose>1) {printf(" not found.\n"); fflush(stdout);}
63 }
64 if(verbose>0 && *dc==DECAY_UNKNOWN) {
65 fprintf(stderr, "Warning: status of decay correction not found.\n");
66 fflush(stdout);
67 }
68 if(verbose>11) printf("decayCorrection := %s\n", decayDescr(*dc));
69 }
70
71 if(isot==NULL) return(0);
72
73 /* Get the radioisotope */
74 if(verbose>1) {printf("trying to find Radiopharmaceutical Information Sequence\n"); fflush(stdout);}
75 dcmTagSet(&tag, 0x0054, 0x0016); iptr=dcmFindTag(d->item, 0, &tag, 0);
76 if(iptr!=NULL) {
77 if(verbose>1) {printf(" found.\n"); fflush(stdout);}
78 iptr=iptr->child_item;
79 } else {
80 if(verbose>1) {printf(" not found.\n"); fflush(stdout);}
81 iptr=d->item;
82 }
83 if(verbose>1) {printf("trying to find Radionuclide Code Sequence\n"); fflush(stdout);}
84 dcmTagSet(&tag, 0x0054, 0x0300); jptr=dcmFindTag(iptr, 0, &tag, 0);
85 if(jptr!=NULL) {
86 if(verbose>1) {printf(" found.\n"); fflush(stdout);}
87 if(verbose>1) {printf("trying to find Code Meaning\n"); fflush(stdout);}
88 dcmTagSet(&tag, 0x0008, 0x0104); jptr=dcmFindDownTag(jptr, 0, &tag, 0);
89 if(jptr!=NULL) {
90 char *buf=dcmValueString(jptr);
91 if(verbose>12) {printf("Code Meaning := '%s'\n", buf); fflush(stdout);}
92 *isot=isotopeIdentify(buf);
93 free(buf);
94 if(verbose>11) {printf("isotope := %s\n", isotopeName(*isot)); fflush(stdout);}
95 } else if(verbose>1) {printf(" not found.\n"); fflush(stdout);}
96 } else {
97 if(verbose>1) {printf(" not found.\n"); fflush(stdout);}
98 }
99 if(*isot==ISOTOPE_UNKNOWN) {
100 if(verbose>1) {printf("trying to find Radionuclide Half Life\n"); fflush(stdout);}
101 dcmTagSet(&tag, 0x0018, 0x1075); jptr=dcmFindTag(iptr, 0, &tag, 0);
102 if(jptr!=NULL) {
103 char *buf=dcmValueString(jptr); strCleanSpaces(buf);
104 if(verbose>12) {printf("Radionuclide Half Life := '%s'\n", buf); fflush(stdout);}
105 double hl=atofVerified(buf);
106 *isot=isotopeIdentifyHalflife(hl/60.0);
107 free(buf);
108 if(verbose>11) {printf("isotope := %s\n", isotopeName(*isot)); fflush(stdout);}
109 } else if(verbose>1) {printf(" not found.\n"); fflush(stdout);}
110 }
111
112 return(0);
113}
dcmFindDownTag
DCMITEM * dcmFindDownTag(DCMITEM *d, const short int omit, DCMTAG *tag, const int verbose)
Definition dcmdata.c:428
dcmValueString
char * dcmValueString(DCMITEM *d)
Definition dcmdata.c:141
decayDescr
char * decayDescr(decaycorrection d)
Definition decay.c:32
atofVerified
double atofVerified(const char *s)
Definition decpoint.c:75
isotopeName
char * isotopeName(int isotope_code)
Definition isotope.c:101
isotopeIdentifyHalflife
int isotopeIdentifyHalflife(double halflife)
Definition isotope.c:121
isotopeIdentify
int isotopeIdentify(const char *isotope)
Definition isotope.c:145
strCleanSpaces
int strCleanSpaces(char *s)
Definition stringext.c:300
strcasestr
char * strcasestr(const char *haystack, const char *needle)
Definition stringext.c:155
DCMITEM::child_item
struct DCMITEM * child_item
Definition tpcdcm.h:143
DECAY_UNKNOWN
@ DECAY_UNKNOWN
Not known; usually assumed that data is corrected.
Definition tpcisotope.h:79
DECAY_CORRECTED_ADMIN
@ DECAY_CORRECTED_ADMIN
Data is corrected for physical decay to radioligand administration time.
Definition tpcisotope.h:83
DECAY_NOTCORRECTED
@ DECAY_NOTCORRECTED
Data is not corrected for physical decay.
Definition tpcisotope.h:80
DECAY_CORRECTED
@ DECAY_CORRECTED
Data is corrected for physical decay.
Definition tpcisotope.h:81
DECAY_CORRECTED_START
@ DECAY_CORRECTED_START
Data is corrected for physical decay to scan start time.
Definition tpcisotope.h:82
ISOTOPE_UNKNOWN
@ ISOTOPE_UNKNOWN
Unknown.
Definition tpcisotope.h:51

Referenced by imgReadDICOM().

◆ dcmImgOrient()

int dcmImgOrient ( DCMFILE * d,
double * iop,
const int verbose )

Get DICOM image orientation (Patient).

See also
dcmImgPos, dcmImgDim, dcmImgPxlsize, dcmImgXform, dcmXformToQuatern
Returns
0 if successful, otherwise >0.
Parameters
dPointer to DCMFILE.
iopPointer to iop array of six values, to be filled here.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 387 of file dcmimage.c.

394 {
395 if(verbose>0) {printf("%s(dcmfile, iop)\n", __func__); fflush(stdout);}
396 if(d==NULL || d->item==NULL || iop==NULL) return(1);
397 for(int i=0; i<6; i++) iop[i]=nan("");
398
399 DCMTAG tag; DCMITEM *iptr;
400
401 if(verbose>1) {printf("reading Image Orientation (Patient)\n"); fflush(stdout);}
402 dcmTagSet(&tag, 0x0020, 0x0037);
403 iptr=dcmFindTag(d->item, 0, &tag, 0);
404 if(iptr==NULL) {
405 if(verbose>0) {fprintf(stderr, "Error: no Image Orientation (Patient)\n"); fflush(stderr);}
406 return(11);
407 }
408 char *buf=dcmValueString(iptr);
409 if(verbose>2) printf(" image_orientation := '%s'\n", buf);
410 int n=sscanf(buf, "%lf\\%lf\\%lf\\%lf\\%lf\\%lf", &iop[0], &iop[1], &iop[2],
411 &iop[3], &iop[4], &iop[5]);
412 free(buf);
413 if(n!=6) {
414 if(verbose>0) {fprintf(stderr, "Error: invalid Image Orientation (Patient)\n"); fflush(stderr);}
415 return(12);
416 }
417
418 return(0);
419}

Referenced by imgReadDICOM().

◆ dcmImgPos()

int dcmImgPos ( DCMFILE * d,
double * imgpos,
const int verbose )

Get DICOM Image Position (Patient).

See also
dcmImgOrient, dcmImgDim, dcmImgPxlsize, dcmImgXform, dcmXformToQuatern
Returns
0 if successful, otherwise >0.
Parameters
dPointer to DCMFILE.
imgposPointer to imgpos array of x, y, and z coordinates of the upper left hand corner (centre of the first voxel) of the image, in mm, to be filled here.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 121 of file dcmimage.c.

129 {
130 if(verbose>0) {printf("%s(dcmfile, imgpos)\n", __func__); fflush(stdout);}
131 if(d==NULL || d->item==NULL || imgpos==NULL) return(1);
132 for(int i=0; i<3; i++) imgpos[i]=nan("");
133
134 DCMTAG tag; DCMITEM *iptr, *jptr;
135
136 if(verbose>1) {printf("trying to find Per Frame Functional Groups Sequence\n"); fflush(stdout);}
137 dcmTagSet(&tag, 0x5200, 0x9230); // Per Frame Functional Groups Sequence
138 iptr=dcmFindTag(d->item, 0, &tag, 0);
139 if(iptr!=NULL) {
140 if(verbose>1) {printf(" found.\n"); fflush(stdout);}
141
142 /* search for Frame Content Sequences until we find one with position number 1 */
143 int posnr=0; dcmTagSet(&tag, 0x0020, 0x9111);
144 jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
145 while(jptr!=NULL) {
146 /* Get In Stack Position Number */
147 DCMITEM *kptr; DCMTAG spostag; dcmTagSet(&spostag, 0x0020, 0x9057);
148 kptr=dcmFindDownTag(jptr->child_item, 0, &spostag, 0);
149 if(kptr!=NULL) {
150 posnr=dcmitemGetInt(kptr); if(verbose>3) printf(" posnr=%d\n", posnr);
151 if(posnr==1) break; // found stack position number 1
152 }
153 jptr=dcmFindDownTag(jptr->next_item, 0, &tag, 0);
154 }
155 /* If stack position 1 was found, then get plane position sequence */
156 if(posnr==1 && jptr!=NULL) {
157 dcmTagSet(&tag, 0x0020, 0x9113);
158 jptr=dcmFindDownTag(jptr->next_item, 0, &tag, 0);
159 if(jptr!=NULL) {
160 /* Get Image Position (Patient) under it */
161 DCMITEM *kptr; DCMTAG ipostag; dcmTagSet(&ipostag, 0x0020, 0x0032);
162 kptr=dcmFindDownTag(jptr->child_item, 0, &ipostag, 0);
163 if(kptr!=NULL) {
164 char *buf=dcmValueString(kptr);
165 if(verbose>4) printf(" image_position := '%s'\n", buf);
166 int n=sscanf(buf, "%lf\\%lf\\%lf", &imgpos[0], &imgpos[1], &imgpos[2]);
167 free(buf);
168 if(n!=3) {
169 if(verbose>0) fprintf(stderr, "Error: invalid Image Position (Patient)\n");
170 return(131);
171 }
172 } else {
173 if(verbose>0) fprintf(stderr, "Error: no Image Position (Patient)\n");
174 return(121);
175 }
176 } else {
177 if(verbose>0) fprintf(stderr, "Error: no Plane Position Sequence\n");
178 return(111);
179 }
180 } else {
181 if(verbose>0) {fprintf(stderr, "Error: stack position 1 not found.\n"); fflush(stderr);}
182 return(101);
183 }
184 return(0);
185 }
186 if(verbose>1) {printf(" not found.\n"); fflush(stdout);}
187
188 if(verbose>1) {printf("find the smallest Instance (image) Number\n"); fflush(stdout);}
189 dcmTagSet(&tag, 0x0020, 0x0013);
190 iptr=dcmFindTag(d->item, 0, &tag, 0);
191 if(iptr==NULL) {
192 if(verbose>0) {fprintf(stderr, "Error: Instance (image) Number not found.\n"); fflush(stderr);}
193 return(201);
194 }
195 int mininr=dcmitemGetInt(iptr); if(verbose>5) printf(" inr=%d\n", mininr);
196 jptr=iptr;
197 while(jptr!=NULL) {
198 jptr=dcmFindDownTag(jptr->next_item, 0, &tag, 0);
199 if(jptr==NULL) break;
200 int inr=dcmitemGetInt(jptr); if(verbose>5) printf(" inr=%d\n", inr);
201 if(inr<mininr) {mininr=inr; iptr=jptr;}
202 }
203 if(verbose>2) {printf(" smallest Instance (image) Number: %d\n", mininr); fflush(stdout);}
204
205 /* Get Image Position (Patient) after the smallest Instance Number */
206 dcmTagSet(&tag, 0x0020, 0x0032);
207 jptr=dcmFindDownTag(iptr, 0, &tag, 0);
208 if(jptr==NULL) {
209 if(verbose>0) fprintf(stderr, "Error: no Image Position (Patient)\n");
210 return(211);
211 } else {
212 char *buf=dcmValueString(jptr);
213 if(verbose>2) printf(" image_position := '%s'\n", buf);
214 int n=sscanf(buf, "%lf\\%lf\\%lf", &imgpos[0], &imgpos[1], &imgpos[2]);
215 free(buf);
216 if(n!=3) {
217 if(verbose>0) fprintf(stderr, "Error: invalid Image Position (Patient)\n");
218 return(212);
219 }
220 }
221
222 return(0);
223}
DCMITEM::next_item
struct DCMITEM * next_item
Definition tpcdcm.h:147

Referenced by imgReadDICOM().

◆ dcmImgPxlsize()

int dcmImgPxlsize ( DCMFILE * d,
double * pxlsize,
const int verbose )

Get DICOM Image voxel size.

See also
dcmImgDim, dcmImgOrient
Returns
0 if successful, otherwise >0.
Parameters
dPointer to DCMFILE.
pxlsizePointer to pxlsize array of x, y, and z sizes of image pixel, or actually, mm distance between neighbour voxel centres in x, y, and z directions.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 320 of file dcmimage.c.

329 {
330 if(verbose>0) {printf("%s(dcmfile, pxlsize)\n", __func__); fflush(stdout);}
331 if(d==NULL || d->item==NULL || pxlsize==NULL) return(1);
332 for(int i=0; i<3; i++) pxlsize[i]=nan("");
333
334 DCMTAG tag; DCMITEM *iptr, *jptr, *kptr;
335
336 if(verbose>1) {printf("trying to find Pixel Measures Sequence\n"); fflush(stdout);}
337 dcmTagSet(&tag, 0x0028, 0x9110);
338 iptr=dcmFindTag(d->item, 0, &tag, 0);
339 if(iptr==NULL) {
340 kptr=d->item; // search from root
341 if(verbose>3) {fprintf(stdout, "Note: missing Pixel Measures Sequence.\n"); fflush(stdout);}
342 } else {
343 kptr=iptr->child_item; // search under this sequence
344 }
345 dcmTagSet(&tag, 0x0018, 0x0050);
346 jptr=dcmFindDownTag(kptr, 0, &tag, 0);
347 if(jptr==NULL) {
348 if(verbose>0) {fprintf(stderr, "Error: missing Slice Thickness.\n"); fflush(stderr);}
349 return(11);
350 }
351 pxlsize[2]=dcmitemGetReal(jptr);
352
353 if(verbose>1) {printf("trying to find Spacing Between Slices\n"); fflush(stdout);}
354 dcmTagSet(&tag, 0x0018, 0x0088);
355 jptr=dcmFindDownTag(kptr, 0, &tag, 0);
356 if(jptr==NULL) {
357 if(verbose>1) {fprintf(stderr, "Warning: missing Spacing Between Slices.\n"); fflush(stderr);}
358 } else {
359 double v=dcmitemGetReal(jptr);
360 if(v>pxlsize[2]) pxlsize[2]=v;
361 }
362
363 if(verbose>1) {printf("trying to find Pixel Spacing/Size\n"); fflush(stdout);}
364 dcmTagSet(&tag, 0x0028, 0x0030);
365 jptr=dcmFindDownTag(kptr, 0, &tag, 0);
366 if(jptr==NULL) {
367 if(verbose>0) {fprintf(stderr, "Error: missing Pixel Spacing/Size.\n"); fflush(stderr);}
368 return(21);
369 }
370 char *buf=dcmValueString(jptr);
371 int n=sscanf(buf, "%lf\\%lf", &pxlsize[1], &pxlsize[0]);
372 free(buf);
373 if(n!=2) {
374 if(verbose>0) {fprintf(stderr, "Error: invalid Pixel Spacing/Size\n"); fflush(stderr);}
375 return(22);
376 }
377
378 return(0);
379}
dcmitemGetReal
double dcmitemGetReal(DCMITEM *d)
Definition dcmdata.c:331

Referenced by imgReadDICOM().

◆ dcmImgXform()

int dcmImgXform ( double * iop,
double * xyzMM,
double * imgPos,
double * xform,
const int verbose )

Calculate xform matrix from DICOM header information.

Based on dcm2niix xform_mat(), which is based on dicm2nii. Last vector xform[3][0..4] is not yet computed.

See also
dcmImgOrient, dcmImgPxlsize, dcmImgPos, dcmImgDim
Returns
0 if successful, otherwise >0.
Parameters
iopArray of 6 Image Orientation (Patient) parameters.
xyzMMArray of 3 voxel sizes in mm (x, y, z).
imgPosArray of 3 Image Position (Patient) parameters.
xformPointer to xform 4x4 matrix stored as an array of 16 doubles, to be filled here.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 431 of file dcmimage.c.

442 {
443 if(verbose>0) {printf("%s(iop[], xyzMM[], imgPos[], xform)\n", __func__); fflush(stdout);}
444 if(iop==NULL || xyzMM==NULL || imgPos==NULL || xform==NULL) return(1);
445 for(int i=0; i<16; i++) xform[i]=nan("");
446
447
448 double readV[3], phaseV[3];
449 readV[0]=iop[0]; readV[1]=iop[1]; readV[2]=iop[2];
450 phaseV[0]=iop[3]; phaseV[1]=iop[4]; phaseV[2]=iop[5];
451 if(verbose>3) {
452 printf("readV := %g", readV[0]); for(int i=1; i<3; i++) printf(", %g", readV[i]);
453 printf("\n");
454 printf("phaseV := %g", phaseV[0]); for(int i=1; i<3; i++) printf(", %g", phaseV[i]);
455 printf("\n");
456 }
457
458 /* Cross-product of vectors */
459 double sliceV[3];
460 sliceV[0]= readV[1]*phaseV[2] - phaseV[1]*readV[2];
461 sliceV[1]=-readV[0]*phaseV[2] + phaseV[0]*readV[2];
462 sliceV[2]= readV[0]*phaseV[1] - phaseV[0]*readV[1];
463 if(verbose>3) {
464 printf("sliceV := %g", sliceV[0]); for(int i=1; i<3; i++) printf(", %g", sliceV[i]);
465 printf("\n");
466 }
467
468 /* Copy vectors into a 3x3 matrix */
469 double rM[3][3];
470 for(int i=0; i<3; i++) rM[0][i]=readV[i];
471 for(int i=0; i<3; i++) rM[1][i]=phaseV[i];
472 for(int i=0; i<3; i++) rM[2][i]=sliceV[i];
473 /* Transpose the matrix */
474 {
475 double tM[3][3];
476 for(int i=0; i<3; i++) for(int j=0; j<3; j++) tM[i][j]=rM[i][j];
477 for(int i=0; i<3; i++) for(int j=0; j<3; j++) rM[i][j]=tM[j][i];
478 }
479 if(verbose>3) {
480 printf("rM :=\t%g", rM[0][0]);
481 for(int i=1; i<3; i++) printf("\t%g", rM[0][i]);
482 printf("\n");
483 for(int i=0; i<3; i++) printf("\t%g", rM[1][i]);
484 printf("\n");
485 for(int i=0; i<3; i++) printf("\t%g", rM[2][i]);
486 printf("\n");
487 }
488 /* Max (left-over) indices */
489 unsigned short int xyz[3];
490 xyz[0]=0;
491 if(fabs(rM[1][0])>fabs(rM[0][0]) && fabs(rM[1][0])>fabs(rM[2][0])) xyz[0]=1;
492 else if(fabs(rM[2][0])>fabs(rM[0][0]) && fabs(rM[1][0])>fabs(rM[1][0])) xyz[0]=2;
493 if(xyz[0]==0) {if(fabs(rM[2][1])>fabs(rM[1][1])) xyz[1]=2; else xyz[1]=1;}
494 else if(xyz[0]==1) {if(fabs(rM[2][1])>fabs(rM[0][1])) xyz[1]=2; else xyz[1]=0;}
495 else {if(fabs(rM[1][1])>fabs(rM[0][1])) xyz[1]=1; else xyz[1]=0;}
496 xyz[2]=3-xyz[0]-xyz[1];
497 if(verbose>3) {
498 printf("xyz := %u", xyz[0]); for(int i=1; i<3; i++) printf(", %u", xyz[i]);
499 printf("\n(0, 1, 2 for Sag/Cor/Tra slice)\n");
500 }
501 double cosSL=rM[xyz[2]][2];
502 if(verbose>3) {printf("cosSL := %g\n", cosSL); fflush(stdout);}
503
504 /* Multiply matrix with voxel sizes */
505 if(verbose>3) {
506 printf("xyzMM := %g", xyzMM[0]); for(int i=1; i<3; i++) printf(", %g", xyzMM[i]);
507 printf("\n");
508 }
509 {
510 double tM[3][3], sM[3][3];
511 for(int i=0; i<3; i++) for(int j=0; j<3; j++) tM[i][j]=rM[i][j];
512 for(int i=0; i<3; i++) for(int j=0; j<3; j++) sM[i][j]=0.0;
513 sM[0][0]=xyzMM[0]; sM[1][1]=xyzMM[1]; sM[2][2]=xyzMM[2];
514 for(int i=0; i<3; i++)
515 for(int j=0; j<3; j++)
516 rM[i][j]= tM[i][0]*sM[0][j] + tM[i][1]*sM[1][j] + tM[i][2]*sM[2][j];
517 }
518
519 /* Fill most part of xform */
520 for(int i=0; i<3; i++)
521 for(int j=0; j<3; j++)
522 xform[j+4*i]=rM[i][j];
523 for(int i=0; i<3; i++)
524 xform[3+4*i]=imgPos[i];
525
526 /* Slice normalization vector not yet done */
527
528 return(0);
529}

Referenced by imgReadDICOM().

◆ dcmXformToQuatern()

int dcmXformToQuatern ( double * xform,
double * quatern,
double * qoffset,
const int verbose )

Calculate quatern parameters for NIfTI from xform matrix.

Based on dcm2niix, but results are compatible with TPC ImageConverter.

See also
dcmImgXform, dcmImgOrient, dcmImgPxlsize, dcmImgPos, dcmImgDim
Returns
0 if successful, otherwise >0.
Parameters
xformPointer to xform 4x4 matrix stored as an array of 16 doubles; not modified.
quaternPointer to quatern parameter b, c, and d array (size of 3); filled here; enter NULL, if not needed.
qoffsetPointer to qoffset parameter x, y, and z array (size of 3); filled here; enter NULL, if not needed.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 540 of file dcmimage.c.

551 {
552 if(verbose>0) {printf("%s(xform[16], quatern[3], qoffset[3])\n", __func__); fflush(stdout);}
553 if(xform==NULL) return(1);
554 if(quatern!=NULL) for(int i=0; i<3; i++) quatern[i]=nan("");
555 if(qoffset!=NULL) for(int i=0; i<3; i++) qoffset[i]=nan("");
556
557 /* Copy array to matrix */
558 double xformM[4][4];
559 for(int i=0; i<4; i++)
560 for(int j=0; j<4; j++)
561 xformM[i][j]=xform[j+4*i];
562 /* DICOM LPS to NIfTI RAS */
563 for(int i=0; i<2; i++)
564 for(int j=0; j<4; j++)
565 xformM[i][j]=-xformM[i][j];
566
567 /* Set quatern offsets */
568 if(qoffset!=NULL) for(int i=0; i<3; i++) qoffset[i]=xformM[i][3];
569
570 /* Ready, if quatern parameters not needed */
571 if(quatern==NULL) return(0);
572
573 // compute lengths
574 double mt[3][3];
575 for(int i=0; i<3; i++)
576 for(int j=0; j<3; j++)
577 mt[i][j]=xformM[i][j];
578 double xd, yd, zd;
579 xd=sqrt(mt[0][0]*mt[0][0] + mt[1][0]*mt[1][0] + mt[2][0]*mt[2][0]);
580 yd=sqrt(mt[0][1]*mt[0][1] + mt[1][1]*mt[1][1] + mt[2][1]*mt[2][1]);
581 zd=sqrt(mt[0][2]*mt[0][2] + mt[1][2]*mt[1][2] + mt[2][2]*mt[2][2]);
582 if(!(xd>0.0)) {mt[0][0]=1.0; mt[1][0]=mt[2][0]=0.0; xd=1.0;}
583 if(!(yd>0.0)) {mt[1][1]=1.0; mt[0][1]=mt[2][1]=0.0; yd=1.0;}
584 if(!(zd>0.0)) {mt[2][2]=1.0; mt[0][2]=mt[1][2]=0.0; zd=1.0;}
585 // normalize
586 for(int i=0; i<3; i++) mt[i][0]/=xd;
587 for(int i=0; i<3; i++) mt[i][1]/=yd;
588 for(int i=0; i<3; i++) mt[i][2]/=zd;
589 // check orthogonality (to be added)
590
591 // determinant
592 double dt;
593 dt= mt[0][0]*mt[1][1]*mt[2][2] - mt[0][0]*mt[2][1]*mt[1][2] - mt[1][0]*mt[0][1]*mt[2][2]
594 +mt[1][0]*mt[2][1]*mt[0][2] + mt[2][0]*mt[0][1]*mt[1][2] - mt[2][0]*mt[1][1]*mt[0][2];
595 if(verbose>3) printf("determinant=%g\n", dt);
596 if(dt<0.0) for(int i=0; i<3; i++) mt[i][2]=-mt[i][2];
597
598 // Calculate quatern parameters
599 double a, b, c, d;
600 a=1.0+mt[0][0]+mt[1][1]+mt[2][2];
601 if(a>0.5) {
602 a=0.5*sqrt(a);
603 b=0.25*(mt[2][1]-mt[1][2])/a;
604 c=0.25*(mt[0][2]-mt[2][0])/a;
605 d=0.25*(mt[1][0]-mt[0][1])/a;
606 } else {
607 xd=1.0+mt[0][0]-(mt[1][1]+mt[2][2]);
608 yd=1.0+mt[1][1]-(mt[0][0]+mt[2][2]);
609 zd=1.0+mt[2][2]-(mt[0][0]+mt[1][1]);
610 if(xd>1.0) {
611 b=0.50*sqrt(xd);
612 c=0.25*(mt[0][1]+mt[1][0])/b;
613 d=0.25*(mt[0][2]+mt[2][0])/b;
614 a=0.25*(mt[2][1]-mt[1][2])/b;
615 } else if(yd>1.0) {
616 c=0.50*sqrt(yd);
617 b=0.25*(mt[0][1]+mt[1][0])/c;
618 d=0.25*(mt[1][2]+mt[2][0])/c;
619 a=0.25*(mt[0][2]-mt[2][0])/c;
620 } else {
621 d=0.50*sqrt(zd);
622 b=0.25*(mt[0][2]+mt[2][0])/d;
623 c=0.25*(mt[1][2]+mt[2][1])/d;
624 a=0.25*(mt[1][0]-mt[0][1])/d;
625 }
626 }
627 if(a<0.0) {b=-b; c=-c; d=-d;}
628 quatern[0]=b;
629 quatern[1]=c;
630 quatern[2]=d;
631
632 return(0);
633}

Referenced by imgReadDICOM().

Generated on Fri Sep 19 2025 13:32:19 for TPCCLIB by doxygen 1.12.0