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

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

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