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

Linux Cross Reference
Tina4/src/tools/nmr-segment/p_gamma_fit.c

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

  1 /*
  2   p_gamma_fit.c
  3 
  4   Xiaoping
  5   a.lacey   13.5.98
  6   NAT 
  7 */
  8 
  9 #include <stdio.h>
 10 #include <math.h>
 11 #include <tina/sys.h>
 12 #include <tina/sysfuncs.h>
 13 #include <tina/math.h>
 14 #include <tina/mathfuncs.h>
 15 #include <tina/vision.h>
 16 #include <tina/visionfuncs.h>
 17 #include <tina/tv.h>
 18 #include <tina/tvfuncs.h>
 19 #include <tina/draw.h>
 20 #include <tina/drawfuncs.h>
 21 #include <tina/seqoral.h>
 22 #include <tina/seqpro.h>
 23 #include <tina/stimdef.h>
 24 
 25 #define NP 3
 26 #define MIN_BOLUS_WIDTH NP+1
 27 #define SEC_PEAK  2.0
 28 #define T0SCAN  10 
 29 #define TTPSCAN  15
 30 #define RR_INIT  3.6
 31 #define MIN_TTP_FRAC  0.5
 32 #define TTP_FRAC_STEP 0.1
 33 
 34 #define FTOL 1.0e-6
 35 
 36 static float SCL = 10000.0;
 37 static float MIN_T0 = 18.0;
 38 static float AVE_TTP = 7.0;
 39 static float RECIRC_CUT = 0.2;
 40 static float RECIRC_PERIOD = 10.0;
 41 
 42 static float MAX_TTP;
 43 static Perfusion *perf_global = NULL;
 44 static float *our_time;
 45 
 46 extern void ***coreg_slice_init(Vec3 *ptr);
 47 extern void ***coreg_limits(int *lz, int *uz, int *ly, int *uy, int *lx, int *ux);
 48 extern float this_pixel(void ***rasptrs, int z, int y, int x);
 49 
 50 static   int     imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
 51 
 52 
 53 float  *s_2d_mask_mean(void ***imptrs, Imrect *mask, Perfusion *p)
 54 {
 55   float   *mean, area_sum;
 56   int      i, j, k, total;
 57   float  *row;
 58   double val;
 59   
 60 
 61   if ((mean = fvector_alloc(p->n1, p->n2)) == NULL)
 62     return(NULL);
 63 
 64   for (k = p->n1; k < p->n2; k++)
 65   {
 66     area_sum = 0.0;
 67     total = 0;
 68 
 69     for (i = imptrly; i < imptruy; i++)
 70     {
 71       for (j = imptrlx; j < imptrux; j++)
 72       {
 73          IM_PIX_GET(mask, i, j, val); 
 74         if (val == 1)
 75         {
 76           area_sum += this_pixel(imptrs, k, i, j);
 77           total++;
 78         }
 79       }
 80     }
 81     mean[k] = area_sum/(float)total;
 82   }
 83 
 84  return(mean);
 85 }
 86 
 87 
 88 
 89 float  *s_2d_pixel_get(void ***imptrs, Vec2 v, Perfusion *p)
 90 {
 91   float   *data;
 92   int      i, j, k;
 93   
 94 
 95   if ((data = fvector_alloc(p->n1, p->n2)) == NULL)
 96     return(NULL);
 97 
 98   j = tina_int(vec2_x(v));
 99   i = tina_int(vec2_y(v));
100 
101   for (k = p->n1; k < p->n2-1; k++)
102   { 
103        data[k] = 0.5*(this_pixel(imptrs, k, i, j)
104                   +this_pixel(imptrs, k+1, i, j));
105   }
106   data[k] = this_pixel(imptrs, k, i, j);
107 
108   return(data);
109 }
110   
111 int find_bolus_time(Perfusion *p)
112 {
113    int     peak_time, s0_time;
114    int k;
115    int bolus_time;
116    float max_val, min_val, perc35_val;
117    float *mean = p->st_1d;
118 
119    min_val = fvector_min(mean, p->n1+MIN_T0, p->n2-AVE_TTP);
120 
121    for (k = p->n1; k < p->n2-AVE_TTP; k++)
122    {
123      if (mean[k] <= min_val)
124         peak_time = k;
125    }
126 
127    max_val = fvector_max(mean, p->n1 + MIN_T0, peak_time);
128 
129    s0_time = peak_time - AVE_TTP;
130 
131    perc35_val = max_val-0.35*(max_val-min_val);
132 
133    for(k=peak_time; k>p->n1+MIN_T0; k--){
134       if((mean[k]>=perc35_val)&&(mean[k-1]>=perc35_val)){
135          s0_time=k;
136          break;
137       }
138    }
139 
140    if (s0_time < MIN_T0 || s0_time > p->n2-AVE_TTP)
141    {
142        s0_time = (p->n2 -p->n1)/2-AVE_TTP;
143    }
144 
145    return(s0_time);
146 }
147   
148 
149 float ave_s0_pre_bolus(Perfusion *p)
150 {
151   float *st_1d = p->st_1d;
152   int    k , bolus_time = p->bolus_time;
153   float  noisesum = 0.0;
154 
155   for (k = p->n1; k < bolus_time-2 && k < p->n2; k++) 
156     noisesum += st_1d[k];
157 
158 
159   return (noisesum/(bolus_time-2-p->n1));
160 }
161 
162 float conv_r2_to_st(float r2, float s0, float te)
163 {
164     float st;
165     st = s0*exp(-r2*te/SCL);
166     return (st);
167 }
168 
169 float *conv_st_to_r2(Perfusion *p)
170 {
171   float  *st_1d = p->st_1d;
172   int     n1 = p->n1;
173   int     n2 = p->n2;
174   float   s0 = p->s0_av;
175   float   te = p->te;
176   float  *r2, *new_r2, val;
177   int     imax, i,j,k;
178 
179   if ((r2 = fvector_alloc(n1, n2)) == NULL)
180     {
181       error("conv_st_to_r2: memory alloc error ", warning);
182       return;
183     }
184 
185     for (k = n1; k < n2; k++)
186     {
187       val = SCL*(-log(st_1d[k]/s0))/te;
188       if (isnan(val))
189         {
190           r2[k] = 0.0;
191                                         printf("%f\n", val);
192         } 
193       else
194         r2[k] = val;
195       val = conv_r2_to_st(val, s0, te);
196     }
197 
198   return(r2);
199 }
200 
201 double gamma_func( double ka0, double ka1, double ka2,
202                    float t, float termi)
203 {
204    double t1 = t - 0.5 *termi ,t2 = t+0.5*termi;
205    double con=0.0, weight;
206 
207    if (t1<=0.0)
208    {
209       if (t2<=0.0)
210       {
211          con = 0.0;
212       }
213       else
214       {
215          weight = t2/termi;
216          con = 0.5*weight*ka0*pow(t2,ka1)*exp(-t2/ka2);
217       }
218    }
219    else
220    {
221       con += 0.5*ka0*pow(t1,ka1)*exp(-t1/ka2);
222       con += 0.5*ka0*pow(t2,ka1)*exp(-t2/ka2);
223    }
224    return(con);
225 }
226 
227 double cut_off(double ka1, double ka2, float cut)
228 {
229    double ans1 , ans2 = 50;
230 
231    while (fabs(ans1 - ans2) > 0.001)
232    {
233        ans1 = ans2;
234        ans2 = ka2*(ka1*log(ans1)-log(cut));
235    }
236    if (ans1>0.0)
237       return(ans1);
238    else
239       return(0.0);
240 }
241 
242 void r2_ther_stats(Perfusion *perf_data, int offset,
243                    float *norm1, float *norm2, float *dot, double *chi)
244 {
245    double ka0,ka1,ka2, t0;
246    int n1 = perf_data->n1;
247    int n2 = perf_data->n2;
248    float termi = perf_data->tscale;
249    float *data = perf_data->r2;
250    float *r2_ther = perf_data->r2_ther;
251    float t;
252    double cut;
253    int i, icut, tcut;
254 
255    ka1 = perf_data->rr;  /* b  */
256    ka2 = perf_data->bb;  /* r  */
257    ka0 = perf_data->qq;  /* q  */
258    t0 = perf_data->bolus_time*termi-perf_data->tt;
259    tcut = perf_data->bolus_time;
260    icut = perf_data->bolus_end;
261 
262    if (offset == 0)
263    {
264       if (perf_data->r2_ther == NULL)
265            perf_data->r2_ther = fvector_alloc(n1, n2);
266 /*
267       else
268       {
269            for (i=n1;i<n2;i++)
270                perf_data->r2_ther[i] = 0.0;
271       }
272 */
273       r2_ther = perf_data->r2_ther;
274       *norm1 = 0.0;
275       for (i=tcut;
276            i<icut&&i<n2;
277            i++)
278       {
279          t = termi*i-t0;
280          r2_ther[i] = gamma_func(ka0,ka1,ka2,t,termi);
281          *norm1 += r2_ther[i]*r2_ther[i];
282       }
283       *norm1 = sqrt(*norm1);
284    }
285    *dot = 0.0;
286    *norm2 = 0.0;
287    for (i=tcut;
288         i<icut&&i+offset<n2;
289         i++)
290    {
291       *dot += r2_ther[i]*data[i+offset];
292       *norm2 += data[i+offset]*data[i+offset];
293    }
294    *chi = (*norm2)*(*dot/(*norm1*sqrt(*norm2)) -1.0)/(float)(i-NP-tcut);
295 }
296 
297 static double gamma_chi2(int ndim, double *x)
298 /* minimise errors on measured data */
299 {
300   float *data = perf_global->st_1d;
301   float termi = perf_global->tscale;
302   double ka0,ka1,ka2, tt;
303   double chi;
304   float dot,norm1,norm2;
305 
306   ka1=x[0]; ka2=x[1]; tt=x[2];
307 
308   if (ka1 > 0.0 && ka2 > 1.0e-6 && ka1*ka2 < MAX_TTP
309      && tt < 2.0*termi && tt > -2.0*termi );
310     /* protect against NaN and Infty etc */
311   else
312   {
313     perf_global->qq = 0.0;
314     return(MAXFLOAT);
315   }
316 
317   perf_global->qq = 1.0;
318   perf_global->bb = ka2;
319   perf_global->rr = ka1;
320   perf_global->tt = x[2];
321   r2_ther_stats(perf_global,0,&norm1,&norm2,&dot,&chi); 
322   if (dot>1.0 && norm1 > 0.0)
323   {
324     perf_global->qq *= dot/(norm1*norm1);
325     return(-chi);
326   }
327   else
328   {
329     perf_global->qq = 0.0;
330     return(MAXFLOAT);
331   }
332 }
333 
334 static gamma_simplex(Perfusion *p)
335 {
336   double   *x;
337   float termi = p->tscale;
338   double   *lambda;
339   int ndim = NP;
340   float cut,maxval;
341   
342   perf_global = p; /* access to parameters for gamma_chi() */
343   x = (double *)ralloc(ndim*sizeof(double));
344   x[0] = p->rr;
345   x[1] = p->bb;
346   x[2] = p->tt;
347 
348   lambda = (double *)ralloc(ndim*sizeof(double));
349   lambda[0] = 0.1;
350   lambda[1] = 0.1;
351   lambda[2] = 2.0;
352   simplexmin2(ndim, x, lambda, gamma_chi2, FTOL, /*(void (*) ()) format */ NULL);
353   p->rr=x[0];
354   p->bb=x[1];
355   p->tt=x[2];
356   gamma_chi2(ndim,x);
357   rfree(x);
358   rfree(lambda);
359 }
360 
361 float  emap(Perfusion *p)
362 {
363     int i, icut, tcut;
364     float cut, theory, ttot=0.0, dof, t;
365     double chi, con, chitot=0.0, frac, weight;
366     double ka0,  ka1,  ka2, t0;
367     float termi = p->tscale;
368     float *data = p->r2;
369 
370     ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
371     cut = cut_off(ka1,ka2,RECIRC_CUT*gamma_func(1.0,ka1,ka2,ka1*ka2,termi)); 
372     icut = tina_int(cut/termi); 
373     tcut = tina_int(t0/termi);
374     for (i=tcut,dof=-4.0;
375        i<(tcut+icut)&&i<p->n2;
376        dof+=1.0,i++)
377    {
378       t = termi*i-t0;
379       con= gamma_func(ka0,ka1,ka2,t,termi);
380       chi= data[i] - con;
381       chitot += chi*chi;
382    }
383    if (dof>0.0)
384       return(chitot/dof);
385    else
386       return(0.0);
387 }
388 
389 float  rcbv(Perfusion *p)
390 {
391    int i, tcut;
392    float cut, ttot=0.0, t;
393    double ka0, ka1, ka2, t0;
394    float termi = p->tscale;
395    double con, frac;
396     
397    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
398    tcut = tina_int(t0/termi);
399    for (i=tcut; i<p->n2; i++)
400    {
401       t = termi*i-t0;
402       con= gamma_func(ka0,ka1,ka2,t,termi);
403       ttot += con;
404    }
405    return(ttot);
406 }
407 
408 float  recirc(Perfusion *p)
409 {
410    int i, icut, tcut, rcut;
411    float cut, ttot=0.0, t;
412    double ka0, ka1, ka2, t0;
413    float termi = p->tscale;
414    float *data = p->r2;
415    double con, frac;
416    
417     ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
418     cut = cut_off(ka1,ka2,RECIRC_CUT*gamma_func(1.0,ka1,ka2,ka1*ka2,termi)); 
419     icut = tina_int(cut/termi); 
420     tcut = tina_int(t0/termi);
421     rcut = tina_int(RECIRC_PERIOD/termi);
422     for (i=tcut+icut; i<tcut+icut+rcut && i<p->n2; i++)
423     {
424        t = termi*i-t0;
425        con= gamma_func(ka0,ka1,ka2,t,termi);
426        ttot += data[i] - con;
427     }
428     return(ttot);
429 }
430 
431 float  bolus_mtt(Perfusion *p)
432 {
433    int i, tcut;
434    float cut, t;
435    double ka0, ka1, ka2, t0;
436    float termi = p->tscale;
437    double con, frac;
438    double mean=0.0;
439     
440    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
441    tcut = tina_int(t0/termi);
442    for (i=tcut; i<p->n2; i++)
443    {
444       t = termi*i-t0;
445       con= gamma_func(ka0,ka1,ka2,t,termi);
446       mean += con*t;
447    }
448    
449    return(mean);
450 }
451 
452 float  bolus_m0(Perfusion *p)
453 {
454    int i, tcut;
455    float cut, t;
456    double ka0, ka1, ka2, t0;
457    float termi = p->tscale;
458    double con, frac;
459    double mean=0.0;
460     
461    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
462    tcut = tina_int(t0/termi);
463    for (i=tcut; i<p->n2; i++)
464    {
465       t = termi*i-t0;
466       con= gamma_func(ka0,ka1,ka2,t,termi);
467       mean += con*(i+0.5)*termi;
468    }
469    
470    return(mean);
471 }
472 
473 float  bolus_m1(Perfusion *p)
474 {
475    int i, tcut;
476    float cut, t;
477    double ka0, ka1, ka2, t0;
478    float termi = p->tscale;
479    double con, frac;
480    double mean=0.0;
481 
482    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
483    tcut = tina_int(t0/termi);
484    for (i=tcut; i<p->n2; i++)
485    {
486       t = termi*i-t0;
487       con= gamma_func(ka0,ka1,ka2,t,termi);
488       mean += con*termi*(i+0.5)*termi*(i+0.5);
489    }
490    return(mean);
491 }
492 
493 float   *recon_thre_r2(float termi, Perfusion *p)
494 {
495   int     i, n;
496   float  *r2_ther, td , t0;
497 
498   if ((r2_ther = fvector_alloc(p->n1, p->n2)) == NULL)
499     {
500       error("recon_thre_r2: memory alloc error", warning);
501       return(NULL);
502     } 
503 
504   t0 = p->bolus_time*termi - p->tt;
505   for (i = p->n1; i < p->n2; i++) 
506   {
507     td = termi*i -t0;
508     r2_ther[i] = gamma_func(p->qq,p->rr,p->bb,td,termi);
509   }
510 
511   return(r2_ther);
512 }
513 
514 search_chi_fits(Perfusion *perf_data, float **chi2array,float **dotarray,
515                 float *norm1)
516 {
517    void get_quadf(float **pixval, int j1, int i1, float *y, float *x);
518    float kmax, jmax, dot, dotmax, norm1_max, norm2, chi;
519    int j,k,test;
520    float chimax;
521 
522    dotmax = -MAXFLOAT;
523    chimax = -MAXFLOAT;
524    kmax = T0SCAN-1;
525    for (j=0;j<TTPSCAN;j++)
526    {
527       for (k=0;k<T0SCAN;k++)
528       {
529          test = 0;
530          if (k>0                && chi2array[j][k]<chi2array[j][k-1]  )  continue; 
531          if (k>0 && j >0        && chi2array[j][k]<chi2array[j-1][k-1])  continue; 
532          if (k>0 && j<TTPSCAN-1 && chi2array[j][k]<chi2array[j+1][k-1])  continue;
533          if (k<T0SCAN-1        && chi2array[j][k]<chi2array[j][k+1]  )  continue;
534          if (k<T0SCAN-1 && j>0 && chi2array[j][k]<chi2array[j-1][k+1])  continue;
535          if (k<T0SCAN-1 && j<TTPSCAN-1 && chi2array[j][k]<chi2array[j+1][k+1])  continue;
536          if (j>0                && chi2array[j][k]<chi2array[j-1][k]  )  continue;
537          if (j<TTPSCAN-1        && chi2array[j][k]<chi2array[j+1][k]  )  continue;
538 
539          if ( ( dotarray[j][k] > dotmax && k < kmax) ||
540               ( dotarray[j][k] > SEC_PEAK*dotmax && k > kmax))
541          {
542             chimax = chi2array[j][k];
543             dotmax = dotarray[j][k];
544             norm1_max = norm1[j];
545             kmax = k;
546             jmax = j;
547          }
548       }
549    }
550    if (dotmax > 0)
551    {
552       float t0,cut;
553       int icut,tcut;
554       perf_data->qq *= dotmax/norm1_max;
555       perf_data->tt = -(kmax)*perf_data->tscale;
556       perf_data->bb = (MIN_TTP_FRAC+(TTP_FRAC_STEP*jmax))*AVE_TTP/perf_data->rr;
557       t0 = perf_data->tscale*perf_data->bolus_time-perf_data->tt;
558 
559       cut = cut_off(perf_data->rr,perf_data->bb,
560             RECIRC_CUT*gamma_func(1.0,perf_data->rr,perf_data->bb,
561             perf_data->rr*perf_data->bb,perf_data->tscale));
562 
563       icut = tina_int(cut/perf_data->tscale);
564       tcut = tina_int(t0/perf_data->tscale);
565       perf_data->bolus_time -= tina_int(perf_data->tt/perf_data->tscale);
566       perf_data->tt -= tina_int(perf_data->tt/perf_data->tscale)*perf_data->tscale;
567       perf_data->refit = 1;
568       perf_data->bolus_end = tcut+icut;
569       if (perf_data->bolus_end - perf_data->bolus_time < MIN_BOLUS_WIDTH)
570       {
571          perf_data->bolus_end = perf_data->bolus_time+MIN_BOLUS_WIDTH;
572 /*
573          perf_data->qq *= 0.0;
574          perf_data->tt = -T0SCAN*perf_data->tscale/2.0;
575          perf_data->bb = AVE_TTP/perf_data->rr;
576 */
577       }
578    }
579    else
580    {
581       perf_data->qq *= 0.0;
582       perf_data->tt = -T0SCAN*perf_data->tscale/2.0;
583       perf_data->bb = AVE_TTP/perf_data->rr;
584    }
585 }
586 
587 void gamma_match(Perfusion *perf_data)
588 {
589    int i,j,k;
590    float kmax,jmax;
591    double ka1,ka2;
592    float termi = perf_data->tscale;
593    float norm2,norm1_max;
594    float *norm1;
595    float dot,dotmax;
596    double chi;
597    float **chi2array;
598    float **dotarray;
599    double cut;
600    int   tcut;
601 
602    ka1 = perf_data->rr = RR_INIT;  /* b  */
603    perf_data->bb = T0SCAN/RR_INIT;  /* r  */
604    perf_data->qq = 1.0;  /* q  */
605    perf_data->tt = 0.0;
606    norm1 = fvector_alloc(0,TTPSCAN);
607    dotarray = farray_alloc(0,0,TTPSCAN,T0SCAN);
608    chi2array = farray_alloc(0,0,TTPSCAN,T0SCAN);
609 
610    for (ka2 = MIN_TTP_FRAC*AVE_TTP/ka1, j=0; j<TTPSCAN ;
611                        ka2+=TTP_FRAC_STEP*AVE_TTP/ka1,j++)
612    {
613       cut = cut_off(ka1,ka2,RECIRC_CUT*gamma_func(1.0,ka1,ka2,ka1*ka2,termi));
614       perf_data->bolus_end = perf_data->bolus_time + tina_int(cut/termi);
615 
616       perf_data->bb = ka2;
617       for (k=0;k<T0SCAN;k++)
618       {
619          r2_ther_stats(perf_data,k,&norm1[j],&norm2,&dot,&chi);
620          dotarray[j][k] = dot/(norm1[j]);
621          chi2array[j][k] = chi;
622       }
623    }
624    search_chi_fits(perf_data, chi2array,dotarray,norm1);
625    farray_free(chi2array,0,0,TTPSCAN,T0SCAN);
626    farray_free(dotarray,0,0,TTPSCAN,T0SCAN);
627    fvector_free(norm1,0);
628 }
629 
630 int   gamma_fit_pixel(Perfusion *perf_data, Vec2 v, void ***imptrs)
631 {
632   void get_perf_params(float *ave_ttp, float *min_t0, float *recirc_cut, 
633                        float *recirc_period, float *scl);
634   List   *store = (List *)get_start_el();
635   Imrect   *im = NULL;
636   Tv       *tv = (Tv *)imcalc_graph_tv_get();
637   float    *te;
638   int       i, j, total, type;
639 
640   if (tv == NULL || tv->tv_screen == NULL)
641     return(-1);
642   if ((store == NULL) || (store->to == NULL))
643     return(-1);
644   im = (Imrect *)(store->to);
645 
646   if ((our_time = (float *)prop_get(im->props, DYNSTIME)) == NULL)
647     {
648       error("gamma_fit: timing info missing", warning);
649       return(-1);
650     }
651 
652   perf_data->tscale = our_time[2]-our_time[1];
653   get_perf_params(&AVE_TTP, &MIN_T0, &RECIRC_CUT, &RECIRC_PERIOD, &SCL);
654   MIN_T0 /= perf_data->tscale;
655   MAX_TTP = AVE_TTP*(TTPSCAN*TTP_FRAC_STEP + MIN_TTP_FRAC);
656   if ((te = (float *)prop_get(im->props, TE_DATA)) == NULL)
657     {
658       error("gamma_fit: echo time not found", warning);
659       return(-1);
660     }
661 
662   coreg_limits(&perf_data->n1, &perf_data->n2, &imptrly, &imptruy, &imptrlx, &imptrux);
663 
664   perf_data->te = *te;
665   if (perf_data->st_1d!=NULL) fvector_free(perf_data->st_1d, perf_data->n1);
666   perf_data->st_1d = s_2d_pixel_get(imptrs,v,perf_data);
667   perf_data->bolus_time = MIN_T0;
668   perf_data->s0_av = ave_s0_pre_bolus(perf_data);
669   if (perf_data->r2!=NULL) fvector_free(perf_data->r2, perf_data->n1);
670   perf_data->r2 = conv_st_to_r2(perf_data);
671   gamma_match(perf_data);
672   if (perf_data->qq > 0.0 ) gamma_simplex(perf_data);
673   if (perf_data->r2_ther!=NULL) fvector_free(perf_data->r2_ther, perf_data->n1);
674   perf_data->r2_ther = recon_thre_r2(perf_data->tscale, perf_data);
675 
676   return (0);
677 }
678 
679 double  gamma_fit_region(Perfusion *perf_data)
680 {
681   List     *tv_poly_get();
682   void get_perf_params(float *ave_ttp, float *min_t0, float *recirc_cut,
683                        float *recirc_period, float *scl);
684   List     *poly = NULL;
685   List   *store = (List *)get_start_el();
686   Imrect   *mask = NULL, *tmask = NULL, *fmask = NULL, *im;
687   Tv       *tv = (Tv *)imcalc_graph_tv_get();
688   float    *te;
689   int       i, j, total, type;
690   void   ***imptrs=NULL;
691   double t0;
692 
693 
694   if (tv == NULL || tv->tv_screen == NULL)
695     return;
696   if ((store == NULL) || (store->to == NULL))
697     return;
698   im = (Imrect *)(store->to);
699 
700   if ((our_time = (float *)prop_get(im->props, DYNSTIME)) == NULL)
701     {
702       error("plot_aveflow_proc: timing info missing", warning);
703       return;
704     }
705   perf_data->tscale = our_time[2]-our_time[1];
706   get_perf_params(&AVE_TTP, &MIN_T0, &RECIRC_CUT, &RECIRC_PERIOD, &SCL);
707   MIN_T0 /= perf_data->tscale;
708   MAX_TTP = AVE_TTP*(TTPSCAN*TTP_FRAC_STEP + MIN_TTP_FRAC);
709 
710   if ((te = (float *)prop_get(im->props, TE_DATA)) == NULL)
711     {
712       error("gamma_fit: echo time not found", warning);
713       return;
714     }
715 
716   if (stack_check_types(IMRECT, NULL) == false)
717     {
718       error("gamma_fit: wrong type on stack", warning);
719       return;
720     }
721   mask = (Imrect *) stack_pop(&type);
722 
723   if ((tmask = im_alloc(mask->height, mask->width, NULL, char_v)) == NULL)
724     {
725       error("gamma_fit_region: memory alloc error 2", warning);
726       stack_push((void *)mask, IMRECT, im_free);
727       return;
728     }
729 
730   for (i = mask->region->ly; i < mask->region->uy; i++)
731     for (j = mask->region->lx; j < mask->region->ux; j++)
732       IM_CHAR(tmask, i, j) = 1;
733 /*
734   if ((poly = tv_poly_get()) != NULL)
735     im_poly_crop(tmask, poly);
736 */
737   fmask = im_prod(tmask, mask);
738 
739   imptrs = coreg_slice_init((Vec3 *)NULL);
740   coreg_limits(&perf_data->n1, &perf_data->n2, &imptrly, &imptruy, &imptrlx, &imptrux);
741 
742   perf_data->te = *te;
743   if (perf_data->st_1d!=NULL) fvector_free(perf_data->st_1d, perf_data->n1);
744   perf_data->st_1d = s_2d_mask_mean(imptrs, fmask, perf_data);
745   perf_data->bolus_time = MIN_T0;
746   perf_data->s0_av = ave_s0_pre_bolus(perf_data);
747   if (perf_data->r2!=NULL) fvector_free(perf_data->r2, perf_data->n1);
748   perf_data->r2 = conv_st_to_r2(perf_data); 
749   perf_data->tscale = our_time[2]-our_time[1];
750   gamma_match(perf_data);
751   if (perf_data->qq > 0.0 ) gamma_simplex(perf_data);
752   if (perf_data->r2_ther !=NULL) fvector_free(perf_data->r2_ther, perf_data->n1);
753   perf_data->r2_ther = recon_thre_r2(perf_data->tscale, perf_data); 
754 
755   stack_push((void *)mask, IMRECT, im_free);
756 
757   im_free(fmask);
758   im_free(tmask);
759   
760   t0 = perf_data->bolus_time*perf_data->tscale-perf_data->tt;
761   return (t0);
762 }
763 

~ [ 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.