~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/medical/medStim_gamma.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU Lesser General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU Lesser General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU Lesser General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  **********
 20  *
 21  * Program :    TINA
 22  * File    :  $Source: medStim_gamma.c $
 23  * Date    :  $Date: 2013/01/28 11:10 $
 24  * Version :  $Revision: 1.14 $
 25  *
 26  * Author  : Legacy TINA modified NAT/HR
 27  *
 28  * Notes:
 29  * Xiaoping
 30  * a.lacey   13.5.98
 31  * NAT       24.5.98->11.11.2002 
 32  *
 33  *  Gio Buonaccorsi altered function gfit_measure_image to take mask as 
 34  *  parameter and moved stack manipulation to gfit_pixels_proc.
 35  *  Also altered function gamma_fit_region to take mask as parameter
 36  *  and moved stack manipulation to pixel_gfit_region.
 37  *  Done to improve lib / tool separation 19 Feb 2003.
 38  *
 39  *********
 40 */
 41 
 42 
 43 #if HAVE_CONFIG_H
 44 #   include <config.h>
 45 #endif
 46 
 47 #include <stdio.h>
 48 #include <math.h>
 49 #include <float.h>
 50 #include <tina/sys/sysDef.h>
 51 #include <tina/sys/sysPro.h>
 52 #include <tina/math/mathDef.h>
 53 #include <tina/math/mathPro.h>
 54 #include <tina/image/imgDef.h>
 55 #include <tina/image/imgPro.h>
 56 #include <tina/medical/med_StimDef.h>
 57 #include <tina/medical/medStim_alloc.h>
 58 #include <tina/medical/medStim_tdata.h>
 59 #include <tina/medical/medStim_gamma.h>  /* HR */
 60 
 61 
 62 #define NP 3
 63 #define MIN_BOLUS_WIDTH NP+1
 64 #define SEC_PEAK  2.0
 65 #define T0SCAN  10 
 66 #define TTPSCAN  15
 67 #define RR_INIT  3.6
 68 #define MIN_TTP_FRAC  0.5
 69 #define TTP_FRAC_STEP 0.1
 70 
 71 #define FTOL 1.0e-6
 72 
 73 static double SCL = 10000.0;                            /* static data! */
 74 static double MIN_T0 = 18.0;                            /* static data! */
 75 static double AVE_MTT = 7.0;                            /* static data! */
 76 static double RECIRC_CUT = 0.2;                         /* static data! */
 77 static double RECIRC_PERIOD = 10.0;                     /* static data! */
 78 static Imrect *perf_TTP=NULL;                           /* static data! */
 79 static Imrect *perf_CBV=NULL;                           /* static data! */
 80 static Imrect *perf_MTT=NULL;                           /* static data! */
 81 static Imrect *perf_ERR=NULL;                           /* static data! */
 82 static Imrect *perf_RCC=NULL;                           /* static data! */
 83 
 84 static Imrect *diff_S0 =NULL;                           /* static data! */
 85 static Imrect *diff_ADC=NULL;                           /* static data! */
 86 
 87 static float MAX_TTP;                                   /* static data! */
 88 static Perfusion *perf_global = NULL;                   /* static data! */
 89 
 90 static Vec2 mask_centre;                             /* HR */
 91 static float adc_centre, s0_centre;                  /* HR */
 92 static float b0=0.0, b1=100.0, b2=500.0, b3=900;     /* HR: default */
 93 static int no_data=0, init_flag=0;                   /* HR: 22/01/2013 */
 94 static float *b_vals=NULL;                           /* HR: 22/01/2013 */
 95 
 96 
 97 double *scl_get()
 98 {
 99    return(&SCL);
100 }
101 
102 double *min_t0_get()
103 {
104    return(&MIN_T0);
105 }
106 
107 double *ave_mtt_get()
108 {
109    return(&AVE_MTT);
110 }
111 
112 double *recirc_cut_get()
113 {
114    return(&RECIRC_CUT);
115 }
116 
117 double *recirc_period_get()
118 {
119    return(&RECIRC_PERIOD);
120 }
121 
122 static float ave_s0_pre_bolus(Perfusion *p)
123 {
124   float *st_1d = p->st_1d;
125   int    k , bolus_time = p->bolus_time;
126   float  noisesum = 0.0;
127 
128   if (bolus_time-2 <= p->n1) bolus_time = p->n1+3;
129 
130   for (k = p->n1; k < bolus_time-2 && k < p->n2; k++)
131     noisesum += st_1d[k];
132 
133 
134   return (noisesum/(bolus_time-2-p->n1));
135 }
136 
137 static float conv_r2_to_st(float r2, float s0, float te)
138 {
139     float st;
140     st = s0*exp(-r2*te/SCL);
141     return (st);
142 }
143 
144 static float *conv_st_to_r2(Perfusion *p)
145 {
146   float  *st_1d = p->st_1d;
147   int     n1 = p->n1;
148   int     n2 = p->n2;
149   float   s0 = p->s0_av;
150   float   te = p->te;
151   float  *r2, val;
152   int      k;
153 
154   if ((r2 = fvector_alloc(n1, n2)) == NULL)
155     {
156       error("conv_st_to_r2: memory alloc error ", warning);
157       return(NULL);
158     }
159 
160     for (k = n1; k < n2; k++)
161     {
162       val = SCL*(-log(st_1d[k]/s0))/te;
163       if (isnan(val))
164         {
165           r2[k] = 0.0;
166                                         printf("%f\n", val);
167         } 
168       else
169         r2[k] = val;
170       val = conv_r2_to_st(val, s0, te);
171     }
172 
173   return(r2);
174 }
175 
176 static double gamma_func( double ka0, double ka1, double ka2,
177                    float t, float termi)
178 {
179    double t1 = t - 0.5 *termi ,t2 = t+0.5*termi;
180    double con=0.0, weight;
181 
182    if (t1<=0.0)
183    {
184       if (t2<=0.0)
185       {
186          con = 0.0;
187       }
188       else
189       {
190          weight = t2/termi;
191          con = 0.5*weight*ka0*pow(t2,ka1)*exp(-t2/ka2);
192       }
193    }
194    else
195    {
196       con += 0.5*ka0*pow(t1,ka1)*exp(-t1/ka2);
197       con += 0.5*ka0*pow(t2,ka1)*exp(-t2/ka2);
198    }
199    return(con);
200 }
201 
202 static double cut_off(double ka1, double ka2, float cut)
203 {
204    double ans1=0, ans2=50;
205 
206    while (fabs(ans1 - ans2) > 0.001)
207    {
208        ans1 = ans2;
209        ans2 = ka2*(ka1*log(ans1)-log(cut));
210    }
211    if (ans1>0.0)
212       return(ans1);
213    else
214       return(0.0);
215 }
216 
217 static void r2_ther_stats(Perfusion *perf_data, int offset,
218                           float *norm1, float *norm2, float *dot, double *chi)
219 {
220    double ka0,ka1,ka2, t0;
221    int n1 = perf_data->n1;
222    int n2 = perf_data->n2;
223    float termi = perf_data->tscale;
224    float *data = perf_data->r2;
225    float *r2_ther = perf_data->r2_ther;
226    float t;
227    int i, icut, tcut;
228 
229    ka1 = perf_data->rr;  /* b  */
230    ka2 = perf_data->bb;  /* r  */
231    ka0 = perf_data->qq;  /* q  */
232    t0 = perf_data->bolus_time*termi-perf_data->tt;
233    tcut = perf_data->bolus_time;
234    icut = perf_data->bolus_end;
235 
236    if (offset == 0)
237    {
238       if (perf_data->r2_ther == NULL)
239            perf_data->r2_ther = fvector_alloc(n1, n2);
240       r2_ther = perf_data->r2_ther;
241       *norm1 = 0.0;
242       for (i=tcut;
243            i<icut&&i<n2;
244            i++)
245       {
246          t = termi*i-t0;
247          r2_ther[i] = gamma_func(ka0,ka1,ka2,t,termi);
248          *norm1 += r2_ther[i]*r2_ther[i];
249       }
250       *norm1 = sqrt(*norm1);
251    }
252    *dot = 0.0;
253    *norm2 = 0.0;
254    for (i=tcut;
255         i<icut&&i+offset<n2;
256         i++)
257    {
258       *dot += r2_ther[i]*data[i+offset];
259       *norm2 += data[i+offset]*data[i+offset];
260    }
261    *chi = (*norm2)*(*dot/(*norm1*sqrt(*norm2)) -1.0)/(float)(i-NP-tcut);
262 }
263 
264 static double gamma_chi2(int ndim, double *x)
265 /* minimise errors on measured data */
266 {
267   float termi = perf_global->tscale;
268   double ka1,ka2, tt;
269   double chi;
270   float dot,norm1,norm2;
271 
272   ka1=x[0]; ka2=x[1]; tt=x[2];
273 
274   if (ka1 > 0.0 && ka2 > 1.0e-6 && ka1*ka2 < MAX_TTP
275      && tt < 2.0*termi && tt > -2.0*termi );
276     /* protect against NaN and Infty etc */
277   else
278   {
279     perf_global->qq = 0.0;
280     return(FLT_MAX);
281   }
282 
283   perf_global->qq = 1.0;
284   perf_global->bb = ka2;
285   perf_global->rr = ka1;
286   perf_global->tt = x[2];
287   r2_ther_stats(perf_global,0,&norm1,&norm2,&dot,&chi); 
288   if (dot>1.0 && norm1 > 0.0)
289   {
290     perf_global->qq *= dot/(norm1*norm1);
291     return(-chi);
292   }
293   else
294   {
295     perf_global->qq = 0.0;
296     return(FLT_MAX);
297   }
298 }
299 
300 static void gamma_simplex(Perfusion *p)
301 {
302   double   *x;
303   double   *lambda;
304   int ndim = NP;
305   
306   perf_global = p; /* access to parameters for gamma_chi() */
307   x = (double *)ralloc(ndim*sizeof(double));
308   x[0] = p->rr;
309   x[1] = p->bb;
310   x[2] = p->tt;
311 
312   lambda = (double *)ralloc(ndim*sizeof(double));
313   lambda[0] = 0.1;
314   lambda[1] = 0.1;
315   lambda[2] = 2.0;
316   simplexmin2(ndim, x, lambda, gamma_chi2, NULL, FTOL, /*(void (*) ()) format */ NULL);
317   p->rr=x[0];
318   p->bb=x[1];
319   p->tt=x[2];
320   gamma_chi2(ndim,x);
321   rfree(x);
322   rfree(lambda);
323 }
324 
325 float mono_exp_func(double s0, float b, double adc)              /* HR: new */
326 {
327     float exp_func;
328 
329     exp_func = s0 * exp(-b * adc);
330 
331     return(exp_func);
332 }
333 
334 /* minimise errors on measured data */
335 static double diff_chi2(int ndim, double *x, void *data)             /* HR: modified: 28/01/2013 */
336 {
337     int    i;
338     float  fit_func;
339     double chisq;
340 
341     Diffusion *dif_data;
342 
343     if ((dif_data = (Diffusion *)data) == NULL)
344         return; /*error code*/
345 
346     chisq = 0.0;
347 
348     for(i=0; i<no_data; i++)
349     {
350         fit_func = mono_exp_func(x[0], dif_data->b_values[i], x[1]);
351 
352         chisq += (fit_func - dif_data->signals[i]) * (fit_func - dif_data->signals[i]) / ((float) (no_data-1));
353     }
354 
355     return(chisq);
356 }
357 
358 static void init_params_s0adc(double *paramts, int bx, int by, float sx, float sy) /* HR: new: 28/01/2013 */
359 {
360     double so_param, adc_param;
361 
362     if ( (sx == 0) || (sy == 0) || (bx == by) )
363     {
364         adc_param = 0.0;
365     }
366     else
367     {
368         adc_param = (double) ((1.0 / (by-bx)) * log(sx/sy));
369     }
370 
371     if (bx == 0)
372     {
373         so_param = (double) sx;
374     }
375     else
376     {
377         so_param = (double) (sx * exp(bx * adc_param));
378     }
379 
380     paramts[0] = so_param;
381     paramts[1] = adc_param;
382 }
383 
384 /* inputs:  a number of values (typically 4) of diffusion for one cylinder */
385 /* outputs: 2 parameter estimates for the mono exponential fit */
386 static void diff_simplex(Diffusion *diff_data)                       /* HR: modified: 28/01/2013 */
387 {
388     double *lambda, *params, chisq_dist;
389     int     no_prms = 2;
390 
391     params = (double *)ralloc(no_prms*sizeof(double));
392     lambda = (double *)ralloc(no_prms*sizeof(double));
393 
394     lambda[0] = 0.5;
395     lambda[1] = 0.1;
396 
397     init_params_s0adc(params, diff_data->b_values[0], diff_data->b_values[1], diff_data->signals[0], diff_data->signals[1]);
398 
399     if( (params[0]*params[1]) != 0 )
400     {
401         chisq_dist = simplexmin2(no_prms, params, lambda, diff_chi2, diff_data, FTOL, /*(void (*) ()) format*/ NULL);
402     }
403     else
404     {
405         params[0] = 0.0;
406         params[1] = 0.0;                                             /* do not run simplex for this pixel */
407     }
408 
409     diff_data->fit_params = params;
410 
411     rfree(lambda);
412 }
413 
414 static float  emap(Perfusion *p)
415 {
416     int i, icut, tcut;
417     float cut, dof, t;
418     double chi, con, chitot=0.0;
419     double ka0,  ka1,  ka2, t0;
420     float termi = p->tscale;
421     float *data = p->r2;
422 
423     ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
424     cut = cut_off(ka1,ka2,RECIRC_CUT*gamma_func(1.0,ka1,ka2,ka1*ka2,termi)); 
425     icut = tina_int(cut/termi); 
426     tcut = tina_int(t0/termi);
427     for (i=tcut,dof=-4.0;
428        i<(tcut+icut)&&i<p->n2;
429        dof+=1.0,i++)
430    {
431       t = termi*i-t0;
432       con= gamma_func(ka0,ka1,ka2,t,termi);
433       chi= data[i] - con;
434       chitot += chi*chi;
435    }
436    if (dof>0.0)
437       return(chitot/dof);
438    else
439       return(0.0);
440 }
441 
442 static float  rcbv(Perfusion *p)
443 {
444    int i, tcut;
445    float ttot=0.0, t;
446    double ka0, ka1, ka2, t0;
447    float termi = p->tscale;
448    double con;
449     
450    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
451    tcut = tina_int(t0/termi);
452    for (i=tcut; i<p->n2; i++)
453    {
454       t = termi*i-t0;
455       con= gamma_func(ka0,ka1,ka2,t,termi);
456       ttot += con;
457    }
458    return(ttot);
459 }
460 
461 static float  recirc(Perfusion *p)
462 {
463    int i, icut, tcut, rcut;
464    float cut, ttot=0.0, t;
465    double ka0, ka1, ka2, t0;
466    float termi = p->tscale;
467    float *data = p->r2;
468    double con;
469    
470     ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
471     cut = cut_off(ka1,ka2,RECIRC_CUT*gamma_func(1.0,ka1,ka2,ka1*ka2,termi)); 
472     icut = tina_int(cut/termi); 
473     tcut = tina_int(t0/termi);
474     rcut = tina_int(RECIRC_PERIOD/termi);
475     for (i=tcut+icut; i<tcut+icut+rcut && i<p->n2; i++)
476     {
477        t = termi*i-t0;
478        con= gamma_func(ka0,ka1,ka2,t,termi);
479        ttot += data[i] - con;
480     }
481     return(ttot);
482 }
483 
484 static float  bolus_m0(Perfusion *p)
485 {
486    int i, tcut;
487    float t;
488    double ka0, ka1, ka2, t0;
489    float termi = p->tscale;
490    double con;
491    double mean=0.0;
492     
493    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
494    tcut = tina_int(t0/termi);
495    for (i=tcut; i<p->n2; i++)
496    {
497       t = termi*i-t0;
498       con= gamma_func(ka0,ka1,ka2,t,termi);
499       mean += con*(i+0.5)*termi;
500    }
501    
502    return(mean);
503 }
504 
505 static float  bolus_m1(Perfusion *p)
506 {
507    int i, tcut;
508    float t;
509    double ka0, ka1, ka2, t0;
510    float termi = p->tscale;
511    double con;
512    double mean=0.0;
513 
514    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
515    tcut = tina_int(t0/termi);
516    for (i=tcut; i<p->n2; i++)
517    {
518       t = termi*i-t0;
519       con= gamma_func(ka0,ka1,ka2,t,termi);
520       mean += con*termi*(i+0.5)*termi*(i+0.5);
521    }
522    return(mean);
523 }
524 
525 static float   *recon_ther_r2(float termi, Perfusion *p)
526 {
527   int     i;
528   float  *r2_ther, td , t0;
529 
530   if ((r2_ther = fvector_alloc(p->n1, p->n2)) == NULL)
531     {
532       error("recon_ther_r2: memory alloc error", warning);
533       return(NULL);
534     } 
535 
536   t0 = p->bolus_time*termi - p->tt;
537   for (i = p->n1; i < p->n2; i++) 
538   {
539     td = termi*i -t0;
540     r2_ther[i] = gamma_func(p->qq,p->rr,p->bb,td,termi);
541   }
542 
543   return(r2_ther);
544 }
545 
546 static void search_chi_fits(Perfusion *perf_data, float **chi2array,float **dotarray,
547                             float *norm1)
548 {
549    float kmax, jmax=0, dotmax, norm1_max=0;
550    int j,k,test;
551    float chimax;
552 
553    dotmax = -FLT_MAX;
554    chimax = -FLT_MAX;
555    kmax = T0SCAN-1;
556    for (j=0;j<TTPSCAN;j++)
557    {
558       for (k=0;k<T0SCAN;k++)
559       {
560          test = 0;
561          if (k>0                && chi2array[j][k]<chi2array[j][k-1]  )  continue; 
562          if (k>0 && j >0        && chi2array[j][k]<chi2array[j-1][k-1])  continue; 
563          if (k>0 && j<TTPSCAN-1 && chi2array[j][k]<chi2array[j+1][k-1])  continue;
564          if (k<T0SCAN-1        && chi2array[j][k]<chi2array[j][k+1]  )  continue;
565          if (k<T0SCAN-1 && j>0 && chi2array[j][k]<chi2array[j-1][k+1])  continue;
566          if (k<T0SCAN-1 && j<TTPSCAN-1 && chi2array[j][k]<chi2array[j+1][k+1])  continue;
567          if (j>0                && chi2array[j][k]<chi2array[j-1][k]  )  continue;
568          if (j<TTPSCAN-1        && chi2array[j][k]<chi2array[j+1][k]  )  continue;
569 
570          if ( ( dotarray[j][k] > dotmax && k < kmax) ||
571               ( dotarray[j][k] > SEC_PEAK*dotmax && k > kmax))
572          {
573             chimax = chi2array[j][k];
574             dotmax = dotarray[j][k];
575             norm1_max = norm1[j];
576             kmax = k;
577             jmax = j;
578          }
579       }
580    }
581    if (dotmax > 0)
582    {
583       float t0,cut;
584       int icut,tcut;
585       perf_data->qq *= dotmax/norm1_max;
586       perf_data->tt = -(kmax)*perf_data->tscale;
587       perf_data->bb = (MIN_TTP_FRAC+(TTP_FRAC_STEP*jmax))*AVE_MTT/perf_data->rr;
588       t0 = perf_data->tscale*perf_data->bolus_time-perf_data->tt;
589 
590       cut = cut_off(perf_data->rr,perf_data->bb,
591             RECIRC_CUT*gamma_func(1.0,perf_data->rr,perf_data->bb,
592             perf_data->rr*perf_data->bb,perf_data->tscale));
593 
594       icut = tina_int(cut/perf_data->tscale);
595       tcut = tina_int(t0/perf_data->tscale);
596       perf_data->bolus_time -= tina_int(perf_data->tt/perf_data->tscale);
597       perf_data->tt -= tina_int(perf_data->tt/perf_data->tscale)*perf_data->tscale;
598       perf_data->refit = 1;
599       perf_data->bolus_end = tcut+icut;
600       if (perf_data->bolus_end - perf_data->bolus_time < MIN_BOLUS_WIDTH)
601       {
602          perf_data->bolus_end = perf_data->bolus_time+MIN_BOLUS_WIDTH;
603 /*
604          perf_data->qq *= 0.0;
605          perf_data->tt = -T0SCAN*perf_data->tscale/2.0;
606          perf_data->bb = AVE_MTT/perf_data->rr;
607 */
608       }
609    }
610    else
611    {
612       perf_data->qq *= 0.0;
613       perf_data->tt = -T0SCAN*perf_data->tscale/2.0;
614       perf_data->bb = AVE_MTT/perf_data->rr;
615    }
616 }
617 
618 static void gamma_match(Perfusion *perf_data)
619 {
620    int j,k;
621    double ka1,ka2;
622    float termi = perf_data->tscale;
623    float norm2;
624    float *norm1;
625    float dot;
626    double chi;
627    float **chi2array;
628    float **dotarray;
629    double cut;
630 
631    ka1 = perf_data->rr = RR_INIT;  /* b  */
632    perf_data->bb = T0SCAN/RR_INIT;  /* r  */
633    perf_data->qq = 1.0;  /* q  */
634    perf_data->tt = 0.0;
635    norm1 = fvector_alloc(0,TTPSCAN);
636    dotarray = farray_alloc(0,0,TTPSCAN,T0SCAN);
637    chi2array = farray_alloc(0,0,TTPSCAN,T0SCAN);
638 
639    for (ka2 = MIN_TTP_FRAC*AVE_MTT/ka1, j=0; j<TTPSCAN ;
640                        ka2+=TTP_FRAC_STEP*AVE_MTT/ka1,j++)
641    {
642       cut = cut_off(ka1,ka2,RECIRC_CUT*gamma_func(1.0,ka1,ka2,ka1*ka2,termi));
643       perf_data->bolus_end = perf_data->bolus_time + tina_int(cut/termi);
644 
645       perf_data->bb = ka2;
646       for (k=0;k<T0SCAN;k++)
647       {
648          r2_ther_stats(perf_data,k,&norm1[j],&norm2,&dot,&chi);
649          dotarray[j][k] = dot/(norm1[j]);
650          chi2array[j][k] = chi;
651       }
652    }
653    search_chi_fits(perf_data, chi2array,dotarray,norm1);
654    farray_free(chi2array,0,0,TTPSCAN,T0SCAN);
655    farray_free(dotarray,0,0,TTPSCAN,T0SCAN);
656    fvector_free(norm1,0);
657 }
658 
659 int   gamma_fit_pixel(Perfusion *perf_data, Vec2 v, void ***imptrs)
660 {
661   Sequence *seq= (Sequence *)seq_get_current();
662   List   *store = (List *)get_current_seq_start_el();
663   Imrect   *im = NULL;
664   float    *te=NULL;
665   float our_time;
666         
667   if ((store == NULL) || (store->to == NULL))
668     return(-1);
669   im = (Imrect *)(store->to);
670 
671   if ((our_time = seq->dim[3]) == 1.0)
672     {
673       error("gamma_fit: timing info missing", warning);
674       return(-1);
675     }
676 
677   perf_data->tscale = our_time;
678   MAX_TTP = AVE_MTT*(TTPSCAN*TTP_FRAC_STEP + MIN_TTP_FRAC);
679   if ((te = (float *)prop_get(seq->props, TE_DATA)) == NULL)
680     {
681       error("gamma_fit: echo time not found", warning);
682       return(-1);
683     }
684 
685   if (perf_data->st_1d!=NULL) fvector_free(perf_data->st_1d, perf_data->n1);
686   if (perf_data->r2!=NULL) fvector_free(perf_data->r2, perf_data->n1);
687   if (perf_data->r2_ther!=NULL) fvector_free(perf_data->r2_ther, perf_data->n1);
688   perf_data->n1 = seq->offset;
689   perf_data->n2 = get_end_frame(seq);
690 
691   /*mjs 8/11/05 assumes that all TE's in sequence are same as that of first image in sequence */
692   perf_data->te = te[seq->offset];
693   perf_data->st_1d = s_2d_pixel_get(imptrs, v, perf_data);
694   perf_data->bolus_time = MIN_T0/perf_data->tscale;
695   perf_data->s0_av = ave_s0_pre_bolus(perf_data);
696   perf_data->r2 = conv_st_to_r2(perf_data);
697   gamma_match(perf_data);
698   if (perf_data->qq > 0.0 ) gamma_simplex(perf_data);
699   perf_data->r2_ther = recon_ther_r2(perf_data->tscale, perf_data);
700 
701   return (0);
702 }
703 
704 void    set_diffusion_bvalues(float *bvals_vec, int no_bvals)        /* HR: new: 22/01/2013 */
705 {
706     int i;
707 
708     no_data = no_bvals;
709 
710     if ((b_vals = fvector_alloc(0, no_data)) == NULL)
711     {
712         error("set_diffusion_bvalues: memory alloc error ", warning);
713         return;
714     }
715 
716     for(i=0; i<no_data; i++)
717     {
718         b_vals[i] = bvals_vec[i];
719     }
720 }
721 
722 static void init_dif_struc(Diffusion *d)                             /* HR: modified: 22/01/2013 */
723 {
724     int i;
725     float *b_data=NULL;
726 
727 /* HR: for number of b values other than 4, input b values for the number of slices loaded, etc.
728    dif_data->no_b = get_end_frame(seq) -seq->offset +1;
729 */
730 
731     if ( ( no_data == 0 ) || (init_flag == 1) )
732     {
733         init_flag = 1;
734         no_data = 4;
735 
736         if ((b_data = fvector_alloc(0, no_data)) == NULL)
737         {
738             error("init_dif_struc: memory alloc error ", warning);
739             return;
740         }
741 
742         b_data[0] = b0;
743         b_data[1] = b1;
744         b_data[2] = b2;
745         b_data[3] = b3;
746     }
747     else
748     {
749         if ((b_data = fvector_alloc(0, no_data)) == NULL)
750         {
751             error("init_dif_struc: memory alloc error ", warning);
752             return;
753         }
754 
755         for(i=0; i<no_data; i++)
756         {
757             b_data[i] = b_vals[i];
758         }
759     }
760 
761     d->b_values = b_data;
762     d->no_b = no_data;
763 }
764 
765 static void print_diff_data(Diffusion *d)                            /* HR: new */  
766 {
767     printf("number of image slices : %d \n", d->no_b);
768     printf("b values : %f %f %f %f \n", d->b_values[0],d->b_values[1],d->b_values[2],d->b_values[3]);
769     printf("pixel signals : %f %f %f %f \n", d->signals[0],d->signals[1],d->signals[2],d->signals[3]);
770 }
771 
772 
773 int   diff_fit_pixel(Diffusion *dif_data, Vec2 v, void ***imptrs)    /* HR: modified: 22/01/2013 */
774 {
775   Sequence *seq= (Sequence *)seq_get_current();
776   List   *store = (List *)get_current_seq_start_el();
777   Imrect   *im = NULL;
778 
779   if ((store == NULL) || (store->to == NULL))
780     return(-1);
781   im = (Imrect *)(store->to);
782 
783   if (dif_data->b_values!=NULL) fvector_free(dif_data->b_values, 0);
784   if (dif_data->signals!=NULL)  fvector_free(dif_data->signals, 0);
785   if (dif_data->fit_params!=NULL) dvector_free(dif_data->fit_params, 0);
786 
787   init_dif_struc(dif_data);
788 
789   dif_data->signals = diff_2d_pixel_get(imptrs, v, dif_data);
790 /*
791   print_diff_data(dif_data);
792 */
793   diff_simplex(dif_data);
794 
795   return (0);
796 }
797 
798 /*
799  *   Altered to improve lib / tool separation - GAB 19 Feb 2003
800  *   Stack manipulation moved to pixel_gfit_region, mask passed as parameter.
801  */
802 double  gamma_fit_region(Perfusion *perf_data, Imrect *mask)
803 {
804   
805   Sequence *seq= (Sequence *)seq_get_current();
806   List     *store = (List *)get_current_seq_start_el();
807   Imrect   *tmask = NULL, *fmask = NULL, *im;
808   float    *te=NULL;
809   int       i, j;
810   void   ***imptrs=NULL;
811   double    t0;
812   float our_time;       
813 
814   if ((store == NULL) || (store->to == NULL))
815     return(0.0);
816   im = (Imrect *)(store->to);
817 
818   if ((our_time = seq->dim[3]) == 1.0)
819     {
820       error("gamma_fit: timing info missing", warning);
821       return(-1);
822     }
823 
824   perf_data->tscale = our_time;
825   MAX_TTP = AVE_MTT*(TTPSCAN*TTP_FRAC_STEP + MIN_TTP_FRAC);
826 
827   if ((te = (float *)prop_get(seq->props, TE_DATA)) == NULL)
828     {
829       error("gamma_fit: echo time not found", warning);
830       return(0.0);
831     }
832 
833   if ((tmask = im_alloc(mask->height, mask->width, NULL, char_v)) == NULL)
834     {
835       error("gamma_fit_region: memory alloc error 2", warning);
836       return(0.0);
837     }
838 
839   for (i = mask->region->ly; i < mask->region->uy; i++)
840     for (j = mask->region->lx; j < mask->region->ux; j++)
841       IM_CHAR(tmask, i, j) = 1;
842   fmask = im_prod(tmask, mask);
843 
844   if (perf_data->st_1d!=NULL) fvector_free(perf_data->st_1d, perf_data->n1);
845   if (perf_data->r2!=NULL) fvector_free(perf_data->r2, perf_data->n1);
846   if (perf_data->r2_ther !=NULL) fvector_free(perf_data->r2_ther, perf_data->n1);
847   perf_data->n1 = seq->offset;
848   perf_data->n2 = get_end_frame(seq);
849 
850   imptrs = seq_slice_init(seq);
851   /* mjs modified 8/11/05 to access TE data properly, assuming all image TEs equal first in sequence */
852   perf_data->te = te[seq->offset];
853   perf_data->st_1d = s_2d_mask_mean(imptrs, fmask, perf_data);
854   perf_data->bolus_time = MIN_T0/perf_data->tscale;
855   perf_data->s0_av = ave_s0_pre_bolus(perf_data);
856   perf_data->r2 = conv_st_to_r2(perf_data); 
857   perf_data->tscale = our_time;
858   gamma_match(perf_data);
859   if (perf_data->qq > 0.0 ) gamma_simplex(perf_data);
860   perf_data->r2_ther = recon_ther_r2(perf_data->tscale, perf_data); 
861 
862   im_free(fmask);
863   im_free(tmask);
864   
865   t0 = perf_data->bolus_time*perf_data->tscale-perf_data->tt;
866   return (t0);
867 }
868 
869 /*
870  *   Altered to improve lib / tool separation - GAB 19 Feb 2003
871  *   Stack manipulation moved to gfit_pixels_proc, mask passed as parameter.
872  */
873 void    gfit_measure_image(double err_thresh, Imrect *mask)
874 {
875   Perfusion *p = NULL;
876   Sequence  *seq= (Sequence *)seq_get_current();
877   List      *store = (List *)get_current_seq_start_el();
878   Imrect    *im = NULL;
879   Imregion  *roi;
880   float      cbv, rcc, mtt, err, ttp;
881   Vec2       pos;
882   int        i, j;
883   void    ***imptrs=NULL;
884 
885   im = (Imrect *)(store->to);
886   roi = im->region;
887 
888   if (perf_TTP!=NULL) im_free(perf_TTP);
889   perf_TTP = im_alloc(im->height, im->width, roi, float_v);
890   if (perf_CBV!=NULL) im_free(perf_CBV);
891   perf_CBV = im_alloc(im->height, im->width, roi, float_v);
892   if (perf_MTT!=NULL) im_free(perf_MTT);
893   perf_MTT = im_alloc(im->height, im->width, roi, float_v);
894   if (perf_ERR!=NULL) im_free(perf_ERR);
895   perf_ERR = im_alloc(im->height, im->width, roi, float_v);
896   if (perf_RCC!=NULL) im_free(perf_RCC);
897   perf_RCC = im_alloc(im->height, im->width, roi, float_v);
898 
899 /*   Stack manipulation was here - GAB 19 Feb 2003   */
900   if (mask != NULL)
901     roi = mask->region;
902 
903   imptrs = seq_slice_init(seq);
904 
905   p = perfusion_alloc();
906   for (i = roi->ly; i < roi->uy; i++)
907   {
908     p->refit = 0;
909     for (j = roi->lx; j < roi->ux; j++)
910     {
911        if (mask == NULL || IM_CHAR(mask, i, j) > 0)
912        {
913           vec2_y(pos) = i + 0.5;
914           vec2_x(pos) = j + 0.5;
915           if (gamma_fit_pixel(p, pos,imptrs) != -1)
916           {
917              err = emap(p);
918              cbv = rcbv(p);
919 /*
920              if (cbv > 0.0 && sqrt(err/cbv) < err_thresh)
921                 p->refit = 1;
922              else
923              {
924                 p->refit = 0;
925                 gamma_fit_pixel(p, pos,imptrs);
926              }
927 */
928              err = emap(p);
929              rcc = recirc(p);
930              cbv = rcbv(p);
931              if (cbv > 0.0 && sqrt(err/cbv) < err_thresh)
932                 p->refit = 1;
933 
934              if(cbv > 0.0)
935              {
936                 ttp = bolus_m0(p)/cbv;
937                 mtt = sqrt(ttp*ttp - 2.0*ttp*bolus_m0(p)/cbv + bolus_m1(p)/cbv);
938              }
939              else
940              {
941                 ttp = 0.0;
942                 mtt = 0.0;
943              }
944              IM_PIX_SET(perf_CBV, i, j, cbv);
945              IM_PIX_SET(perf_TTP, i, j, ttp - p->n1*p->tscale);
946              IM_PIX_SET(perf_MTT, i, j, mtt);
947              IM_PIX_SET(perf_ERR, i, j, err);
948              IM_PIX_SET(perf_RCC, i, j, rcc);
949           }
950        }
951      }
952      if (i%10 == 0) format("raster %d done \n",i);
953    }
954    perfusion_free(p);
955 
956   return;
957 }
958 
959 
960 void    dfit_measure_image(double err_thresh, Imrect *mask)          /* HR: new */
961 {
962   Diffusion *d = NULL;
963   Sequence  *seq= (Sequence *)seq_get_current();
964   List      *store = (List *)get_current_seq_start_el();
965   Imrect    *im = NULL;
966   Imregion  *roi;
967   float      s0, adc, s0_norm=1.0, adc_norm=100000.0, std_dev, adc_range=128.0, S0_range=2048.0; /* adc_range depends on the liquids used in the five cylinders of the phantom */
968   Vec2       pos;
969   int        i, j, hnum, xmask_centre, ymask_centre;
970   void    ***imptrs=NULL;
971   shistogram **hists;
972 
973   hists = (shistogram **)hist_vec();
974 
975   if( hists[150] !=NULL) hfree(hists[150]);
976   hists[150] = hbook1(150, "adc (8 bins)\0", 0.0, adc_range, 8);          /* course-to-fine aproach for adc */
977   if( hists[151] !=NULL) hfree(hists[151]);
978   hists[151] = hbook1(151, "adc (16 bins)\0", 0.0, adc_range, 16);
979   if( hists[152] !=NULL) hfree(hists[152]);
980   hists[152] = hbook1(152, "adc (32 bins)\0", 0.0, adc_range, 32);
981   if( hists[153] !=NULL) hfree(hists[153]);
982   hists[153] = hbook1(153, "adc (64 bins)\0", 0.0, adc_range, 64);
983   if( hists[154] !=NULL) hfree(hists[154]);
984   hists[154] = hbook1(154, "adc (128 bins)\0", 0.0, adc_range, 128);
985   if( hists[155] !=NULL) hfree(hists[155]);
986   hists[155] = hbook1(155, "adc (256 bins)\0", 0.0, adc_range, 256);
987 
988   if( hists[170] !=NULL) hfree(hists[170]);
989   hists[170] = hbook1(170, "S0 (8 bins)\0", 0.0, S0_range, 8);          /* course-to-fine aproach for S0 */
990   if( hists[171] !=NULL) hfree(hists[171]);
991   hists[171] = hbook1(171, "S0 (16 bins)\0", 0.0, S0_range, 16);
992   if( hists[172] !=NULL) hfree(hists[172]);
993   hists[172] = hbook1(172, "S0 (32 bins)\0", 0.0, S0_range, 32);
994   if( hists[173] !=NULL) hfree(hists[173]);
995   hists[173] = hbook1(173, "S0 (64 bins)\0", 0.0, S0_range, 64);
996   if( hists[174] !=NULL) hfree(hists[174]);
997   hists[174] = hbook1(174, "S0 (128 bins)\0", 0.0, S0_range, 128);
998   if( hists[175] !=NULL) hfree(hists[175]);
999   hists[175] = hbook1(175, "S0 (256 bins)\0", 0.0, S0_range, 256);
1000 
1001   im = (Imrect *)(store->to);
1002   roi = im->region;
1003 
1004   if (diff_S0!=NULL) im_free(diff_S0);
1005   diff_S0 = im_alloc(im->height, im->width, roi, float_v);
1006   if (diff_ADC!=NULL) im_free(diff_ADC);
1007   diff_ADC = im_alloc(im->height, im->width, roi, float_v);
1008 
1009   ymask_centre = (int) vec2_y(mask_centre);
1010   xmask_centre = (int) vec2_x(mask_centre);
1011 
1012   if (mask != NULL)
1013     roi = mask->region;
1014 
1015   imptrs = seq_slice_init(seq);
1016 
1017   d = diffusion_alloc();
1018   for (i = roi->ly; i < roi->uy; i++)
1019   {
1020     for (j = roi->lx; j < roi->ux; j++)
1021     {
1022        if (mask == NULL || IM_FLOAT(mask, i, j) > 0.5)
1023        {
1024           vec2_y(pos) = i + 0.5;
1025           vec2_x(pos) = j + 0.5;
1026 
1027           if (diff_fit_pixel(d, pos,imptrs) != -1)
1028           {
1029              s0  = s0_norm  * d->fit_params[0];
1030              adc = adc_norm * d->fit_params[1];
1031 
1032              IM_PIX_SET(diff_S0, i, j, s0);
1033              IM_PIX_SET(diff_ADC, i, j, adc);
1034 
1035              hfill1s(hists[150], adc , 1.0);
1036              hfill1s(hists[151], adc , 1.0);
1037              hfill1s(hists[152], adc , 1.0);
1038              hfill1s(hists[153], adc , 1.0);
1039              hfill1s(hists[154], adc , 1.0);
1040              hfill1s(hists[155], adc , 1.0);
1041 
1042              hfill1s(hists[170], s0 , 1.0);
1043              hfill1s(hists[171], s0 , 1.0);
1044              hfill1s(hists[172], s0 , 1.0);
1045              hfill1s(hists[173], s0 , 1.0);
1046              hfill1s(hists[174], s0 , 1.0);
1047              hfill1s(hists[175], s0 , 1.0);
1048 
1049              if( (j==xmask_centre) && (i==ymask_centre) )
1050              {
1051                  adc_centre = adc;
1052                  s0_centre  = s0;
1053              }
1054           }
1055        }
1056      }
1057 /*
1058      if (i%10 == 0) format("raster %d done \n",i);
1059 */
1060    }
1061 
1062    diffusion_free(d);
1063 }
1064 
1065 void set_mask_centre(Vec2 *msk_centre)           /* HR: new: set the centre coords to a static global variable */
1066 {
1067     vec2_y(mask_centre) = vec2_y(*msk_centre);
1068     vec2_x(mask_centre) = vec2_x(*msk_centre);
1069 }
1070 
1071 float get_mask_centre_adc()                      /* HR: new: get the adc value from static global variable */
1072 {
1073     return(adc_centre);
1074 }
1075 
1076 float get_mask_centre_s0()                      /* HR: new: get the s0 value from static global variable */
1077 {
1078     return(s0_centre);
1079 }
1080 
1081 Imrect *gfit_get_image(int perf_type)
1082 {
1083   if (perf_type ==0)
1084       return(im_copy(perf_TTP));
1085   else if (perf_type ==1)
1086       return(im_copy(perf_CBV));
1087   else if (perf_type ==2)
1088       return(im_copy(perf_MTT));
1089   else if (perf_type ==3)
1090       return(im_copy(perf_ERR));
1091   else if (perf_type ==4)
1092       return(im_copy(perf_RCC));
1093 
1094   return(NULL);
1095 }
1096 
1097 Imrect *dfit_get_image(int diff_type)                    /* HR: new */
1098 {
1099   if (diff_type ==0)
1100       return(im_copy(diff_S0));
1101   else if (diff_type ==1)
1102       return(im_copy(diff_ADC));
1103 
1104   return(NULL);
1105 }
1106 
1107 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.