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

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

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