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

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

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