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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.