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

Linux Cross Reference
Tina5/tina-libs/tina/medical/medStim_ct_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    :  
 23  * Date    :  
 24  * Version :  
 25  * CVS Id  :  
 26  *
 27  * Notes: mjs copied from medStim_gamma.c 4/1/06
 28  *
 29  *
 30  *  
 31  *
 32  *********
 33 */
 34 
 35 #if HAVE_CONFIG_H
 36 #   include <config.h>
 37 #endif
 38 
 39 #include "medStim_ct_gamma.h"
 40 
 41 #include <stdio.h>
 42 #include <math.h>
 43 #include <float.h>
 44 #include <tina/sys/sysDef.h>
 45 #include <tina/sys/sysPro.h>
 46 #include <tina/math/mathDef.h>
 47 #include <tina/math/mathPro.h>
 48 #include <tina/image/imgDef.h>
 49 #include <tina/image/imgPro.h>
 50 #include <tina/medical/med_StimDef.h>
 51 #include <tina/medical/medStim_alloc.h>
 52 #include <tina/medical/medStim_tdata.h>
 53 
 54 #define NP 3
 55 #define MIN_BOLUS_WIDTH NP+1
 56 #define SEC_PEAK  2.0
 57 #define T0SCAN  10 
 58 #define TTPSCAN  15
 59 #define RR_INIT  3.6
 60 #define MIN_TTP_FRAC  0.5
 61 #define TTP_FRAC_STEP 0.1
 62 
 63 #define FTOL 1.0e-6
 64 
 65 static double SCL = 10000.0;
 66 static double MIN_T0 = 18.0;
 67 static double AVE_MTT = 7.0;
 68 static double RECIRC_CUT = 0.2;
 69 static double RECIRC_PERIOD = 10.0;
 70 static Imrect *perf_TTP=NULL;
 71 static Imrect *perf_CBV=NULL;
 72 static Imrect *perf_MTT=NULL;
 73 static Imrect *perf_ERR=NULL;
 74 static Imrect *perf_RCC=NULL;
 75 
 76 static float MAX_TTP;
 77 static Perfusion *perf_global = NULL;
 78 
 79 
 80 double *ct_scl_get()
 81 {
 82    return(&SCL);
 83 }
 84 
 85 double *ct_min_t0_get()
 86 {
 87    return(&MIN_T0);
 88 }
 89 
 90 double *ct_ave_mtt_get()
 91 {
 92    return(&AVE_MTT);
 93 }
 94 
 95 double *ct_recirc_cut_get()
 96 {
 97    return(&RECIRC_CUT);
 98 }
 99 
100 double *ct_recirc_period_get()
101 {
102    return(&RECIRC_PERIOD);
103 }
104 
105 static float ct_ave_s0_pre_bolus(Perfusion *p)
106 {
107   float *st_1d = p->st_1d;
108   int    k , bolus_time = p->bolus_time;
109   float  noisesum = 0.0;
110 
111   if (bolus_time-2 <= p->n1) bolus_time = p->n1+3;
112 
113   for (k = p->n1; k < bolus_time-2 && k < p->n2; k++)
114     noisesum += st_1d[k];
115 
116 
117   return (noisesum/(bolus_time-2-p->n1));
118 }
119 
120 /*static float ct_conv_r2_to_st(float r2, float s0, float te)
121 {
122     float st;
123     st = s0*exp(-r2*te/SCL);
124     return (st);
125 }
126 */
127 
128 static float *ct_conv_st_to_r2(Perfusion *p)
129 {
130   float  *st_1d = p->st_1d;
131   int     n1 = p->n1;
132   int     n2 = p->n2;
133   float   s0 = p->s0_av;
134   float  *r2, val;
135   int      k;
136 
137   if ((r2 = fvector_alloc(n1, n2)) == NULL)
138     {
139       error("conv_st_to_r2: memory alloc error ", warning);
140       return(NULL);
141     }
142 
143     for (k = n1; k < n2; k++)
144     {
145       /*  float   te = p->te;
146          val = SCL*(-log(st_1d[k]/s0))/te;
147     */
148       val = st_1d[k]-s0; 
149       if (isnan(val))
150         {
151           r2[k] = 0.0;
152           printf("%f\n", val);
153         } 
154       else
155         r2[k] = val;
156       /*val = conv_r2_to_st(val, s0, te);*/
157     }
158 
159   return(r2);
160 }
161 
162 
163 static double ct_gamma_func( double ka0, double ka1, double ka2,
164                    float t, float termi)
165 {
166    double t1 = t - 0.5 *termi ,t2 = t+0.5*termi;
167    double con=0.0, weight;
168 
169    if (t1<=0.0)
170    {
171       if (t2<=0.0)
172       {
173          con = 0.0;
174       }
175       else
176       {
177          weight = t2/termi;
178          con = 0.5*weight*ka0*pow(t2,ka1)*exp(-t2/ka2);
179       }
180    }
181    else
182    {
183       con += 0.5*ka0*pow(t1,ka1)*exp(-t1/ka2);
184       con += 0.5*ka0*pow(t2,ka1)*exp(-t2/ka2);
185    }
186    return(con);
187 }
188 
189 /* grr. don't really know about this */
190 static double ct_cut_off(double ka1, double ka2, float cut)
191 {
192    double ans1=0, ans2=50;
193 
194    while (fabs(ans1 - ans2) > 0.001)
195    {
196        ans1 = ans2;
197        ans2 = ka2*(ka1*log(ans1)-log(cut));
198    }
199    if (ans1>0.0)
200       return(ans1);
201    else
202       return(0.0);
203 }
204 
205 static void ct_r2_ther_stats(Perfusion *perf_data, int offset,
206                           float *norm1, float *norm2, float *dot, double *chi)
207 {
208    double ka0,ka1,ka2, t0;
209    int n1 = perf_data->n1;
210    int n2 = perf_data->n2;
211    float termi = perf_data->tscale;
212    float *data = perf_data->r2;
213    /*float *data = perf_data->st_1d;*/
214    float *r2_ther = perf_data->r2_ther;
215    float t;
216    int i, icut, tcut;
217 
218    ka1 = perf_data->rr;  /* b  */
219    ka2 = perf_data->bb;  /* r  */
220    ka0 = perf_data->qq;  /* q  */
221    t0 = perf_data->bolus_time*termi-perf_data->tt;
222    tcut = perf_data->bolus_time;
223    icut = perf_data->bolus_end;
224 
225    if (offset == 0)
226    {
227       if (perf_data->r2_ther == NULL)
228            perf_data->r2_ther = fvector_alloc(n1, n2);
229       r2_ther = perf_data->r2_ther;
230       *norm1 = 0.0;
231       for (i=tcut;
232            i<icut&&i<n2;
233            i++)
234       {
235          t = termi*i-t0;
236          r2_ther[i] = ct_gamma_func(ka0,ka1,ka2,t,termi);
237          *norm1 += r2_ther[i]*r2_ther[i];
238       }
239       *norm1 = sqrt(*norm1);
240    }
241    *dot = 0.0;
242    *norm2 = 0.0;
243    for (i=tcut;
244         i<icut&&i+offset<n2;
245         i++)
246    {
247       *dot += r2_ther[i]*data[i+offset];
248       *norm2 += data[i+offset]*data[i+offset];
249    }
250    *chi = (*norm2)*(*dot/(*norm1*sqrt(*norm2)) -1.0)/(float)(i-NP-tcut);
251 }
252 
253 static double ct_gamma_chi2(int ndim, double *x)
254 /* minimise errors on measured data */
255 {
256   float termi = perf_global->tscale;
257   double ka1,ka2, tt;
258   double chi;
259   float dot,norm1,norm2;
260 
261   ka1=x[0]; ka2=x[1]; tt=x[2];
262 
263   if (ka1 > 0.0 && ka2 > 1.0e-6 && ka1*ka2 < MAX_TTP
264      && tt < 2.0*termi && tt > -2.0*termi );
265     /* protect against NaN and Infty etc */
266   else
267   {
268     perf_global->qq = 0.0;
269     return(FLT_MAX);
270   }
271 
272   perf_global->qq = 1.0;
273   perf_global->bb = ka2;
274   perf_global->rr = ka1;
275   perf_global->tt = x[2];
276   ct_r2_ther_stats(perf_global,0,&norm1,&norm2,&dot,&chi); 
277   /* it's life jim but not as we know it */
278   if (dot>1.0 && norm1 > 0.0)
279   {
280     perf_global->qq *= dot/(norm1*norm1);
281     return(-chi);
282   }
283   else
284   {
285     perf_global->qq = 0.0;
286     return(FLT_MAX);
287   }
288 }
289 
290 static void ct_gamma_simplex(Perfusion *p)
291 {
292   double   *x;
293   double   *lambda;
294   int ndim = NP;
295   
296   perf_global = p; /* access to parameters for gamma_chi() */
297   x = (double *)ralloc(ndim*sizeof(double));
298   x[0] = p->rr;
299   x[1] = p->bb;
300   x[2] = p->tt;
301 
302   lambda = (double *)ralloc(ndim*sizeof(double));
303   lambda[0] = 0.1;
304   lambda[1] = 0.1;
305   lambda[2] = 2.0;
306   simplexmin2(ndim, x, lambda, ct_gamma_chi2, NULL, FTOL, /*(void (*) ()) format */ NULL);
307   p->rr=x[0];
308   p->bb=x[1];
309   p->tt=x[2];
310   ct_gamma_chi2(ndim,x);
311   rfree(x);
312   rfree(lambda);
313 }
314 
315 
316 static float  ct_emap(Perfusion *p)
317 {
318     int i, icut, tcut;
319     float cut, dof, t;
320     double chi, con, chitot=0.0;
321     double ka0,  ka1,  ka2, t0;
322     float termi = p->tscale;
323     float *data = p->r2;
324 
325     ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
326     cut = ct_cut_off(ka1,ka2,RECIRC_CUT*ct_gamma_func(1.0,ka1,ka2,ka1*ka2,termi)); 
327     icut = tina_int(cut/termi); 
328     tcut = tina_int(t0/termi);
329     for (i=tcut,dof=-4.0;
330        i<(tcut+icut)&&i<p->n2;
331        dof+=1.0,i++)
332    {
333       t = termi*i-t0;
334       con= ct_gamma_func(ka0,ka1,ka2,t,termi);
335       chi= data[i] - con;
336       chitot += chi*chi;
337    }
338    if (dof>0.0)
339       return(chitot/dof);
340    else
341       return(0.0);
342 }
343 
344 static float  ct_rcbv(Perfusion *p)
345 {
346    int i, tcut;
347    float ttot=0.0, t;
348    double ka0, ka1, ka2, t0;
349    float termi = p->tscale;
350    double con;
351     
352    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
353    tcut = tina_int(t0/termi);
354    for (i=tcut; i<p->n2; i++)
355    {
356       t = termi*i-t0;
357       con= ct_gamma_func(ka0,ka1,ka2,t,termi);
358       ttot += con;
359    }
360    return(ttot);
361 }
362 
363 static float  ct_recirc(Perfusion *p)
364 {
365    int i, icut, tcut, rcut;
366    float cut, ttot=0.0, t;
367    double ka0, ka1, ka2, t0;
368    float termi = p->tscale;
369    float *data = p->r2;
370    double con;
371    
372     ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
373     cut = ct_cut_off(ka1,ka2,RECIRC_CUT*ct_gamma_func(1.0,ka1,ka2,ka1*ka2,termi)); 
374     icut = tina_int(cut/termi); 
375     tcut = tina_int(t0/termi);
376     rcut = tina_int(RECIRC_PERIOD/termi);
377     for (i=tcut+icut; i<tcut+icut+rcut && i<p->n2; i++)
378     {
379        t = termi*i-t0;
380        con= ct_gamma_func(ka0,ka1,ka2,t,termi);
381        ttot += data[i] - con;
382     }
383     return(ttot);
384 }
385 
386 static float  ct_bolus_m0(Perfusion *p)
387 {
388    int i, tcut;
389    float t;
390    double ka0, ka1, ka2, t0;
391    float termi = p->tscale;
392    double con;
393    double mean=0.0;
394     
395    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
396    tcut = tina_int(t0/termi);
397    for (i=tcut; i<p->n2; i++)
398    {
399       t = termi*i-t0;
400       con= ct_gamma_func(ka0,ka1,ka2,t,termi);
401       mean += con*(i+0.5)*termi;
402    }
403    
404    return(mean);
405 }
406 
407 static float  ct_bolus_m1(Perfusion *p)
408 {
409    int i, tcut;
410    float t;
411    double ka0, ka1, ka2, t0;
412    float termi = p->tscale;
413    double con;
414    double mean=0.0;
415 
416    ka0=p->qq; ka1=p->rr; ka2=p->bb; t0=termi*p->bolus_time-p->tt;
417    tcut = tina_int(t0/termi);
418    for (i=tcut; i<p->n2; i++)
419    {
420       t = termi*i-t0;
421       con= ct_gamma_func(ka0,ka1,ka2,t,termi);
422       mean += con*termi*(i+0.5)*termi*(i+0.5);
423    }
424    return(mean);
425 }
426 
427 static float   *ct_recon_ther_r2(float termi, Perfusion *p)
428 {
429   int     i;
430   float  *r2_ther, td , t0;
431 
432   if ((r2_ther = fvector_alloc(p->n1, p->n2)) == NULL)
433     {
434       error("recon_ther_r2: memory alloc error", warning);
435       return(NULL);
436     } 
437 
438   t0 = p->bolus_time*termi - p->tt;
439   for (i = p->n1; i < p->n2; i++) 
440   {
441     td = termi*i -t0;
442     r2_ther[i] = ct_gamma_func(p->qq,p->rr,p->bb,td,termi);
443   }
444 
445   return(r2_ther);
446 }
447 
448 static void ct_search_chi_fits(Perfusion *perf_data, float **chi2array,float **dotarray,
449                             float *norm1)
450 {
451    float kmax, jmax=0, dotmax, norm1_max=0;
452    int j,k,test;
453    float chimax;
454 
455    dotmax = -FLT_MAX;
456    chimax = -FLT_MAX;
457    kmax = T0SCAN-1;
458    for (j=0;j<TTPSCAN;j++)
459    {
460       for (k=0;k<T0SCAN;k++)
461       {
462          test = 0;
463          if (k>0                && chi2array[j][k]<chi2array[j][k-1]  )  continue; 
464          if (k>0 && j >0        && chi2array[j][k]<chi2array[j-1][k-1])  continue; 
465          if (k>0 && j<TTPSCAN-1 && chi2array[j][k]<chi2array[j+1][k-1])  continue;
466          if (k<T0SCAN-1        && chi2array[j][k]<chi2array[j][k+1]  )  continue;
467          if (k<T0SCAN-1 && j>0 && chi2array[j][k]<chi2array[j-1][k+1])  continue;
468          if (k<T0SCAN-1 && j<TTPSCAN-1 && chi2array[j][k]<chi2array[j+1][k+1])  continue;
469          if (j>0                && chi2array[j][k]<chi2array[j-1][k]  )  continue;
470          if (j<TTPSCAN-1        && chi2array[j][k]<chi2array[j+1][k]  )  continue;
471 
472          if ( ( dotarray[j][k] > dotmax && k < kmax) ||
473               ( dotarray[j][k] > SEC_PEAK*dotmax && k > kmax))
474          {
475             chimax = chi2array[j][k];
476             dotmax = dotarray[j][k];
477             norm1_max = norm1[j];
478             kmax = k;
479             jmax = j;
480          }
481       }
482    }
483    if (dotmax > 0)
484    {
485       float t0,cut;
486       int icut,tcut;
487       perf_data->qq *= dotmax/norm1_max;
488       perf_data->tt = -(kmax)*perf_data->tscale;
489       perf_data->bb = (MIN_TTP_FRAC+(TTP_FRAC_STEP*jmax))*AVE_MTT/perf_data->rr;
490       t0 = perf_data->tscale*perf_data->bolus_time-perf_data->tt;
491 
492       cut = ct_cut_off(perf_data->rr,perf_data->bb,
493             RECIRC_CUT*ct_gamma_func(1.0,perf_data->rr,perf_data->bb,
494             perf_data->rr*perf_data->bb,perf_data->tscale));
495 
496       icut = tina_int(cut/perf_data->tscale);
497       tcut = tina_int(t0/perf_data->tscale);
498       perf_data->bolus_time -= tina_int(perf_data->tt/perf_data->tscale);
499       perf_data->tt -= tina_int(perf_data->tt/perf_data->tscale)*perf_data->tscale;
500       perf_data->refit = 1;
501       perf_data->bolus_end = tcut+icut;
502       if (perf_data->bolus_end - perf_data->bolus_time < MIN_BOLUS_WIDTH)
503       {
504          perf_data->bolus_end = perf_data->bolus_time+MIN_BOLUS_WIDTH;
505 /*
506          perf_data->qq *= 0.0;
507          perf_data->tt = -T0SCAN*perf_data->tscale/2.0;
508          perf_data->bb = AVE_MTT/perf_data->rr;
509 */
510       }
511    }
512    else
513    {
514       perf_data->qq *= 0.0;
515       perf_data->tt = -T0SCAN*perf_data->tscale/2.0;
516       perf_data->bb = AVE_MTT/perf_data->rr;
517    }
518 }
519 
520 static void ct_gamma_match(Perfusion *perf_data)
521 {
522    int j,k;
523    double ka1,ka2;
524    float termi = perf_data->tscale;
525    float norm2;
526    float *norm1;
527    float dot;
528    double chi;
529    float **chi2array;
530    float **dotarray;
531    double cut;
532 
533    ka1 = perf_data->rr = RR_INIT;  /* b  */
534    perf_data->bb = T0SCAN/RR_INIT;  /* r  */
535    perf_data->qq = 1.0;  /* q  */
536    perf_data->tt = 0.0;
537    norm1 = fvector_alloc(0,TTPSCAN);
538    dotarray = farray_alloc(0,0,TTPSCAN,T0SCAN);
539    chi2array = farray_alloc(0,0,TTPSCAN,T0SCAN);
540 
541    for (ka2 = MIN_TTP_FRAC*AVE_MTT/ka1, j=0; j<TTPSCAN ;
542                        ka2+=TTP_FRAC_STEP*AVE_MTT/ka1,j++)
543    {
544       cut = ct_cut_off(ka1,ka2,RECIRC_CUT*ct_gamma_func(1.0,ka1,ka2,ka1*ka2,termi));
545       perf_data->bolus_end = perf_data->bolus_time + tina_int(cut/termi);
546 
547       perf_data->bb = ka2;
548       for (k=0;k<T0SCAN;k++)
549       {
550         /*need to keep this, but it will need modifying as a function.. 
551          it's not really going to be r2 anymore, but there we are*/
552          ct_r2_ther_stats(perf_data,k,&norm1[j],&norm2,&dot,&chi);
553          dotarray[j][k] = dot/(norm1[j]);
554          chi2array[j][k] = chi;
555       }
556    }
557    ct_search_chi_fits(perf_data, chi2array,dotarray,norm1);
558    farray_free(chi2array,0,0,TTPSCAN,T0SCAN);
559    farray_free(dotarray,0,0,TTPSCAN,T0SCAN);
560    fvector_free(norm1,0);
561 }
562 
563 int   ct_gamma_fit_pixel(Perfusion *perf_data, Vec2 v, void ***imptrs)
564 {
565   Sequence *seq= (Sequence *)seq_get_current();
566   List   *store = (List *)get_current_seq_start_el();
567   Imrect   *im = NULL;
568   float our_time;       
569 
570   if ((store == NULL) || (store->to == NULL))
571     return(-1);
572   im = (Imrect *)(store->to);
573 
574   if ((our_time = seq->dim[3]) == 0.0)
575     {
576       error("gamma_fit: timing info missing", warning);
577       return(-1);
578     }
579 
580   perf_data->tscale = our_time;
581   MAX_TTP = AVE_MTT*(TTPSCAN*TTP_FRAC_STEP + MIN_TTP_FRAC);
582   /*
583     float    *te=NULL;
584 
585     if ((te = (float *)prop_get(seq->props, TE_DATA)) == NULL)
586     {
587       error("gamma_fit: echo time not found", warning);
588       return(-1);
589     }
590   */
591   if (perf_data->st_1d!=NULL) fvector_free(perf_data->st_1d, perf_data->n1);
592   if (perf_data->r2!=NULL) fvector_free(perf_data->r2, perf_data->n1);
593   if (perf_data->r2_ther!=NULL) fvector_free(perf_data->r2_ther, perf_data->n1);
594   perf_data->n1 = seq->offset;
595   perf_data->n2 = get_end_frame(seq);
596 
597   /*mjs 8/11/05 assumes that all TE's in sequence are same as that of first image in sequence */
598   /* perf_data->te = te[seq->offset];*/
599   perf_data->st_1d = s_2d_pixel_get(imptrs, v, perf_data);
600   perf_data->bolus_time = MIN_T0/perf_data->tscale;
601   perf_data->s0_av = ct_ave_s0_pre_bolus(perf_data);
602   perf_data->r2 = ct_conv_st_to_r2(perf_data);
603   ct_gamma_match(perf_data);
604   if (perf_data->qq > 0.0 ) ct_gamma_simplex(perf_data);
605   perf_data->r2_ther = ct_recon_ther_r2(perf_data->tscale, perf_data);
606 
607   return (0);
608 }
609 
610 /*
611  *   Altered to improve lib / tool separation - GAB 19 Feb 2003
612  *   Stack manipulation moved to pixel_gfit_region, mask passed as parameter.
613  */
614 double  ct_gamma_fit_region(Perfusion *perf_data, Imrect *mask)
615 {
616   Sequence *seq= (Sequence *)seq_get_current();
617   List     *store = (List *)get_current_seq_start_el();
618   Imrect   *tmask = NULL, *fmask = NULL, *im;
619   int       i, j;
620   void   ***imptrs=NULL;
621   double    t0;
622   float our_time;       
623 
624   if ((store == NULL) || (store->to == NULL))
625     return(0.0);
626   im = (Imrect *)(store->to);
627 
628   if ((our_time = seq->dim[3]) == 0.0)
629     {
630       error("gamma_fit: timing info missing", warning);
631       return(-1);
632     }
633 
634   perf_data->tscale = our_time;
635   MAX_TTP = AVE_MTT*(TTPSCAN*TTP_FRAC_STEP + MIN_TTP_FRAC);
636 
637   /*
638     float    *te=NULL;
639 
640     if ((te = (float *)prop_get(seq->props, TE_DATA)) == NULL)
641     {
642     error("gamma_fit: echo time not found", warning);
643     return(0.0);
644     }
645   */
646 
647   if ((tmask = im_alloc(mask->height, mask->width, NULL, char_v)) == NULL)
648     {
649       error("gamma_fit_region: memory alloc error 2", warning);
650       return(0.0);
651     }
652 
653   for (i = mask->region->ly; i < mask->region->uy; i++)
654     for (j = mask->region->lx; j < mask->region->ux; j++)
655       IM_CHAR(tmask, i, j) = 1;
656   fmask = im_prod(tmask, mask);
657 
658   if (perf_data->st_1d!=NULL) fvector_free(perf_data->st_1d, perf_data->n1);
659   if (perf_data->r2!=NULL) fvector_free(perf_data->r2, perf_data->n1);
660   if (perf_data->r2_ther !=NULL) fvector_free(perf_data->r2_ther, perf_data->n1);
661   perf_data->n1 = seq->offset;
662   perf_data->n2 = get_end_frame(seq);
663 
664   imptrs = seq_slice_init(seq);
665   /*perf_data->te = te[seq->offset];*/
666   perf_data->st_1d = s_2d_mask_mean(imptrs, fmask, perf_data);
667   perf_data->bolus_time = MIN_T0/perf_data->tscale;
668   perf_data->s0_av = ct_ave_s0_pre_bolus(perf_data);
669   perf_data->r2 = ct_conv_st_to_r2(perf_data); 
670   perf_data->tscale = our_time;
671   ct_gamma_match(perf_data);
672   if (perf_data->qq > 0.0 ) ct_gamma_simplex(perf_data);
673   perf_data->r2_ther = ct_recon_ther_r2(perf_data->tscale, perf_data); 
674 
675   im_free(fmask);
676   im_free(tmask);
677   
678   t0 = perf_data->bolus_time*perf_data->tscale-perf_data->tt;
679   format("T_zero: %f\n", t0);
680   return (t0);
681 }
682 
683 /*
684  *   Altered to improve lib / tool separation - GAB 19 Feb 2003
685  *   Stack manipulation moved to gfit_pixels_proc, mask passed as parameter.
686  */
687 void    ct_gfit_measure_image(double err_thresh, Imrect *mask)
688 {
689   Perfusion *p = NULL;
690   Sequence  *seq= (Sequence *)seq_get_current();
691   List      *store = (List *)get_current_seq_start_el();
692   Imrect    *im = NULL;
693   Imregion  *roi;
694   float      cbv, rcc, mtt, err, ttp;
695   Vec2       pos;
696   int        i, j;
697   void    ***imptrs=NULL;
698 
699   im = (Imrect *)(store->to);
700   roi = im->region;
701 
702   if (perf_TTP!=NULL) im_free(perf_TTP);
703   perf_TTP = im_alloc(im->height, im->width, roi, float_v);
704   if (perf_CBV!=NULL) im_free(perf_CBV);
705   perf_CBV = im_alloc(im->height, im->width, roi, float_v);
706   if (perf_MTT!=NULL) im_free(perf_MTT);
707   perf_MTT = im_alloc(im->height, im->width, roi, float_v);
708   if (perf_ERR!=NULL) im_free(perf_ERR);
709   perf_ERR = im_alloc(im->height, im->width, roi, float_v);
710   if (perf_RCC!=NULL) im_free(perf_RCC);
711   perf_RCC = im_alloc(im->height, im->width, roi, float_v);
712 
713 /*   Stack manipulation was here - GAB 19 Feb 2003   */
714   if (mask != NULL)
715     roi = mask->region;
716 
717   imptrs = seq_slice_init(seq);
718 
719   p = perfusion_alloc();
720   for (i = roi->ly; i < roi->uy; i++)
721   {
722     p->refit = 0;
723     for (j = roi->lx; j < roi->ux; j++)
724     {
725        if (mask == NULL || IM_CHAR(mask, i, j) > 0)
726        {
727           vec2_y(pos) = i + 0.5;
728           vec2_x(pos) = j + 0.5;
729           if (ct_gamma_fit_pixel(p, pos,imptrs) != -1)
730           {
731              err = ct_emap(p);
732              cbv = ct_rcbv(p);
733 /*
734              if (cbv > 0.0 && sqrt(err/cbv) < err_thresh)
735                 p->refit = 1;
736              else
737              {
738                 p->refit = 0;
739                 gamma_fit_pixel(p, pos,imptrs);
740              }
741 */
742              err = ct_emap(p);
743              rcc = ct_recirc(p);
744              cbv = ct_rcbv(p);
745              if (cbv > 0.0 && sqrt(err/cbv) < err_thresh)
746                 p->refit = 1;
747 
748              if(cbv > 0.0)
749              {
750                 ttp = ct_bolus_m0(p)/cbv;
751                 mtt = sqrt(ttp*ttp - 2.0*ttp*ct_bolus_m0(p)/cbv + ct_bolus_m1(p)/cbv);
752              }
753              else
754              {
755                 ttp = 0.0;
756                 mtt = 0.0;
757              }
758              IM_PIX_SET(perf_CBV, i, j, cbv);
759              IM_PIX_SET(perf_TTP, i, j, ttp - p->n1*p->tscale);
760              IM_PIX_SET(perf_MTT, i, j, mtt);
761              IM_PIX_SET(perf_ERR, i, j, err);
762              IM_PIX_SET(perf_RCC, i, j, rcc);
763           }
764        }
765      }
766      if (i%10 == 0) format("raster %d done \n",i);
767    }
768    perfusion_free(p);
769 
770   return;
771 }
772 
773 Imrect *ct_gfit_get_image(int perf_type)
774 {
775   if (perf_type ==0)
776       return(im_copy(perf_TTP));
777   else if (perf_type ==1)
778       return(im_copy(perf_CBV));
779   else if (perf_type ==2)
780       return(im_copy(perf_MTT));
781   else if (perf_type ==3)
782       return(im_copy(perf_ERR));
783   else if (perf_type ==4)
784       return(im_copy(perf_RCC));
785 
786   return(NULL);
787 }
788 

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