00001
00002
00003
00004
00005
00006
00007 #include "thd_shear3d.h"
00008
00009 # include "matrix.h"
00010 # include "matrix.c"
00011
00012 #include "afni.h"
00013
00014 #define TINYNUMBER 1E-10
00015 #define SMALLNUMBER 1E-4
00016
00017 static char prefix[THD_MAX_PREFIX] = "DT";
00018 static int datum = MRI_float;
00019 static matrix Dmatrix;
00020 static vector Dvector;
00021
00022 static double *bmatrix;
00023
00024 static float *I0_ptr;
00025 static int automask = 0;
00026
00027 static void DTtoDWI_tsfunc (double tzero, double tdelta, int npts, float ts[], double ts_mean, double ts_slope, void *ud, int nbriks, float *val);
00028 static void Computebmatrix (MRI_IMAGE * grad1Dptr);
00029 static void InitGlobals (int npts);
00030 static void FreeGlobals (void);
00031
00032 int
00033 main (int argc, char *argv[])
00034 {
00035 THD_3dim_dataset *old_dset, *new_dset, *I0_dset;
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
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
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);
00097
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
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
00149
00150 if (nopt >= argc)
00151 {
00152 fprintf (stderr, "*** Error - No input dataset!?\n");
00153 exit (1);
00154 }
00155
00156
00157
00158
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;
00183 nopt++;
00184
00185
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);
00196 data_im = DSET_BRICK (I0_dset, 0);
00197 fac = DSET_BRICK_FACTOR(I0_dset, 0);
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) ;
00207
00208 nopt++;
00209
00210
00211
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
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);
00233 Computebmatrix (grad1Dptr);
00234
00235 if (automask)
00236 {
00237 DSET_mallocize (old_dset);
00238 DSET_load (old_dset);
00239
00240 anat_im = DSET_BRICK (old_dset, 0);
00241 maskptr = mri_automask_image (anat_im);
00242
00243
00244 nvox = DSET_NVOX (old_dset);
00245 sar = (short *) calloc (nvox, sizeof (short));
00246
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
00258 EDIT_add_brick (old_dset, MRI_short, 0.0, sar);
00259
00260
00261 }
00262
00263
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
00270
00271 new_dset = MAKER_4D_to_typed_fbuc (old_dset,
00272 prefix,
00273 datum,
00274 0,
00275 0,
00276 nbriks,
00277 DTtoDWI_tsfunc,
00278 NULL
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 }
00313
00314
00315
00316
00317
00318 static void
00319 DTtoDWI_tsfunc (double tzero, double tdelta,
00320 int npts, float ts[],
00321 double ts_mean, double ts_slope,
00322 void *ud, int nbriks, float *val)
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
00333
00334
00335
00336
00337
00338
00339 if (val == NULL)
00340 {
00341
00342 if (npts > 0)
00343 {
00344 nvox = npts;
00345 ncall = 0;
00346 }
00347 else
00348 {
00349
00350
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
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))
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
00383
00384 if(allzeros||(I0==0.0) || (automask&&(ts[npts] == 0)))
00385 {
00386 for (i = 0; i < nbriks; i++)
00387 val[i] = 0.0;
00388 EXRETURN;
00389 }
00390
00391 #endif
00392
00393 val[0] = I0;
00394 bptr = bmatrix+6;
00395
00396 for(i=1;i<nbriks;i++) {
00397 bptr = bmatrix+(6*i);
00398 bq_d = *bptr++ * ts[0];
00399 bq_d += *bptr++ * ts[1] * 2;
00400 bq_d += *bptr++ * ts[2] * 2;
00401 bq_d += *bptr++ * ts[3];
00402 bq_d += *bptr++ * ts[4] * 2;
00403 bq_d += *bptr++ * ts[5];
00404
00405 val[i] = I0 * exp(-bq_d);
00406
00407 }
00408
00409 EXRETURN;
00410 }
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421 static void
00422 Computebmatrix (MRI_IMAGE * grad1Dptr)
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;
00431 Gxptr = MRI_FLOAT_PTR (grad1Dptr);
00432 Gyptr = Gxptr + n;
00433 Gzptr = Gyptr + n;
00434
00435 bptr = bmatrix;
00436 for (i = 0; i < 6; i++)
00437 *bptr++ = 0.0;
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 }
00458
00459
00460 static void
00461 InitGlobals (int npts)
00462 {
00463 int i;
00464 double *cumulativewtptr;
00465
00466 ENTRY ("InitGlobals");
00467
00468 bmatrix = malloc (npts * 6 * sizeof (double));
00469
00470 EXRETURN;
00471 }
00472
00473
00474 static void
00475 FreeGlobals ()
00476 {
00477
00478 ENTRY ("FreeGlobals");
00479 free (bmatrix);
00480 bmatrix = NULL;
00481 EXRETURN;
00482 }