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

Linux Cross Reference
Tina5/tina-libs/tina/medical/medStim_ct_perm.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_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 

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