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_ct_perm.c,v $
23 * Date : $Date: 2007/03/05 03:32:31 $
24 * Version : $Revision: 1.2 $
25 * CVS Id : $Id: medStim_ct_perm.c,v 1.2 2007/03/05 03:32:31 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_ct_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 #include <tina/medical/medStim_ct_gamma.h>
56
57 #define NP2 4
58 #define FTOL 1.0e-6
59
60 /* FIXME */
61
62 static Permeability *perm_global = NULL;
63 static int imptrlx, imptrux, imptrly, imptruy;
64 static Imrect *perm_TTP=NULL;
65 static Imrect *perm_VP=NULL;
66 static Imrect *perm_VE=NULL;
67 static Imrect *K_trans=NULL;
68 static Imrect *perm_ERR=NULL;
69 static int perm_con=0;
70 static double prebolus=10.0;
71 static float *plasma_conc = NULL; /* it's a static I know mjs 14/12/05*/
72
73
74 double *ct_prebolus_get()
75 {
76 return(&prebolus);
77 }
78
79 void ct_set_perm_con(int val)
80 {
81 perm_con = val;
82 }
83
84 static float ave_con_pre_bolus(Permeability *p)
85 {
86 float *st_1d = p->st_1d;
87 int k , bolus_time = p->bolus_time;
88 float noisesum = 0.0;
89
90 if (bolus_time-2 <= p->n1) bolus_time = p->n1+3;
91
92 for (k = p->n1; k < bolus_time-2 && k < p->n2; k++)
93 noisesum += st_1d[k];
94
95 return (noisesum/(bolus_time-2-p->n1));
96 }
97
98
99 static float *ther_conc1(Permeability *p)
100 {
101 float *conc1_ther = p->conc1_ther;
102 float *AIF=p->AIF;
103 float tscale = p->tscale/60.0;
104 int n1 = p->n1;
105 int n2 = p->n2;
106 float ve = p->ve;
107 float vp = p->vp;
108 float kin = p->kin;
109 float kout = p->kout;
110 float tt=p->tt, t;
111 int tau;
112 float sum, a ,b;
113 double part;
114 int k, i;
115
116 if (conc1_ther == NULL)
117 if ((conc1_ther = fvector_alloc(n1, n2)) == NULL)
118 {
119 error("ther_conc1: memory alloc error ", warning);
120 return(NULL);
121 }
122
123
124 b = tt-tina_int(tt);
125 a = 1.0-b;
126
127 sum = 0.0;
128
129
130 for(k=n1; k<n2; k++)
131 {
132 t = k-tt;
133 sum=0.0;
134 part=0.0;
135 if (t>=0.0)
136 {
137 for(i=tina_int(tt), tau=0; i<= k; i++, tau++)
138 {
139 part= (a*AIF[tau] + b*AIF[tau+1]); /* aif now always starts at index 0*/
140 sum+= part*exp(-kout*tscale*(t-tau)/ve)*tscale;
141
142 }
143 }
144
145
146 if (t>= 0.0)
147 {
148 part = (a*AIF[tina_int(t)] + b*AIF[tina_int(t+1)]);
149 }
150 else
151 {
152 part = 0.0;
153 }
154 conc1_ther[k]=vp*part + kin*sum;
155 }
156
157 return(conc1_ther);
158
159 }
160
161
162
163
164
165 static float *aif_est(Permeability *p)
166 {
167 int n1=p->n1;
168 int n2=p->n2;
169 int length;
170 int i,k;
171 float *AIF_data = NULL;
172
173
174 if (plasma_conc == NULL)
175 {
176 length = n2-n1;
177 if (( AIF_data = fvector_alloc(0, length)) == NULL)
178 {
179 error("AIF_data: memory alloc error ", warning);
180 return NULL;
181 }
182
183 AIF_data[0] = 0.0;
184 AIF_data[1] = 1000.0;
185 AIF_data[2] = 5400.0;
186 AIF_data[3] = 1000.0;
187 for(i=4,k=1;i<n2-1;i++,k++)
188 {
189 AIF_data[i]=(1000.0-k*5.0);
190 }
191 }
192 else
193 {
194 length = n2-n1;
195 if (( AIF_data = fvector_alloc(0, length)) == NULL)
196 {
197 error("AIF_data: memory alloc error ", warning);
198 return NULL;
199 }
200 for (i=0; i<length; i++)
201 {
202 AIF_data[i] = plasma_conc[i];
203 }
204 }
205 return(AIF_data);
206 }
207
208 static float *conv_st_to_conc1(Permeability *p)
209 {
210 float *st_1d = p->st_1d;
211 int n1 = p->n1;
212 int n2 = p->n2;
213 float s0 = p->s0_av;
214 float *conc;
215 int k;
216
217 if ((conc = fvector_alloc(n1, n2)) == NULL)
218 {
219 error("conv_st_to_conc: memory alloc error ", warning);
220 return NULL;
221 }
222
223 for (k = n1; k < n2; k++)
224 {
225 conc[k] = st_1d[k]-s0;
226 }
227 return(conc);
228 }
229
230
231 static double perm_chi2(int ndim, double *x)
232 /* minimise errors on measured data */
233 {
234 int i;
235 int n1 = perm_global->n1;
236 int n2 = perm_global->n2;
237 float *r1_ther;
238 float *r1 = perm_global->conc1;
239 double diff;
240 double chi=0.0;
241
242 if (x[1]<=0.0) return(FLT_MAX);
243 if (x[2]<=0.0) return(FLT_MAX);
244 if (x[3]<=0.0) return(FLT_MAX);
245
246 perm_global->tt = x[0];
247 perm_global->kin = x[1];
248 if (perm_con==1)
249 {
250 perm_global->kout = 0.0;
251 perm_global->ve = FLT_MAX;
252 }
253 else
254 {
255 perm_global->ve = x[2];
256 if (perm_global->ve > 0.0)
257 {
258 perm_global->kout = x[1];
259 }
260 else
261 {
262 perm_global->kout = 0.0; /* actually. i doubt this matters cos it won't contribute to the
263 result anyway if ve is zero. */
264 /* maybe i should take it back out. check... */
265 }
266
267 }
268 if (perm_con==2) perm_global->vp = 0.0;
269 else perm_global->vp = x[3];
270
271 r1_ther = perm_global->conc1_ther = ther_conc1(perm_global);
272 for(i=n1;i<=n2;i++)
273 {
274 diff = (r1[i]-r1_ther[i]);
275 chi += diff*diff;
276 }
277 return(chi);
278
279 }
280
281
282
283 static void perm_simplex(Permeability *p)
284 {
285 double *x;
286 double *lambda;
287 int ndim = NP2;
288
289 perm_global = p;
290 x = (double *)ralloc(ndim*sizeof(double));
291 x[0] = p->tt;
292 x[1] = p->kin;
293 x[2] = p->ve;
294 x[3] = p->vp;
295
296 lambda = (double *)ralloc(ndim*sizeof(double));
297 lambda[0] = 0.8;
298 lambda[1] = 0.003;
299 lambda[2] = 0.04;
300 lambda[3] = 0.04;
301 p->err=simplexmin2(ndim, x, lambda, perm_chi2, NULL, FTOL, NULL);
302 p->tt=x[0];
303 p->kin=x[1];
304
305 p->kout=x[1];
306 if (perm_con==1)
307 {
308 p->ve = FLT_MAX;
309 p->kout = 0.0; /* this was missing.. not that it makes any difference*/
310 }
311 else
312 {
313 p->ve=x[2];
314 }
315 if (perm_con==2)
316 {
317 p->vp = 0.0;
318 }
319 else
320 {
321 p->vp=x[3];
322 }
323
324 /* 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? */
325 /* p->err = perm_chi2(ndim,x);*/
326 /*p->err = 0.0;*/
327 rfree(x);
328 rfree(lambda);
329 }
330
331 void ct_est_plasma_conc(Imrect *mask)
332 {
333 Perfusion *p = NULL;
334 Sequence *seq = seq_get_current();
335 List *store = get_seq_start_el(seq);
336 Imrect *im = NULL;
337 Imregion *roi;
338 void ***imptrs=NULL;
339 int n_start, n_end, count=0, length=0, a;
340 float *plasma_conc_local = NULL, c, b, t_zero=0.0;
341
342 im = (Imrect *)(store->to);
343 roi = im->region;
344
345 if (mask != NULL)
346 {
347 roi = mask->region;
348 if (mask->vtype!=uchar_v)
349 {
350 format("invalid mask image on stack\n");
351 return;
352 }
353 }
354
355 imptrs = seq_slice_init(seq);
356 seq_limits(&n_start, &n_end, &imptrly, &imptruy, &imptrlx, &imptrux);
357
358 if (plasma_conc_local != NULL)
359 fvector_free(plasma_conc_local, n_start);
360 plasma_conc_local = fvector_alloc(n_start, n_end);
361 for (a=n_start; a<n_end; a++)
362 {
363 plasma_conc_local[a]=0.0;
364 }
365
366 p = perfusion_alloc();
367
368 t_zero = ct_gamma_fit_region(p, mask);
369
370 for (a=n_start; a<n_end; a++)
371 {
372 plasma_conc_local[a]=p->r2[a];
373 format("%f\n", plasma_conc_local[a]);
374 }
375 format("\n");
376 perfusion_free(p);
377
378 if (plasma_conc != NULL)
379 fvector_free(plasma_conc, 0);
380 length = n_end - n_start;
381 plasma_conc = fvector_alloc(0, length);
382
383 plasma_conc[0] = 0.0;
384 b = t_zero-tina_int(t_zero);
385 c = 1.0-b;
386
387 for (a=tina_int(t_zero+1), count=1; a<n_end-2; a++, count++)
388 {
389 plasma_conc[count]=(c*plasma_conc_local[a] + b*plasma_conc_local[a+1]);
390 format("%f\n", plasma_conc[count]);
391 }
392
393 /* maybe should change the zeroth point of plasma conc? YES. or not.*/
394
395 for (a=count; a<n_end-1; a++)
396 {
397 plasma_conc[a] = plasma_conc[a-1];
398 format("%f\n", plasma_conc[count]);
399 }
400
401
402 fvector_free(plasma_conc_local, n_start);
403
404 return;
405
406 }
407
408
409 /* mjs 7/11/05 modified perm_fit_pixel to correctly access TE values. */
410 int ct_perm_fit_pixel(Permeability *perm_data, Vec2 v, void ***imptrs)
411 {
412 Sequence *seq= (Sequence *)seq_get_current();
413 List *get_seq_start_el(Sequence *seq);
414 List *store = get_seq_start_el(seq);
415 float our_time;
416
417 if ((store == NULL) || (store->to == NULL))
418 return(-1);
419 /*im = (Imrect *)(store->to); mjs 7/11/05 removed unnecs. line*/
420
421 if ((our_time = seq->dim[3]) == 0.0)
422 {
423 error("perm_fit: timing info missing", warning);
424 return(-1);
425 }
426
427 perm_data->tscale = our_time;
428
429 seq_limits(&perm_data->n1, &perm_data->n2, &imptrly, &imptruy, &imptrlx, &imptrux);
430
431 if (perm_data->st_1d!=NULL)
432 fvector_free(perm_data->st_1d, perm_data->n1);
433 perm_data->st_1d = s_2d_pixel_get(imptrs,v,(Perfusion *)perm_data);
434 perm_data->bolus_time = prebolus/perm_data->tscale;
435 perm_data->s0_av = ave_con_pre_bolus(perm_data);
436 if (perm_data->conc1!=NULL)
437 fvector_free(perm_data->conc1, perm_data->n1);
438 perm_data->conc1 = conv_st_to_conc1(perm_data);
439 if (perm_data->AIF!=NULL)
440 fvector_free(perm_data->AIF, 0);
441 perm_data->AIF = aif_est(perm_data);
442 perm_data->tt = perm_data->bolus_time;
443 perm_data->kin = 0.005;
444 perm_data->kout = 0.005; /* why are these here twice? kin and kout are the same value. mjs 1/12/05 */
445 perm_data->ve = 0.1;
446 perm_data->vp = 0.1;
447 perm_simplex(perm_data);
448 perm_data->conc1_ther = ther_conc1(perm_data);
449
450 return (0);
451 }
452
453
454
455
456 void ct_kfpfit_measure_image(double err_thresh, Imrect *mask)
457 {
458 Permeability *p = NULL;
459 Sequence *seq = seq_get_current();
460 List *store = get_seq_start_el(seq);
461 Imrect *im = NULL;
462 Imregion *roi;
463 Vec2 pos;
464 int i, j;
465 void ***imptrs=NULL;
466 float err;
467
468 im = (Imrect *)(store->to);
469 roi = im->region;
470
471 /*if (mask != NULL)
472 {
473 roi = mask->region;
474
475 if (mask->vtype!=uchar_v)
476 {
477 format("invalid mask image on stack\n");
478 return;
479 }
480 }*/
481
482 if (perm_TTP!=NULL) im_free(perm_TTP);
483 perm_TTP = im_alloc(im->height, im->width, roi, float_v);
484 if (perm_VP!=NULL) im_free(perm_VP);
485 perm_VP = im_alloc(im->height, im->width, roi, float_v);
486 if (perm_VE!=NULL) im_free(perm_VE);
487 perm_VE = im_alloc(im->height, im->width, roi, float_v);
488 if (K_trans!=NULL) im_free(K_trans);
489 K_trans = im_alloc(im->height, im->width, roi, float_v);
490 if (perm_ERR!=NULL) im_free(perm_ERR);
491 perm_ERR = im_alloc(im->height, im->width, roi, float_v);
492
493
494 if (mask != NULL)
495 {
496 roi = mask->region;
497
498 /*if (mask->vtype!=uchar_v)
499 {
500 format("invalid mask image on stack\n");
501 return;
502 }*/
503 }
504
505 imptrs = seq_slice_init(seq);
506
507 p = permeability_alloc();
508 for (i = roi->ly; i < roi->uy; i++)
509 {
510 for (j = roi->lx; j < roi->ux; j++)
511 {
512 if (mask == NULL || IM_CHAR(mask, i, j) > 0 )
513 {
514 vec2_y(pos) = i + 0.5;
515 vec2_x(pos) = j + 0.5;
516 err = 0.0;
517 if (ct_perm_fit_pixel(p, pos,imptrs) != -1)
518 {
519 IM_PIX_SET(perm_TTP, i, j, p->tt*p->tscale);
520 IM_PIX_SET(perm_VP, i, j, p->vp);
521 IM_PIX_SET(perm_VE, i, j, p->ve);
522 IM_PIX_SET(K_trans, i, j, p->kin);
523 IM_PIX_SET(perm_ERR, i, j, p->err);
524 }
525 }
526 }
527 if (i%10 == 0) format("raster %d done \n",i);
528 }
529 permeability_free(p);
530
531 return;
532 }
533
534
535
536
537 Imrect *ct_kfpfit_get_image(int perm_type)
538 {
539 if (perm_type ==0)
540 return(im_copy(perm_TTP));
541 else if (perm_type ==1)
542 return(im_copy(perm_VP));
543 else if (perm_type ==2)
544 return(im_copy(perm_VE));
545 else if (perm_type ==3)
546 return(im_copy(K_trans));
547 else if (perm_type ==4)
548 return(im_copy(perm_ERR));
549 else
550 return(NULL);
551 }
552
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.