1 /*********************** 3dDTtoDWI.c **********************************************/
2 /* Author: Daniel Glen, 17 Nov 2004 */
3 /* compute multiple gradient images from 6 principle direction tensors*/
4 /* and corresponding gradient vector coordinates */
5 /* used for modeling diffusion tensor data */
6
7 // [May, 2017] input option to complement 3dDWItoDT's -scale_out_1000
8 // option. It has the same name, funnily enough!
9
10 // [May, 2017]
11 // + bug fix: had been a mismatch between bmatrix and DT. Now fixed.
12
13 #include "thd_shear3d.h"
14 # include "matrix.h"
15 #include "mrilib.h"
16
17 #define TINYNUMBER 1E-10
18 #define SMALLNUMBER 1E-4
19
20 static char prefix[THD_MAX_PREFIX] = "DWI";
21 static int datum = MRI_float;
22
23 static double *bmatrix; /* b matrix = GiGj for gradient intensities */
24 static int BMAT_NZ = 0; /* non-zero bmatrix supplied: apr,2016: not
25 doing anything here at the moment, just
26 in case we want to start using bmatrices
27 as input instead of grads, such as in
28 3dDWItoDT. */
29 static float MAX_BVAL = 1.;
30 static double *B0list = NULL;
31
32 static float *I0_ptr; /* I0 matrix */
33 static int automask = 0; /* automasking flag - user option */
34
35 static void DTtoDWI_tsfunc (double tzero, double tdelta, int npts, float ts[],
36 double ts_mean, double ts_slope, void *ud,
37 int nbriks, float *val);
38 // apr,2016: copying/pasting this function now from 3dDWItoDT
39 static void Computebmatrix (MRI_IMAGE * grad1Dptr, int NO_ZERO_ROW1);
40 static void InitGlobals (int npts);
41 static void FreeGlobals (void);
42
43 int
main(int argc,char * argv[])44 main (int argc, char *argv[])
45 {
46 THD_3dim_dataset *old_dset, *new_dset, *I0_dset; /* input and output datasets */
47 int nopt, nbriks, nvox;
48 int i;
49 MRI_IMAGE *grad1Dptr = NULL;
50 MRI_IMAGE *anat_im = NULL;
51 MRI_IMAGE *data_im = NULL;
52 double fac;
53 short *sar = NULL, *tempsptr = NULL, tempval;
54 byte *maskptr = NULL, *tempbptr = NULL;
55 char tempstr[25];
56
57 float SCALE_VAL_OUT = -1.0 ; // allow users to scaled physical
58 // values by 1000 easily; will be
59 // set to 1 if negative after
60 // reading in inputs
61
62 /*----- Read command line -----*/
63 if (argc < 2 || strcmp (argv[1], "-help") == 0) {
64 printf (
65 "Usage: 3dDTtoDWI [options] gradient-file I0-dataset DT-dataset\n"
66 "\n"
67 "Computes multiple gradient images from 6 principle direction tensors and\n"
68 " corresponding gradient vector coordinates applied to the I0-dataset.\n"
69 " The program takes three parameters as input : \n"
70 " a 1D file of the gradient vectors with lines of ASCII floats Gxi,Gyi,Gzi.\n"
71 " Only the non-zero gradient vectors are included in this file (no G0 line).\n"
72 " The I0 dataset is a volume without any gradient applied.\n"
73 " The DT dataset is the 6-sub-brick dataset containing the diffusion tensor data,\n"
74 " Dxx, Dxy, Dyy, Dxz, Dyz, Dzz (lower triangular row-wise order)\n"
75 "\n"
76 " Options:\n"
77 " -prefix pname = Use 'pname' for the output dataset prefix name.\n"
78 " [default='DWI']\n"
79 " -automask = mask dataset so that the gradient images\n"
80 " are computed only for high-intensity (presumably\n"
81 " brain) voxels. The intensity level is determined\n"
82 " the same way that 3dClipLevel works.\n\n"
83 " -datum type = output dataset type [float/short/byte] \n"
84 " (default is float).\n"
85 " -help = show this help screen.\n"
86 " -scale_out_1000 = matches with 3dDWItoDT's '-scale_out_1000'\n"
87 " functionality. If the option was used\n"
88 " there, then use it here, too.\n"
89 "\n"
90 " Example:\n"
91 "\n"
92 " 3dDTtoDWI -prefix DWI -automask tensor25.1D 'DT+orig[26]' DT+orig.\n\n"
93 "\n"
94 " The output is a n sub-brick bucket dataset containing computed DWI images.\n"
95 " where n is the number of vectors in the gradient file + 1\n"
96 "\n"
97 );
98 printf ("\n" MASTER_SHORTHELP_STRING);
99 exit (0);
100 }
101
102 mainENTRY ("3dDTtoDWI main");
103 machdep ();
104 AFNI_logger ("3dDTtoDWI", argc, argv);
105 PRINT_VERSION("3dDTtoDWI") ;
106
107 nopt = 1;
108 datum = MRI_float;
109
110
111 while (nopt < argc && argv[nopt][0] == '-')
112 {
113
114 /*-- prefix --*/
115
116 if (strcmp (argv[nopt], "-prefix") == 0)
117 {
118 if (++nopt >= argc) {
119 fprintf (stderr, "*** Error - prefix needs an argument!\n");
120 exit (1);
121 }
122 MCW_strncpy (prefix, argv[nopt], THD_MAX_PREFIX); /* change name from default prefix */
123 /* check file name to be sure not to overwrite - mod drg 12/9/2004 */
124 if (!THD_filename_ok (prefix))
125 {
126 fprintf (stderr, "*** Error - %s is not a valid prefix!\n", prefix);
127 exit (1);
128 }
129 nopt++;
130 continue;
131 }
132
133 /*-- datum --*/
134
135 if (strcmp (argv[nopt], "-datum") == 0)
136 {
137 if (++nopt >= argc) {
138 fprintf (stderr, "*** Error - datum needs an argument!\n");
139 exit (1);
140 }
141 if (strcmp (argv[nopt], "short") == 0) {
142 datum = MRI_short;
143 }
144 else if (strcmp (argv[nopt], "float") == 0) {
145 datum = MRI_float;
146 }
147 else if (strcmp (argv[nopt], "byte") == 0) {
148 datum = MRI_byte;
149 }
150 else {
151 fprintf (stderr, "-datum of type '%s' is not supported!\n",
152 argv[nopt]);
153 exit (1);
154 }
155 nopt++;
156 continue;
157 }
158 if (strcmp (argv[nopt], "-automask") == 0) {
159 automask = 1;
160 nopt++;
161 continue;
162 }
163
164 // can mult output diff values by 1000
165 if (strcmp(argv[nopt], "-scale_out_1000") == 0) {
166 SCALE_VAL_OUT = 1000.;
167 nopt++;
168 continue;
169 }
170
171 fprintf(stderr, "*** Error - unknown option %s\n", argv[nopt]);
172 exit(1);
173 }
174
175 // ==============================================================
176 // ====================== begin processing ======================
177
178 /*----- read input datasets -----*/
179
180 if (nopt >= argc) {
181 fprintf (stderr, "*** Error - No input dataset!?\n");
182 exit (1);
183 }
184
185 if( SCALE_VAL_OUT < 0 )
186 SCALE_VAL_OUT = 1.;
187 else
188 INFO_message("Implementing user-selected scaling of diffusion "
189 "values by %.0f",
190 SCALE_VAL_OUT);
191
192 /* first input dataset - should be gradient vector file of ascii
193 floats Gx,Gy,Gz */
194
195 /* read gradient vector 1D file */
196 grad1Dptr = mri_read_1D (argv[nopt]);
197 if (grad1Dptr == NULL)
198 {
199 fprintf (stderr, "*** Error reading gradient vector file\n");
200 exit (1);
201 }
202
203 if (grad1Dptr->ny != 3) {
204 fprintf (stderr, "*** Error - Only 3 columns of gradient vectors allowed\n");
205 fprintf (stderr, "%d columns found\n", grad1Dptr->nx);
206 mri_free (grad1Dptr);
207 exit (1);
208 }
209
210 if (grad1Dptr->nx < 6) {
211 fprintf (stderr, "*** Error - Must have at least 6 gradient vectors\n");
212 fprintf (stderr, "%d columns found\n", grad1Dptr->nx);
213 mri_free (grad1Dptr);
214 exit (1);
215 }
216
217 nbriks = grad1Dptr->nx + 1; /* number of gradients specified here from file */
218 nopt++;
219
220 // [PT: May, 2017] copied from 3dDWItoDT to apply the scaling
221 int npix;
222 register float *flar=mri_data_pointer(grad1Dptr);
223 npix = grad1Dptr->nvox;
224 for( i=0 ; i<npix ; i++ ){
225 flar[i]/= SCALE_VAL_OUT;
226 }
227
228 /* open I0 dataset - idealized no gradient image */
229 I0_dset = THD_open_dataset (argv[nopt]);
230 CHECK_OPEN_ERROR(I0_dset,argv[nopt]) ;
231
232 DSET_mallocize (I0_dset);
233 DSET_load (I0_dset); /* load dataset */
234 data_im = DSET_BRICK (I0_dset, 0); /* set pointer to the 0th sub-brik of the dataset */
235 fac = DSET_BRICK_FACTOR(I0_dset, 0); /* get scale factor for each sub-brik*/
236 if(fac==0.0) fac=1.0;
237 if((data_im->kind != MRI_float)) {
238 fprintf (stderr, "*** Error - Can only open float datasets. Use 3dcalc to convert.\n");
239 mri_free (grad1Dptr);
240 mri_free (data_im);
241 exit (1);
242 }
243
244 I0_ptr = mri_data_pointer(data_im) ; /* pointer to I0 data */
245 if( I0_ptr == NULL ) {
246 ERROR_exit("Null pointer when trying to make I0 set.");
247 }
248
249
250 nopt++;
251
252 /* Now read in all the MRI volumes for each gradient vector */
253 /* assumes first one is no gradient */
254 old_dset = THD_open_dataset (argv[nopt]);
255 CHECK_OPEN_ERROR(old_dset,argv[nopt]) ;
256
257 /* expect at least 6 values per voxel - 6 sub-briks as input dataset */
258 if (DSET_NVALS (old_dset) <6)
259 {
260 fprintf (stderr,
261 "*** Error - Dataset must have at least 6 sub-briks to describe the diffusion tensor\n");
262 mri_free (grad1Dptr);
263 mri_free (data_im);
264 exit (1);
265 }
266
267
268 InitGlobals (grad1Dptr->nx + 1); /* initialize all the matrices and vectors */
269 Computebmatrix (grad1Dptr, BMAT_NZ); /* compute bij=GiGj */
270 INFO_message("The maximum magnitude of the bmatrix appears to be: %.2f", MAX_BVAL);
271
272 if (automask)
273 {
274 DSET_mallocize (old_dset);
275 DSET_load (old_dset); /* get B0 (anatomical image) from dataset */
276 /*anat_im = THD_extract_float_brick( 0, old_dset ); */
277 anat_im = DSET_BRICK (old_dset, 0); /* set pointer to the 0th sub-brik of the dataset */
278 maskptr = mri_automask_image (anat_im); /* maskptr is a byte pointer for volume */
279
280 /* convert byte mask to same format type as dataset */
281 nvox = DSET_NVOX (old_dset);
282 sar = (short *) calloc (nvox, sizeof (short));
283 /* copy maskptr values to far ptr */
284 tempsptr = sar;
285 tempbptr = maskptr;
286 for (i = 0; i < nvox; i++)
287 {
288 *tempsptr++ = (short) *tempbptr++;
289 tempval = *(tempsptr - 1);
290 }
291
292 free (maskptr);
293
294 /*old_dset->dblk->malloc_type = DATABLOCK_MEM_MALLOC; *//* had to set this? */
295 EDIT_add_brick (old_dset, MRI_short, 0.0, sar); /* add sub-brik to end */
296
297
298 }
299
300 /* temporarily set artificial timing to 1 second interval */
301 EDIT_dset_items (old_dset,
302 ADN_ntt, DSET_NVALS (old_dset),
303 ADN_ttorg, 0.0,
304 ADN_ttdel, 1.0, ADN_tunits, UNITS_SEC_TYPE, NULL);
305
306 /*------------- ready to compute new dataset -----------*/
307
308 new_dset = MAKER_4D_to_typed_fbuc (old_dset, /* input dataset */
309 prefix, /* output prefix */
310 datum, /* output datum */
311 0, /* ignore count */
312 0, /* can't detrend in maker function KRH 12/02 */
313 nbriks, /* number of briks */
314 DTtoDWI_tsfunc, /* timeseries processor */
315 NULL, /* data for tsfunc */
316 NULL, /* mask */
317 0 /* Allow auto scaling of output */
318 );
319
320
321
322 FreeGlobals ();
323 mri_free (grad1Dptr);
324
325
326 if (automask)
327 {
328 mri_free (anat_im);
329 DSET_unload_one (old_dset, 0);
330 sar = NULL;
331 }
332
333 if (new_dset != NULL)
334 {
335 tross_Copy_History (old_dset, new_dset);
336 for(i=0;i<nbriks;i++) {
337 sprintf(tempstr,"grad%3.3d", i);
338 EDIT_dset_items (new_dset, ADN_brick_label_one + i, tempstr, ADN_none);
339 }
340 tross_Make_History ("3dDTtoDWI", argc, argv, new_dset);
341 DSET_write (new_dset);
342 fprintf(stderr,"--- Output dataset %s\n", DSET_BRIKNAME(new_dset));
343 }
344 else
345 {
346 fprintf (stderr, "*** Error - Unable to compute output dataset!\n");
347 exit (1);
348 }
349
350 exit (0);
351 }
352
353 /**********************************************************************
354 Function that does the real work
355 ***********************************************************************/
356
357 /*
358 [PT: May, 2017]
359 !! NB !! This has been copy/pasted *to* 3dDWItoDT.c, so if changes
360 happen there/here, likely make them here/there as well.
361 This inglorious route was chosen rather than deal with
362 compile-time issues with static variables. I'm sure I'll
363 adjust this more properly later. Definitely.
364 */
365 static void
DTtoDWI_tsfunc(double tzero,double tdelta,int npts,float ts[],double ts_mean,double ts_slope,void * ud,int nbriks,float * val)366 DTtoDWI_tsfunc (double tzero, double tdelta,
367 int npts, float ts[],
368 double ts_mean, double ts_slope,
369 void *ud, int nbriks, float *val)
370 {
371 int i;
372 static int nvox2, ncall2;
373 double I0, bq_d;
374 double *bptr;
375 float *tempptr;
376
377 ENTRY ("DTtoDWI_tsfunc");
378 /* ts is input vector data of 6 floating point numbers.*/
379 /* ts should come from data sub-briks in form of Dxx, Dxy, Dyy, Dxz, Dyz, Dzz */
380 /* if automask is turned on, ts has 7 floating point numbers */
381 /* val is output vector of form DWI0, DWI1, DWIn sub-briks */
382 /* where n = number of gradients = nbriks-1 */
383
384 /** is this a "notification"? **/
385 if (val == NULL)
386 {
387
388 if (npts > 0) { /* the "start notification" */
389 nvox2 = npts; /* keep track of */
390 ncall2 = 0; /* number of calls */
391 }
392 else { /* the "end notification" */
393
394 /* nothing to do here */
395 }
396 EXRETURN;
397 }
398
399 ncall2++;
400
401 if (automask)
402 npts = npts - 1;
403
404 tempptr = I0_ptr+ncall2-1;
405 I0 = *tempptr;
406
407
408 val[0] = I0; /* the first sub-brik is the I0 sub-brik */
409 bptr = bmatrix+6; /* start at the first gradient */
410
411 /* NB: what the ordering of bmatrix and DT are:
412 // bmatr
413 [0] *bptr++ = Gx * Gx / gscale;
414 [1] *bptr++ = Gx * Gy / gscale;
415 [2] *bptr++ = Gx * Gz / gscale;
416 [3] *bptr++ = Gy * Gy / gscale;
417 [4] *bptr++ = Gy * Gz / gscale;
418 [5] *bptr++ = Gz * Gz / gscale;
419
420 [0] Dxx,
421 [1] Dxy,
422 [2] Dyy,
423 [3] Dxz,
424 [4] Dyz,
425 [5] Dzz
426 */
427
428 // Need to match bmatrix and DT orderings, which are different.
429 // Yay.
430 for(i=1;i<nbriks;i++) {
431 bptr = bmatrix+(6*i); // start at the first gradient
432 bq_d = *bptr++ * ts[0]; // GxGxDxx
433 bq_d += *bptr++ * ts[1] * 2; // 2GxGyDxy
434 bq_d += *bptr++ * ts[3] * 2; // 2GxGzDxz
435 bq_d += *bptr++ * ts[2]; // GyGyDyy
436 bq_d += *bptr++ * ts[4] * 2; // 2GyGzDyz
437 bq_d += *bptr++ * ts[5]; // GzGzDzz
438
439 // for each gradient,q, Iq = J e -(bq.D)
440 val[i] = I0 * exp(-bq_d);
441 // where bq.D is the large dot product of bq and D
442 }
443
444 EXRETURN;
445 }
446
447
448 /*! compute the diffusion weighting matrix bmatrix for q number of gradients */
449 /* only need to calculate once */
450 /* bq = diffusion weighting matrix of qth gradient */
451 /* GxGx GxGy GxGz
452 GxGy GyGy GyGz
453 GxGz GyGz GzGz */
454 /* b0 is 0 for all 9 elements */
455 /* bmatrix is really stored as 6 x npts array */
456 // ----> apr,2016: now copying the Computebmatrix(...) from 3dDWItoDT,
457 // ----> to have the scale stuff
458 static void
Computebmatrix(MRI_IMAGE * grad1Dptr,int NO_ZERO_ROW1)459 Computebmatrix (MRI_IMAGE * grad1Dptr, int NO_ZERO_ROW1)
460 {
461 int i, n;
462 register double *bptr;
463 register float *Gxptr, *Gyptr, *Gzptr;
464 register float *Bxxptr, *Byyptr, *Bzzptr, *Bxyptr, *Bxzptr, *Byzptr;
465 double Gx, Gy, Gz, Bxx, Byy, Bzz, Bxy, Bxz, Byz;
466 double gscale;
467 int j;
468 ENTRY ("Computebmatrix");
469 n = grad1Dptr->nx; /* number of gradients other than I0 */
470
471 if ( (grad1Dptr->ny == 6) && !NO_ZERO_ROW1 ) { // extra switch to keep OLD, zero-row version
472 /* just read in b matrix */
473 Bxxptr = MRI_FLOAT_PTR (grad1Dptr); /* use simple floating point pointers to get values */
474 Byyptr = Bxxptr + n;
475 Bzzptr = Byyptr + n;
476 Bxyptr = Bzzptr + n;
477 Bxzptr = Bxyptr + n;
478 Byzptr = Bxzptr + n;
479
480 bptr = bmatrix;
481
482 /* B0list[0]= 1;*/ /* keep a record of which volumes have no gradient: first one always assumed */
483
484 for (i = 0; i < n; i++){
485 Bxx = *Bxxptr++;
486 Byy = *Byyptr++;
487 Bzz = *Bzzptr++;
488 Bxy = *Bxyptr++;
489 Bxz = *Bxzptr++;
490 Byz = *Byzptr++;
491 *bptr++ = Bxx;
492 *bptr++ = Bxy;
493 *bptr++ = Bxz;
494 *bptr++ = Byy;
495 *bptr++ = Byz;
496 *bptr++ = Bzz;
497 if(Bxx==0.0 && Byy==0.0 && Bzz==0.0) /* is this a zero gradient volume also? */
498 B0list[i] = 1;
499 else{
500 B0list[i] = 0;
501 // apr,2016: need the MAX_BVAL for scaling
502 gscale = Bxx + Byy + Bzz;
503 if(gscale > MAX_BVAL)
504 MAX_BVAL = gscale;
505 }
506
507 }
508 }
509 else if( (grad1Dptr->ny == 6) && NO_ZERO_ROW1 ) { // very similar to old bmatrix option
510 /* just read in b matrix */
511 Bxxptr = MRI_FLOAT_PTR (grad1Dptr); /* use simple floating point pointers to get values */
512 Byyptr = Bxxptr + n;
513 Bzzptr = Byyptr + n;
514 Bxyptr = Bzzptr + n;
515 Bxzptr = Bxyptr + n;
516 Byzptr = Bxzptr + n;
517
518 bptr = bmatrix;
519
520 // do as grads below
521 for (i = 0; i < 6; i++)
522 *bptr++ = 0.0; /* initialize first 6 elements to 0.0 for the I0 gradient */
523 B0list[0]= 1; /* keep a record of which volumes have no gradient */
524
525
526 for (i = 0; i < n; i++){
527 Bxx = *Bxxptr++;
528 Byy = *Byyptr++;
529 Bzz = *Bzzptr++;
530 Bxy = *Bxyptr++;
531 Bxz = *Bxzptr++;
532 Byz = *Byzptr++;
533 *bptr++ = Bxx;
534 *bptr++ = Bxy;
535 *bptr++ = Bxz;
536 *bptr++ = Byy;
537 *bptr++ = Byz;
538 *bptr++ = Bzz;
539 if(Bxx==0.0 && Byy==0.0 && Bzz==0.0) /* is this a zero gradient volume also? */
540 B0list[i+1] = 1;
541 else{
542 B0list[i+1] = 0;
543 // apr,2016: need the MAX_BVAL for scaling
544 gscale = Bxx + Byy + Bzz;
545 if(gscale > MAX_BVAL)
546 MAX_BVAL = gscale;
547 }
548
549 }
550 }
551 else {
552 Gxptr = MRI_FLOAT_PTR (grad1Dptr); /* use simple floating point pointers to get values */
553 Gyptr = Gxptr + n;
554 Gzptr = Gyptr + n;
555
556 bptr = bmatrix;
557 for (i = 0; i < 6; i++)
558 *bptr++ = 0.0; /* initialize first 6 elements to 0.0 for the I0 gradient */
559
560 B0list[0]= 1; /* keep a record of which volumes have no gradient */
561
562 for (i = 0; i < n; i++)
563 {
564 gscale = 1.; // apr,2016: allow for non-unit gradient magnitudes
565 Gx = *Gxptr++;
566 Gy = *Gyptr++;
567 Gz = *Gzptr++;
568 if((Gx==0.0) && (Gy==0.0) && (Gz==0.0))
569 B0list[i+1] = 1; /* no gradient applied*/
570 else{
571 B0list[i+1] = 0;
572 gscale = sqrt(Gx*Gx + Gy*Gy + Gz*Gz);
573 if(gscale > MAX_BVAL)
574 MAX_BVAL = gscale; // apr,2016
575 }
576
577 *bptr++ = Gx * Gx / gscale; // apr,2016: allow for non-unit gradient magnitudes
578 *bptr++ = Gx * Gy / gscale;
579 *bptr++ = Gx * Gz / gscale;
580 *bptr++ = Gy * Gy / gscale;
581 *bptr++ = Gy * Gz / gscale;
582 *bptr++ = Gz * Gz / gscale;
583
584 for(j=0;j<6;j++)
585 printf(" %15.5f ", *(bmatrix + 6 + (i*6)+ j) );
586 printf("\n");
587 }
588 }
589 EXRETURN;
590 }
591 /* --->apr,2016: This is the **original** function for this program-- now using
592 the one from 3dDWItoDT, for each of duplicity. Or something.
593 {
594 int i, j, n;
595 register double *bptr;
596 register float *Gxptr, *Gyptr, *Gzptr;
597 double Gx, Gy, Gz;
598
599 ENTRY ("Computebmatrix");
600 n = grad1Dptr->nx; // number of gradients other than I0
601 Gxptr = MRI_FLOAT_PTR (grad1Dptr); // use simple floating point pointers to get values
602 Gyptr = Gxptr + n;
603 Gzptr = Gyptr + n;
604
605 bptr = bmatrix;
606 for (i = 0; i < 6; i++)
607 *bptr++ = 0.0; // initialize first 6 elements to 0.0 for the I0 gradient
608
609 for (i = 0; i < n; i++)
610 {
611 Gx = *Gxptr++;
612 Gy = *Gyptr++;
613 Gz = *Gzptr++;
614 *bptr++ = Gx * Gx;
615 *bptr++ = Gx * Gy;
616 *bptr++ = Gy * Gy;
617 *bptr++ = Gx * Gz;
618 *bptr++ = Gy * Gz;
619 *bptr++ = Gz * Gz;
620 for(j=0;j<6;j++)
621 printf("%-13.6g ", *(bmatrix + 6 + (i*6)+ j) );
622 printf("\n");
623 }
624
625
626 EXRETURN;
627 }*/
628
629 /*! allocate all the global matrices and arrays once */
630 static void
InitGlobals(int npts)631 InitGlobals (int npts)
632 {
633
634 ENTRY ("InitGlobals");
635
636 bmatrix = malloc (npts * 6 * sizeof (double));
637 B0list = malloc (npts * sizeof (double));
638
639 EXRETURN;
640 }
641
642 /*! free up all the matrices and arrays */
643 static void
FreeGlobals()644 FreeGlobals ()
645 {
646
647 ENTRY ("FreeGlobals");
648 free (bmatrix);
649 bmatrix = NULL;
650
651 free(B0list);
652 B0list = NULL;
653
654
655 EXRETURN;
656 }
657