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