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_perm.c,v $
23 * Date : $Date: 2007/02/15 01:52:29 $
24 * Version : $Revision: 1.14 $
25 * CVS Id : $Id: medStim_perm.c,v 1.14 2007/02/15 01:52:29 paul Exp $
26 *
27 * Notes:
28 *
29 * FIXME: Statics at the top of the file need to be replaced by either a new
30 * structure of by addition to the existing permiability structure
31 * a.lacey@man.ac.uk 26.05.05
32 *
33 *
34 *********
35 */
36
37 #if HAVE_CONFIG_H
38 # include <config.h>
39 #endif
40
41 #include "medStim_perm.h"
42
43 #include <stdio.h>
44 #include <math.h>
45 #include <float.h>
46 #include <tina/sys/sysDef.h>
47 #include <tina/sys/sysPro.h>
48 #include <tina/math/mathDef.h>
49 #include <tina/math/mathPro.h>
50 #include <tina/image/imgDef.h>
51 #include <tina/image/imgPro.h>
52 #include <tina/medical/med_StimDef.h>
53 #include <tina/medical/medStim_alloc.h>
54 #include <tina/medical/medStim_tdata.h>
55
56 #define NP2 4
57 #define FTOL 1.0e-6
58
59 /* FIXME */
60
61 static Permeability *perm_global = NULL;
62 static int imptrlx, imptrux, imptrly, imptruy;
63 static Imrect *perm_TTP=NULL;
64 static Imrect *perm_VP=NULL;
65 static Imrect *perm_VE=NULL;
66 static Imrect *K_trans=NULL;
67 static Imrect *perm_ERR=NULL;
68 static int perm_con=0;
69 static double tr=5.0, prebolus=60.0, r1cont = 4.61, alpha=35.0;
70 static float *plasma_conc = NULL; /* it's a static I know mjs 14/12/05*/
71
72
73 double *tr_get()
74 {
75 return(&tr);
76 }
77
78 double *prebolus_get()
79 {
80 return(&prebolus);
81 }
82
83 double *r1cont_get()
84 {
85 return(&r1cont);
86 }
87
88 double *alpha_get()
89 {
90 return(&alpha);
91 }
92
93 void set_perm_con(int val)
94 {
95 perm_con = val;
96 }
97
98 static float ave_con_pre_bolus(Permeability *p)
99 {
100 float *st_1d = p->st_1d;
101 int k , bolus_time = p->bolus_time;
102 float noisesum = 0.0;
103
104 if (bolus_time-2 <= p->n1) bolus_time = p->n1+3;
105
106 for (k = p->n1; k < bolus_time-2 && k < p->n2; k++)
107 noisesum += st_1d[k];
108
109 return (noisesum/(bolus_time-2-p->n1));
110 }
111
112
113 static float *ther_conc1(Permeability *p)
114 {
115 float *conc1_ther = p->conc1_ther;
116 float *AIF=p->AIF;
117 float tscale = p->tscale/60.0; /*eh why on earth do you want to convert this into minutes? ahhh cos k are in min-1. aaahhhhhhhh*/
118 int n1 = p->n1;
119 int n2 = p->n2;
120 float ve = p->ve;
121 float vp = p->vp;
122 float kin = p->kin;
123 float kout = p->kout;
124 float tt=p->tt, t;
125 int tau;
126 float sum, a ,b;
127 double part;
128 int k, i;
129
130 if (conc1_ther == NULL)
131 if ((conc1_ther = fvector_alloc(n1, n2)) == NULL)
132 {
133 error("ther_conc1: memory alloc error ", warning);
134 return(NULL);
135 }
136
137
138 b = tt-tina_int(tt);
139 a = 1.0-b;
140
141 sum = 0.0;
142
143
144 for(k=n1; k<n2; k++)
145 {
146 t = k-tt;
147 sum=0.0;
148 part=0.0;
149 if (t>=0.0)
150 {
151 for(i=tina_int(tt), tau=0; i<= k; i++, tau++)
152 {
153 part= (a*AIF[tau] + b*AIF[tau+1]); /* aif now always starts at index 0*/
154 sum+= part*exp(-kout*tscale*(t-tau)/ve)*tscale;
155
156 }
157 }
158
159
160 if (t>= 0.0)
161 {
162 part = (a*AIF[tina_int(t)] + b*AIF[tina_int(t+1)]);
163 }
164 else
165 {
166 part = 0.0;
167 }
168 conc1_ther[k]=vp*part + kin*sum;
169 }
170
171 return(conc1_ther);
172
173 }
174
175
176
177
178
179 static float *aif_est(Permeability *p)
180 {
181 int n1=p->n1;
182 int n2=p->n2;
183 int length;
184 int i,k;
185 float *AIF_data = NULL;
186
187
188
189 /*if (( AIF_data = fvector_alloc(n1-(n2-n1), n2+(n2-n1))) == NULL)
190 {
191 error("AIF_data: memory alloc error ", warning);
192 return NULL;
193 }
194
195 for(i=n1-(n2-n1);i<n2+(n2-n1);i++){
196 AIF_data[i]=0.0;
197 }
198
199 AIF_data[n1-1] = 1000.0; */
200 /* i don't think this is consistent, as you're saying the peak occurs at the first point, whereas i think it might be the upstroke occurs at the first point? according to the equations. really need to sort it out */
201 /* oh yes. neil wanted a normalisation to the area anyway didn't he */
202 /*AIF_data[n1] = 5400.0;
203 AIF_data[n1+1] = 1000.0;
204 for(i=n1+2,k=1;i<n2+(n2-n1);i++,k++){
205 AIF_data[i]+=(1000.0-k*5.0);
206 }
207 */
208
209 if (plasma_conc == NULL)
210 {
211 length = n2-n1;
212 if (( AIF_data = fvector_alloc(0, length)) == NULL)
213 {
214 error("AIF_data: memory alloc error ", warning);
215 return NULL;
216 }
217
218 AIF_data[0] = 0.0;
219 AIF_data[1] = 1000.0;
220 AIF_data[2] = 5400.0;
221 AIF_data[3] = 1000.0;
222 for(i=4,k=1;i<n2-1;i++,k++)
223 {
224 AIF_data[i]=(1000.0-k*5.0);
225 }
226 }
227 else
228 {
229 length = n2-n1;
230 if (( AIF_data = fvector_alloc(0, length)) == NULL)
231 {
232 error("AIF_data: memory alloc error ", warning);
233 return NULL;
234 }
235 for (i=0; i<length; i++)
236 {
237 AIF_data[i] = plasma_conc[i];
238 }
239 }
240 return(AIF_data);
241 }
242
243
244 static float *conv_st_to_conc1(Permeability *p)
245 {
246 float *st_1d = p->st_1d;
247 int n1 = p->n1;
248 int n2 = p->n2;
249 float s0 = p->s0_av;
250 float tr = p->tr;
251 float ralpha = p->alpha*PI/180.0;
252 float NH = p->NH;
253 float t10 = p->t10;
254 float *conc;
255 int k;
256 double a,b,c;
257 double R1gd = p->r1cont * 1.0e-6; /* concentration will be parts per million */
258
259 if ((conc = fvector_alloc(n1, n2)) == NULL)
260 {
261 error("conv_st_to_conc: memory alloc error ", warning);
262 return NULL;
263 }
264
265 b = (1.0 - exp(-tr / t10))
266 /(1.0 - cos(ralpha)*exp(-tr/t10));
267 NH = s0/ (b*sin(ralpha));
268 p->NH = NH;
269 for (k = n1; k < n2; k++)
270 {
271 a = (st_1d[k]-s0)/ (NH* sin(ralpha));
272 c=(1.0- (a+b))/(1.0-cos(ralpha)*(a+b));
273 if(c<0.0001) /*mjs just testing 1/12/05 */
274 {
275 c = 0.0001;
276 }
277 conc[k] = (-log(c)/tr - 1.0/ t10)/R1gd;
278 }
279 return(conc);
280 }
281
282
283 static double perm_chi2(int ndim, double *x)
284 /* minimise errors on measured data */
285 {
286 int i;
287 int n1 = perm_global->n1;
288 int n2 = perm_global->n2;
289 float *r1_ther;
290 float *r1 = perm_global->conc1;
291 double diff;
292 double chi=0.0;
293
294 if (x[1]<=0.0) return(FLT_MAX);
295 if (x[2]<=0.0) return(FLT_MAX);
296 if (x[3]<=0.0) return(FLT_MAX);
297
298 perm_global->tt = x[0];
299 perm_global->kin = x[1];
300 if (perm_con==1)
301 {
302 perm_global->kout = 0.0;
303 perm_global->ve = FLT_MAX;
304 }
305 else
306 {
307 perm_global->ve = x[2];
308 if (perm_global->ve > 0.0)
309 {
310 perm_global->kout = x[1];
311 }
312 else
313 {
314 perm_global->kout = 0.0; /* actually. i doubt this matters cos it won't contribute to the
315 result anyway if ve is zero. */
316 /* maybe i should take it back out. check... */
317 }
318
319 }
320 if (perm_con==2) perm_global->vp = 0.0;
321 else perm_global->vp = x[3];
322
323 r1_ther = perm_global->conc1_ther = ther_conc1(perm_global);
324 for(i=n1;i<=n2;i++)
325 {
326 diff = (r1[i]-r1_ther[i]);
327 chi += diff*diff;
328 }
329 return(chi);
330
331 }
332
333
334
335 static void perm_simplex(Permeability *p)
336 {
337 double *x;
338 double *lambda;
339 int ndim = NP2;
340
341 perm_global = p;
342 x = (double *)ralloc(ndim*sizeof(double));
343 x[0] = p->tt;
344 x[1] = p->kin;
345 x[2] = p->ve;
346 x[3] = p->vp;
347
348 lambda = (double *)ralloc(ndim*sizeof(double));
349 lambda[0] = 0.8;
350 lambda[1] = 0.003;
351 lambda[2] = 0.04;
352 lambda[3] = 0.04;
353 p->err=simplexmin2(ndim, x, lambda, perm_chi2, NULL, FTOL, NULL);
354 p->tt=x[0];
355 p->kin=x[1];
356
357 p->kout=x[1];
358 if (perm_con==1)
359 {
360 p->ve = FLT_MAX;
361 p->kout = 0.0; /* this was missing.. not that it makes any difference*/
362 }
363 else
364 {
365 p->ve=x[2];
366 }
367 if (perm_con==2)
368 {
369 p->vp = 0.0;
370 }
371 else
372 {
373 p->vp=x[3];
374 }
375
376 /* doesn't this muck everything up, cos it runs it through perm_chi2 again, after it's done all the optimisation and setting things to zero and stuff.the return value of the simplex is the Chi squared value? */
377 /* p->err = perm_chi2(ndim,x);*/
378 /*p->err = 0.0;*/
379 rfree(x);
380 rfree(lambda);
381 }
382
383
384 void est_plasma_conc(Imrect *mask, Imrect *t1im)
385 {
386 Permeability *p = NULL;
387 Sequence *seq = seq_get_current();
388 List *store = get_seq_start_el(seq);
389 Imrect *im = NULL;
390 Imregion *roi;
391 Vec2 pos;
392 int i, j, a;
393 void ***imptrs=NULL;
394 float err, t_zero=0.0;
395 int n_start, n_end, count=0, length=0;
396 float *plasma_conc_local = NULL;
397 float c, b;
398 float our_time;
399
400 im = (Imrect *)(store->to);
401 roi = im->region;
402
403 if (mask != NULL)
404 {
405 roi = mask->region;
406 if (mask->vtype!=uchar_v)
407 {
408 format("invalid mask image on stack\n");
409 return;
410 }
411 }
412
413 imptrs = seq_slice_init(seq);
414 seq_limits(&n_start, &n_end, &imptrly, &imptruy, &imptrlx, &imptrux);
415
416 if (plasma_conc_local != NULL)
417 fvector_free(plasma_conc_local, n_start);
418 plasma_conc_local = fvector_alloc(n_start, n_end);
419 for (a=n_start; a<n_end; a++)
420 {
421 plasma_conc_local[a]=0.0;
422 }
423
424 p = permeability_alloc();
425 for (i = roi->ly; i < roi->uy; i++)
426 {
427 for (j = roi->lx; j < roi->ux; j++)
428 {
429 if (mask == NULL || IM_CHAR(mask, i, j) > 0 )
430 {
431 vec2_y(pos) = i + 0.5;
432 vec2_x(pos) = j + 0.5;
433 err = 0.0;
434 /* don't really need to do the perm fitting to get the voxe concentrations, but it does
435 repeat an awful lot of the calculations */
436 if (perm_fit_pixel(p, pos,imptrs, t1im) != -1)
437 {
438 for (a=n_start; a<n_end; a++)
439 {
440 plasma_conc_local[a] += p->conc1[a];
441 }
442 count++;
443 }
444 }
445 }
446 }
447 permeability_free(p);
448
449 for (a=n_start; a<n_end; a++)
450 {
451 plasma_conc_local[a] = plasma_conc_local[a]/count;
452 }
453
454 /* now do the permeability fit to the average concentration over time */
455
456 p = permeability_alloc();
457 if ((our_time = seq->dim[3]) == 1.0)
458 {
459 error("perm_fit: timing info missing", warning);
460 /*return(-1); mjs need to fix this to zero.*/
461 }
462
463 p->n1 = n_start;
464 p->n2 = n_end;
465 p->tscale = our_time;
466 p->bolus_time = prebolus/p->tscale;
467 p->conc1 = fvector_alloc(n_start, n_end);
468 for (a=n_start; a<n_end; a++)
469 {
470 p->conc1[a] = plasma_conc_local[a];
471 }
472 if (p->AIF!=NULL)
473 fvector_free(p->AIF, 0);
474 p->AIF = aif_est(p);
475 p->tt = p->bolus_time;
476 p->kin = 0.005;
477 p->kout = 0.005; /* why are these here twice? kin and kout are the same value. mjs 1/12/05 */
478 p->ve = 0.1;
479 p->vp = 0.1;
480 perm_simplex(p);
481 p->conc1_ther = ther_conc1(p);
482
483 t_zero = p->tt;
484
485
486 /* then need to resample the concentration curve in the same way as above*/
487 format("\nt_zero: %f\n", t_zero);
488
489 for (a=n_start; a<n_end; a++)
490 {
491 format("%f\t%f\n", plasma_conc_local[a], p->conc1_ther[a]);
492 }
493
494 permeability_free(p);
495
496 if (plasma_conc != NULL)
497 fvector_free(plasma_conc, 0);
498 length = n_end - n_start;
499 plasma_conc = fvector_alloc(0, length);
500
501 plasma_conc[0] = 0.0;
502 b = t_zero-tina_int(t_zero);
503 c = 1.0-b;
504
505 for (a=tina_int(t_zero+1), count=1; a<n_end-2; a++, count++)
506 {
507 plasma_conc[count]=(c*plasma_conc_local[a] + b*plasma_conc_local[a+1]);
508 }
509
510 /* maybe should change the zeroth point of plasma conc? YES. or not.*/
511
512 for (a=count; a<n_end-1; a++)
513 {
514 plasma_conc[a] = plasma_conc[a-1];
515 }
516
517 for (a=0; a<length; a++)
518 {
519 format("%f\n", plasma_conc[a]);
520 }
521
522 fvector_free(plasma_conc_local, n_start);
523
524 return;
525
526 }
527
528
529 /* mjs 7/11/05 modified perm_fit_pixel to correctly access TE values. */
530 int perm_fit_pixel(Permeability *perm_data, Vec2 v, void ***imptrs, Imrect *t1im)
531 {
532 Sequence *seq= (Sequence *)seq_get_current();
533 List *get_seq_start_el(Sequence *seq);
534 List *store = get_seq_start_el(seq);
535 /* Imrect *im = NULL; */
536 float te = 0.0, *TE_arr = NULL;
537 float our_time;
538
539
540 if (t1im == NULL)
541 return(-1);
542
543 if ((store == NULL) || (store->to == NULL))
544 return(-1);
545 /*im = (Imrect *)(store->to); mjs 7/11/05 removed unnecs. line*/
546
547 if ((our_time = seq->dim[3]) == 1.0)
548 {
549 error("perm_fit: timing info missing", warning);
550 return(-1);
551 }
552
553 perm_data->tscale = our_time;
554
555 if (((float *)prop_get(seq->props, TE_DATA)) == NULL)
556 {
557 error("perm_fit: echo time not found", warning);
558 return(-1);
559 }
560
561 TE_arr = (float*)prop_get(seq->props, TE_DATA);
562 te = TE_arr[seq->offset]; /* mjs 7/11/05 assumes that the TE of the first image
563 is equal to the TEs of the rest */
564
565 seq_limits(&perm_data->n1, &perm_data->n2, &imptrly, &imptruy, &imptrlx, &imptrux);
566
567 perm_data->te = te;
568 perm_data->tr = tr;
569 perm_data->alpha = alpha;
570 perm_data->r1cont = r1cont;
571 perm_data->t10 = im_get_pixf(t1im, tina_int(v.el[1]), tina_int(v.el[0]));
572 perm_data->NH = 100;
573 if (perm_data->st_1d!=NULL)
574 fvector_free(perm_data->st_1d, perm_data->n1);
575 perm_data->st_1d = s_2d_pixel_get(imptrs,v,(Perfusion *)perm_data);
576 perm_data->bolus_time = prebolus/perm_data->tscale;
577 perm_data->s0_av = ave_con_pre_bolus(perm_data);
578 if (perm_data->conc1!=NULL)
579 fvector_free(perm_data->conc1, perm_data->n1);
580 perm_data->conc1 = conv_st_to_conc1(perm_data);
581 if (perm_data->AIF!=NULL)
582 fvector_free(perm_data->AIF, 0);
583 perm_data->AIF = aif_est(perm_data);
584 perm_data->tt = perm_data->bolus_time;
585 perm_data->kin = 0.005;
586 perm_data->kout = 0.005; /* why are these here twice? kin and kout are the same value. mjs 1/12/05 */
587 perm_data->ve = 0.1;
588 perm_data->vp = 0.1;
589 perm_simplex(perm_data);
590 perm_data->conc1_ther = ther_conc1(perm_data);
591
592 return (0);
593 }
594
595
596
597 void pfit_measure_image(double err_thresh, Imrect *mask, Imrect *t1im)
598 {
599 Permeability *p = NULL;
600 Sequence *seq = seq_get_current();
601 List *store = get_seq_start_el(seq);
602 Imrect *im = NULL;
603 Imregion *roi;
604 Vec2 pos;
605 int i, j;
606 void ***imptrs=NULL;
607 float err;
608
609 im = (Imrect *)(store->to);
610 roi = im->region;
611
612 if (mask != NULL)
613 {
614 roi = mask->region;
615 if (mask->vtype!=uchar_v)
616 {
617 format("invalid mask image on stack\n");
618 return;
619 }
620 }
621
622 if (perm_TTP!=NULL) im_free(perm_TTP);
623 perm_TTP = im_alloc(im->height, im->width, roi, float_v);
624 if (perm_VP!=NULL) im_free(perm_VP);
625 perm_VP = im_alloc(im->height, im->width, roi, float_v);
626 if (perm_VE!=NULL) im_free(perm_VE);
627 perm_VE = im_alloc(im->height, im->width, roi, float_v);
628 if (K_trans!=NULL) im_free(K_trans);
629 K_trans = im_alloc(im->height, im->width, roi, float_v);
630 if (perm_ERR!=NULL) im_free(perm_ERR);
631 perm_ERR = im_alloc(im->height, im->width, roi, float_v);
632
633 imptrs = seq_slice_init(seq);
634
635 p = permeability_alloc();
636 for (i = roi->ly; i < roi->uy; i++)
637 {
638 for (j = roi->lx; j < roi->ux; j++)
639 {
640 if (mask == NULL || IM_CHAR(mask, i, j) > 0 )
641 {
642 vec2_y(pos) = i + 0.5;
643 vec2_x(pos) = j + 0.5;
644 err = 0.0;
645 if (perm_fit_pixel(p, pos,imptrs, t1im) != -1)
646 {
647 IM_PIX_SET(perm_TTP, i, j, p->tt*p->tscale);
648 IM_PIX_SET(perm_VP, i, j, p->vp);
649 IM_PIX_SET(perm_VE, i, j, p->ve);
650 IM_PIX_SET(K_trans, i, j, p->kin);
651 IM_PIX_SET(perm_ERR, i, j, p->err);
652 }
653 }
654 }
655 if (i%10 == 0) format("raster %d done \n",i);
656 }
657 permeability_free(p);
658
659 return;
660 }
661
662
663 Imrect *pfit_get_image(int perm_type)
664 {
665 if (perm_type ==0)
666 return(im_copy(perm_TTP));
667 else if (perm_type ==1)
668 return(im_copy(perm_VP));
669 else if (perm_type ==2)
670 return(im_copy(perm_VE));
671 else if (perm_type ==3)
672 return(im_copy(K_trans));
673 else if (perm_type ==4)
674 return(im_copy(perm_ERR));
675 else
676 return(NULL);
677 }
678
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.