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

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

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