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

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

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