Doxygen Source Code Documentation
3dDTtoDWI.c File Reference
#include "thd_shear3d.h"
#include "matrix.h"
#include "matrix.c"
#include "afni.h"
Go to the source code of this file.
Defines | |
#define | TINYNUMBER 1E-10 |
#define | SMALLNUMBER 1E-4 |
Functions | |
void | DTtoDWI_tsfunc (double tzero, double tdelta, int npts, float ts[], double ts_mean, double ts_slope, void *ud, int nbriks, float *val) |
void | Computebmatrix (MRI_IMAGE *grad1Dptr) |
void | InitGlobals (int npts) |
void | FreeGlobals (void) |
int | main (int argc, char *argv[]) |
Variables | |
char | prefix [THD_MAX_PREFIX] = "DT" |
int | datum = MRI_float |
matrix | Dmatrix |
vector | Dvector |
double * | bmatrix |
float * | I0_ptr |
int | automask = 0 |
Define Documentation
|
Definition at line 15 of file 3dDTtoDWI.c. |
|
Definition at line 14 of file 3dDTtoDWI.c. |
Function Documentation
|
compute the diffusion weighting matrix bmatrix for q number of gradients Definition at line 422 of file 3dDTtoDWI.c. References bmatrix, ENTRY, i, MRI_FLOAT_PTR, and MRI_IMAGE::nx.
00423 { 00424 int i, j, n; 00425 register double *bptr; 00426 register float *Gxptr, *Gyptr, *Gzptr; 00427 double Gx, Gy, Gz; 00428 00429 ENTRY ("Computebmatrix"); 00430 n = grad1Dptr->nx; /* number of gradients other than I0 */ 00431 Gxptr = MRI_FLOAT_PTR (grad1Dptr); /* use simple floating point pointers to get values */ 00432 Gyptr = Gxptr + n; 00433 Gzptr = Gyptr + n; 00434 00435 bptr = bmatrix; 00436 for (i = 0; i < 6; i++) 00437 *bptr++ = 0.0; /* initialize first 6 elements to 0.0 for the I0 gradient */ 00438 00439 for (i = 0; i < n; i++) 00440 { 00441 Gx = *Gxptr++; 00442 Gy = *Gyptr++; 00443 Gz = *Gzptr++; 00444 *bptr++ = Gx * Gx; 00445 *bptr++ = Gx * Gy; 00446 *bptr++ = Gx * Gz; 00447 *bptr++ = Gy * Gy; 00448 *bptr++ = Gy * Gz; 00449 *bptr++ = Gz * Gz; 00450 for(j=0;j<6;j++) 00451 printf("%-13.6g ", *(bmatrix + 6 + (i*6)+ j) ); 00452 printf("\n"); 00453 } 00454 00455 00456 EXRETURN; 00457 } |
|
Definition at line 319 of file 3dDTtoDWI.c. References automask, bmatrix, ENTRY, i, and I0_ptr. Referenced by main().
00323 { 00324 int i, j; 00325 static int nvox, ncall; 00326 int allzeros; 00327 double I0, bq_d; 00328 double *bptr; 00329 float *tempptr; 00330 00331 ENTRY ("DTtoDWI_tsfunc"); 00332 /* ts is input vector data of 6 floating point numbers.*/ 00333 /* ts should come from data sub-briks in form of Dxx, Dxy, Dxz, Dyy, Dyz, Dzz */ 00334 /* if automask is turned on, ts has 7 floating point numbers */ 00335 /* val is output vector of form DWI0, DWI1, DWIn sub-briks */ 00336 /* where n = number of gradients = nbriks-1 */ 00337 00338 /** is this a "notification"? **/ 00339 if (val == NULL) 00340 { 00341 00342 if (npts > 0) 00343 { /* the "start notification" */ 00344 nvox = npts; /* keep track of */ 00345 ncall = 0; /* number of calls */ 00346 } 00347 else 00348 { /* the "end notification" */ 00349 00350 /* nothing to do here */ 00351 } 00352 EXRETURN; 00353 } 00354 00355 ncall++; 00356 00357 if (automask) 00358 npts = npts - 1; 00359 00360 tempptr = I0_ptr+ncall-1; 00361 I0 = *tempptr; 00362 00363 #if 0 00364 /* check for reasons not to calculate anything and just return zeros */ 00365 if(ncall==1) 00366 printf("checking for all zeros\n"); 00367 allzeros = 0; 00368 i=0; 00369 00370 while ((i<6)&&(I0!=0.0)) /* check if all the DT values are 0 */ 00371 { 00372 if(ts[i++]!=0) 00373 allzeros = 0; 00374 } 00375 00376 00377 for(i=0;i<nbriks;i++) 00378 val[i] = 0.0; 00379 EXRETURN; 00380 00381 00382 /* return zeros if all the DT values are 0, if the I0 value is 0 or 00383 the mask is off at this voxel*/ 00384 if(allzeros||(I0==0.0) || (automask&&(ts[npts] == 0))) 00385 { 00386 for (i = 0; i < nbriks; i++) 00387 val[i] = 0.0; /* return 0 for all DWIn points */ 00388 EXRETURN; 00389 } 00390 00391 #endif 00392 00393 val[0] = I0; /* the first sub-brik is the I0 sub-brik */ 00394 bptr = bmatrix+6; /* start at the first gradient */ 00395 00396 for(i=1;i<nbriks;i++) { 00397 bptr = bmatrix+(6*i); /* start at the first gradient */ 00398 bq_d = *bptr++ * ts[0]; /* GxGxDxx */ 00399 bq_d += *bptr++ * ts[1] * 2; /* 2GxGyDxy */ 00400 bq_d += *bptr++ * ts[2] * 2; /* 2GxGzDxz */ 00401 bq_d += *bptr++ * ts[3]; /* GyGyDyy */ 00402 bq_d += *bptr++ * ts[4] * 2; /* 2GyGzDyz */ 00403 bq_d += *bptr++ * ts[5]; /* GzGzDzz */ 00404 00405 val[i] = I0 * exp(-bq_d); /* for each gradient,q, Iq = J e -(bq.D) */ 00406 /* where bq.D is the large dot product of bq and D */ 00407 } 00408 00409 EXRETURN; 00410 } |
|
free up all the matrices and arrays Definition at line 475 of file 3dDTtoDWI.c. References bmatrix, ENTRY, and free. Referenced by main().
|
|
allocate all the global matrices and arrays once Definition at line 461 of file 3dDTtoDWI.c. References bmatrix, ENTRY, i, and malloc.
|
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 33 of file 3dDTtoDWI.c. References ADN_brick_label_one, ADN_none, ADN_ntt, ADN_ttdel, ADN_ttorg, ADN_tunits, AFNI_logger(), argc, automask, calloc, Computebmatrix(), datum, DSET_BRICK, DSET_BRICK_FACTOR, DSET_BRIKNAME, DSET_load, DSET_mallocize, DSET_NVALS, DSET_NVOX, DSET_unload_one, DSET_write, DTtoDWI_tsfunc(), EDIT_add_brick(), EDIT_dset_items(), free, FreeGlobals(), i, I0_ptr, InitGlobals(), ISVALID_DSET, MRI_IMAGE::kind, machdep(), mainENTRY, MAKER_4D_to_typed_fbuc(), MASTER_SHORTHELP_STRING, MCW_strncpy, mri_automask_image(), mri_data_pointer(), mri_free(), mri_read_1D(), MRI_IMAGE::nx, MRI_IMAGE::ny, prefix, PRINT_VERSION, THD_filename_ok(), THD_MAX_PREFIX, THD_open_dataset(), tross_Copy_History(), tross_Make_History(), and UNITS_SEC_TYPE.
00034 { 00035 THD_3dim_dataset *old_dset, *new_dset, *I0_dset; /* input and output datasets */ 00036 int nopt, nbriks, nvox; 00037 int i, eigs_brik; 00038 MRI_IMAGE *grad1Dptr = NULL; 00039 MRI_IMAGE *anat_im = NULL; 00040 MRI_IMAGE *data_im = NULL; 00041 double fac; 00042 short *sar = NULL, *tempsptr = NULL, tempval; 00043 byte *maskptr = NULL, *tempbptr = NULL; 00044 char tempstr[25]; 00045 00046 /*----- Read command line -----*/ 00047 if (argc < 2 || strcmp (argv[1], "-help") == 0) 00048 { 00049 printf ("Usage: 3dDTtoDWI [options] gradient-file I0-dataset DT-dataset\n" 00050 "Computes multiple gradient images from 6 principle direction tensors and\n" 00051 " corresponding gradient vector coordinates applied to the I0-dataset.\n" 00052 " The program takes three parameters as input : \n" 00053 " a 1D file of the gradient vectors with lines of ASCII floats Gxi,Gyi,Gzi.\n" 00054 " Only the non-zero gradient vectors are included in this file (no G0 line).\n" 00055 " The I0 dataset is a volume without any gradient applied.\n" 00056 " The DT dataset is the 6-sub-brick dataset containing the diffusion tensor data,\n" 00057 " Dxx, Dxy, Dxz, Dyy, Dyz, Dzz\n" 00058 " Options:\n" 00059 " -prefix pname = Use 'pname' for the output dataset prefix name.\n" 00060 " [default='DWI']\n" 00061 " -automask = mask dataset so that the gradient images are computed only for\n" 00062 " high-intensity (presumably brain) voxels. The intensity level is\n" 00063 " determined the same way that 3dClipLevel works.\n\n" 00064 " -datum type = output dataset type [float/short/byte] (default is float).\n" 00065 " -help = show this help screen.\n" 00066 " Example:\n" 00067 " 3dDTtoDWI -prefix DWI -automask tensor25.1D 'DT+orig[26]' DT+orig.\n\n" 00068 " The output is a n sub-brick bucket dataset containing computed DWI images.\n" 00069 " where n is the number of vectors in the gradient file + 1\n" 00070 "\n"); 00071 printf ("\n" MASTER_SHORTHELP_STRING); 00072 exit (0); 00073 } 00074 00075 mainENTRY ("3dDTtoDWI main"); 00076 machdep (); 00077 AFNI_logger ("3dDTtoDWI", argc, argv); 00078 PRINT_VERSION("3dDTtoDWI") ; 00079 00080 nopt = 1; 00081 datum = MRI_float; 00082 00083 00084 while (nopt < argc && argv[nopt][0] == '-') 00085 { 00086 00087 /*-- prefix --*/ 00088 00089 if (strcmp (argv[nopt], "-prefix") == 0) 00090 { 00091 if (++nopt >= argc) 00092 { 00093 fprintf (stderr, "*** Error - prefix needs an argument!\n"); 00094 exit (1); 00095 } 00096 MCW_strncpy (prefix, argv[nopt], THD_MAX_PREFIX); /* change name from default prefix */ 00097 /* check file name to be sure not to overwrite - mod drg 12/9/2004 */ 00098 if (!THD_filename_ok (prefix)) 00099 { 00100 fprintf (stderr, "*** Error - %s is not a valid prefix!\n", prefix); 00101 exit (1); 00102 } 00103 nopt++; 00104 continue; 00105 } 00106 00107 /*-- datum --*/ 00108 00109 if (strcmp (argv[nopt], "-datum") == 0) 00110 { 00111 if (++nopt >= argc) 00112 { 00113 fprintf (stderr, "*** Error - datum needs an argument!\n"); 00114 exit (1); 00115 } 00116 if (strcmp (argv[nopt], "short") == 0) 00117 { 00118 datum = MRI_short; 00119 } 00120 else if (strcmp (argv[nopt], "float") == 0) 00121 { 00122 datum = MRI_float; 00123 } 00124 else if (strcmp (argv[nopt], "byte") == 0) 00125 { 00126 datum = MRI_byte; 00127 } 00128 else 00129 { 00130 fprintf (stderr, "-datum of type '%s' is not supported!\n", 00131 argv[nopt]); 00132 exit (1); 00133 } 00134 nopt++; 00135 continue; 00136 } 00137 if (strcmp (argv[nopt], "-automask") == 0) 00138 { 00139 automask = 1; 00140 nopt++; 00141 continue; 00142 } 00143 00144 fprintf(stderr, "*** Error - unknown option %s\n", argv[nopt]); 00145 exit(1); 00146 } 00147 00148 /*----- read input datasets -----*/ 00149 00150 if (nopt >= argc) 00151 { 00152 fprintf (stderr, "*** Error - No input dataset!?\n"); 00153 exit (1); 00154 } 00155 00156 /* first input dataset - should be gradient vector file of ascii floats Gx,Gy,Gz */ 00157 00158 /* read gradient vector 1D file */ 00159 grad1Dptr = mri_read_1D (argv[nopt]); 00160 if (grad1Dptr == NULL) 00161 { 00162 fprintf (stderr, "*** Error reading gradient vector file\n"); 00163 exit (1); 00164 } 00165 00166 if (grad1Dptr->ny != 3) 00167 { 00168 fprintf (stderr, "*** Error - Only 3 columns of gradient vectors allowed\n"); 00169 fprintf (stderr, "%d columns found\n", grad1Dptr->nx); 00170 mri_free (grad1Dptr); 00171 exit (1); 00172 } 00173 00174 if (grad1Dptr->nx < 6) 00175 { 00176 fprintf (stderr, "*** Error - Must have at least 6 gradient vectors\n"); 00177 fprintf (stderr, "%d columns found\n", grad1Dptr->nx); 00178 mri_free (grad1Dptr); 00179 exit (1); 00180 } 00181 00182 nbriks = grad1Dptr->nx + 1; /* number of gradients specified here from file */ 00183 nopt++; 00184 00185 /* open I0 dataset - idealized no gradient image */ 00186 I0_dset = THD_open_dataset (argv[nopt]); 00187 if (!ISVALID_DSET (I0_dset)) 00188 { 00189 fprintf (stderr, "*** Error - Can not open dataset %s\n", argv[nopt]); 00190 mri_free (grad1Dptr); 00191 exit (1); 00192 } 00193 00194 DSET_mallocize (I0_dset); 00195 DSET_load (I0_dset); /* load dataset */ 00196 data_im = DSET_BRICK (I0_dset, 0); /* set pointer to the 0th sub-brik of the dataset */ 00197 fac = DSET_BRICK_FACTOR(I0_dset, 0); /* get scale factor for each sub-brik*/ 00198 if(fac==0.0) fac=1.0; 00199 if((data_im->kind != MRI_float) || (fac!=1.0)) { 00200 fprintf (stderr, "*** Error - Can only open float datasets with scale factors of 1\n"); 00201 mri_free (grad1Dptr); 00202 mri_free (data_im); 00203 exit (1); 00204 } 00205 00206 I0_ptr = mri_data_pointer(data_im) ; /* pointer to I0 data */ 00207 00208 nopt++; 00209 00210 /* Now read in all the MRI volumes for each gradient vector */ 00211 /* assumes first one is no gradient */ 00212 old_dset = THD_open_dataset (argv[nopt]); 00213 00214 if (!ISVALID_DSET (old_dset)) 00215 { 00216 fprintf (stderr, "*** Error - Can not open dataset %s\n", argv[nopt]); 00217 mri_free (grad1Dptr); 00218 exit (1); 00219 } 00220 00221 /* expect at least 6 values per voxel - 6 sub-briks as input dataset */ 00222 if (DSET_NVALS (old_dset) <6) 00223 { 00224 fprintf (stderr, 00225 "*** Error - Dataset must have at least 6 sub-briks to describe the diffusion tensor\n"); 00226 mri_free (grad1Dptr); 00227 mri_free (data_im); 00228 exit (1); 00229 } 00230 00231 00232 InitGlobals (grad1Dptr->nx + 1); /* initialize all the matrices and vectors */ 00233 Computebmatrix (grad1Dptr); /* compute bij=GiGj */ 00234 00235 if (automask) 00236 { 00237 DSET_mallocize (old_dset); 00238 DSET_load (old_dset); /* get B0 (anatomical image) from dataset */ 00239 /*anat_im = THD_extract_float_brick( 0, old_dset ); */ 00240 anat_im = DSET_BRICK (old_dset, 0); /* set pointer to the 0th sub-brik of the dataset */ 00241 maskptr = mri_automask_image (anat_im); /* maskptr is a byte pointer for volume */ 00242 00243 /* convert byte mask to same format type as dataset */ 00244 nvox = DSET_NVOX (old_dset); 00245 sar = (short *) calloc (nvox, sizeof (short)); 00246 /* copy maskptr values to far ptr */ 00247 tempsptr = sar; 00248 tempbptr = maskptr; 00249 for (i = 0; i < nvox; i++) 00250 { 00251 *tempsptr++ = (short) *tempbptr++; 00252 tempval = *(tempsptr - 1); 00253 } 00254 00255 free (maskptr); 00256 00257 /*old_dset->dblk->malloc_type = DATABLOCK_MEM_MALLOC; *//* had to set this? */ 00258 EDIT_add_brick (old_dset, MRI_short, 0.0, sar); /* add sub-brik to end */ 00259 00260 00261 } 00262 00263 /* temporarily set artificial timing to 1 second interval */ 00264 EDIT_dset_items (old_dset, 00265 ADN_ntt, DSET_NVALS (old_dset), 00266 ADN_ttorg, 0.0, 00267 ADN_ttdel, 1.0, ADN_tunits, UNITS_SEC_TYPE, NULL); 00268 00269 /*------------- ready to compute new dataset -----------*/ 00270 00271 new_dset = MAKER_4D_to_typed_fbuc (old_dset, /* input dataset */ 00272 prefix, /* output prefix */ 00273 datum, /* output datum */ 00274 0, /* ignore count */ 00275 0, /* can't detrend in maker function KRH 12/02 */ 00276 nbriks, /* number of briks */ 00277 DTtoDWI_tsfunc, /* timeseries processor */ 00278 NULL /* data for tsfunc */ 00279 ); 00280 00281 00282 00283 FreeGlobals (); 00284 mri_free (grad1Dptr); 00285 00286 00287 if (automask) 00288 { 00289 mri_free (anat_im); 00290 DSET_unload_one (old_dset, 0); 00291 sar = NULL; 00292 } 00293 00294 if (new_dset != NULL) 00295 { 00296 tross_Copy_History (old_dset, new_dset); 00297 for(i=0;i<nbriks;i++) { 00298 sprintf(tempstr,"grad%3.3d", i); 00299 EDIT_dset_items (new_dset, ADN_brick_label_one + i, tempstr, ADN_none); 00300 } 00301 tross_Make_History ("3dDTtoDWI", argc, argv, new_dset); 00302 DSET_write (new_dset); 00303 fprintf(stderr,"--- Output dataset %s\n", DSET_BRIKNAME(new_dset)); 00304 } 00305 else 00306 { 00307 fprintf (stderr, "*** Error - Unable to compute output dataset!\n"); 00308 exit (1); 00309 } 00310 00311 exit (0); 00312 } |
Variable Documentation
|
Definition at line 25 of file 3dDTtoDWI.c. Referenced by DTtoDWI_tsfunc(), and main(). |
|
Definition at line 22 of file 3dDTtoDWI.c. Referenced by Computebmatrix(), DTtoDWI_tsfunc(), FreeGlobals(), and InitGlobals(). |
|
Definition at line 18 of file 3dDTtoDWI.c. Referenced by main(). |
|
Definition at line 19 of file 3dDTtoDWI.c. |
|
Definition at line 20 of file 3dDTtoDWI.c. |
|
Definition at line 24 of file 3dDTtoDWI.c. Referenced by DTtoDWI_tsfunc(), and main(). |
|
Definition at line 17 of file 3dDTtoDWI.c. Referenced by main(). |