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

Linux Cross Reference
Tina4/src/tools/nmr-segment/stim_results.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /*
  2   stim_results.c
  3 
  4   n.thacker
  5   a.lacey   3.2.98
  6 */
  7 
  8 #include <stdio.h>
  9 #include <math.h>
 10 #include <tina/sys.h>
 11 #include <tina/sysfuncs.h>
 12 #include <tina/math.h>
 13 #include <tina/mathfuncs.h>
 14 #include <tina/vision.h>
 15 #include <tina/visionfuncs.h>
 16 #include <tina/tv.h>
 17 #include <tina/tvfuncs.h>
 18 #include <tina/draw.h>
 19 #include <tina/drawfuncs.h>
 20 #include <tina/seqoral.h>
 21 #include <tina/seqpro.h>
 22 #include <tina/stimdef.h>
 23 
 24 
 25 /*
 26   Prototypes
 27 */
 28 
 29 extern void init_interp(int nblx, int nbux, int nbly, int nbuy,
 30                         int nblz, int nbuz);                        
 31 extern void ***coreg_limits(int *lz, int *uz, int *ly,
 32                             int *uy, int *lx, int *ux);
 33 extern void ***coreg_slice_init(Vec3 *v);
 34 extern float   nearest_pixel(void ***rasptrs, Vec3 post);
 35 extern float  *get_stimulus();
 36 
 37 
 38 static void ***imptrs = NULL;
 39 static Imrect *perf_TTP=NULL;
 40 static Imrect *perf_CBV=NULL;
 41 static Imrect *perf_MTT=NULL;
 42 static Imrect *perf_ERR=NULL;
 43 static Imrect *perf_RCC=NULL;
 44 
 45 void  plot_stimfunc(void)
 46 {
 47   Sequence *seq = NULL;
 48   Pl_flow  *plot_data = pl_flow_alloc();
 49   List   *store = (List *)get_start_el();
 50   Tv       *tv = (Tv *)imcalc_graph_tv_get();
 51   Imrect   *im;
 52   float    *time, *flow;
 53   int       imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
 54   int       length,  k;
 55   int       type;
 56   float    *row;
 57   float    *stim;
 58 
 59   if (tv == NULL || tv->tv_screen == NULL)
 60     return;
 61   if ((store == NULL) || (store->to == NULL))
 62     return;
 63   if ((stim=get_stimulus()) == NULL)
 64     return;
 65 
 66   seq = (Sequence *)get_seq_ptr();
 67   im = (Imrect *)(store->to);
 68 
 69   coreg_slice_init((Vec3 *)NULL);
 70   imptrs = coreg_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx, 
 71                         &imptrux);
 72   
 73   length = imptruz-imptrlz;    
 74   flow = (float *)fvector_alloc(0, length);
 75   time = (float *)fvector_alloc(0, length); 
 76 
 77   for (k = 0; k < length; k++)
 78     {
 79       if (!isnan((double)stim[k]))
 80         flow[k] = stim[k];
 81       else
 82         flow[k] = 0.0;
 83       time[k] = (float)k;
 84     }
 85 
 86   tv_erase(tv);
 87   plot(PL_INIT, PL_TV, tv,
 88        PL_AXIS_COLOR, black,
 89        PL_X_TEXT, true, 
 90        PL_Y_TEXT, true,
 91        PL_TITLE, "Stim Signal v Frame",
 92        PL_STYLE, PL_CROSS,
 93        PL_COLOR, red,
 94        PL_GRAPH_DATA, length, time, flow,
 95        PL_PLOT,
 96        NULL);
 97 
 98   plot_data->mask = NULL;
 99   plot_data->y1 = flow;
100   plot_data->x = time;
101   plot_data->y1_n2 = plot_data->x_n2 = length;
102   plot_data->y1_n1 = plot_data->x_n1 = 0;
103   set_pl_flow(plot_data);
104   
105   return;
106 }
107 
108 
109 void  plot_flow(void)
110 {
111   Sequence *seq = NULL;
112   Pl_flow  *plot_data = pl_flow_alloc();
113   Vec3     *scale = NULL;
114   List   *store = (List *)get_start_el();
115   Imrect   *mask = NULL, *fmask = NULL, *im;
116   Tv       *tv = (Tv *)imcalc_graph_tv_get();
117   float    *time, *flow;
118   int       imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
119   int       length,  i, j, k;
120   int       type;
121 
122 
123   if (tv == NULL || tv->tv_screen == NULL)
124     return;
125   if ((store == NULL) || (store->to == NULL))
126     return;
127 
128   seq = (Sequence *)get_seq_ptr();
129   im = (Imrect *)(store->to);
130 
131   if ((time = (float *)prop_get(im->props, DYNSTIME)) == NULL)
132     {
133       error("plot_flow: timing info missing", warning);
134       return;
135     }
136     
137   if ((scale = (Vec3 *)prop_get(im->props, VOXELS)) == NULL)
138     {
139       error("plot_flow: scale info missing", warning);
140       return;
141     }
142                  
143   if (stack_check_types(IMRECT, NULL) == false)
144     {
145       error("plot_flow: wrong type on stack", warning);
146       return;
147     }
148   mask = (Imrect *) stack_pop(&type);
149   fmask = im_copy(mask);
150 
151   coreg_slice_init((Vec3 *)NULL);
152   imptrs = coreg_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx, 
153                         &imptrux);
154   init_interp(imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz);
155   
156   length = imptruz-imptrlz;    
157   flow = (float *)fvector_alloc(0, length);
158  
159   for (k = 0; k < length; k++)
160     {
161       flow[k] = 0.0;
162 
163       for (i = fmask->region->ly; i < fmask->region->uy; i++)
164         {
165           for (j = fmask->region->lx; j < fmask->region->ux; j++)
166             {
167               if (IM_CHAR(fmask, i, j) > 0)
168                  flow[k] += nearest_pixel(imptrs, vec3((double)j, (double)i, (double)k+imptrlz));
169             }
170         }
171       flow[k] *= ((scale->el[0])*(scale->el[1])/100);
172     }
173 
174 
175   tv_erase(tv);
176   plot(PL_INIT, PL_TV, tv,
177        PL_AXIS_COLOR, black,
178        PL_X_TEXT, true,
179        PL_Y_TEXT, true,
180        PL_TITLE, "Signal v Time",
181        PL_STYLE, PL_CROSS,
182        PL_COLOR, red,
183        PL_GRAPH_DATA, length, time, flow,
184        PL_PLOT,
185        NULL);
186 
187   stack_push((void *)mask, IMRECT, im_free);
188 
189   plot_data->mask = fmask;
190   plot_data->y1 = flow;
191   plot_data->x = fvector_copy((void *)time, 0, length);
192   plot_data->y1_n2 = plot_data->x_n2 = length;
193   plot_data->y1_n1 = plot_data->x_n1 = 0;
194   set_pl_flow(plot_data);
195   
196   return;
197 }
198 
199 
200 void   plot_ct_flow(void)
201 {
202   float   *s_2d_mask_mean(void ***imptrs, Imrect *mask, Perfusion *p);
203   float    *conv_st_to_r2(Perfusion *p); 
204   float     ave_s0_pre_bolus(Perfusion *p);
205   int       find_bolus_time(Perfusion *p);
206   Perfusion *p;
207   Pl_flow  *plot_data = pl_flow_alloc();
208   Sequence *seq = NULL;
209   Vec3     *scale = NULL;
210   List   *store = (List *)get_start_el();
211   Imrect   *mask = NULL, *fmask = NULL, *im = NULL;
212   Tv       *tv = (Tv *)imcalc_graph_tv_get();
213   float    *TE;
214   float    *flow, *r2, *time;
215   int      *stime;
216   int       imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
217   int       length,  i, j, k;
218   int       type;
219 
220  
221   if (tv == NULL || tv->tv_screen == NULL)
222     return;
223   if ((store == NULL) || (store->to == NULL))
224     return;
225 
226   seq = (Sequence *)get_seq_ptr();
227   im = (Imrect *)(store->to);
228   if ((time = (float *)prop_get(im->props, DYNSTIME)) == NULL)
229     {
230       error("plot_ct_flow: timing info missing", warning);
231       return;
232     }
233 
234   if ((TE = (float *)prop_get(im->props, TE_DATA)) == NULL)
235     {
236       error("plot_flow: echo time info missing", warning);
237       return;
238     }
239 
240   if (stack_check_types(IMRECT, NULL) == false)
241     {
242       error("plot_ctflow_proc: wrong type on stack", warning);
243       return;
244     }
245 
246   mask = (Imrect *) stack_pop(&type);
247   fmask = im_copy(mask);
248 
249   coreg_slice_init((Vec3 *)NULL);
250   imptrs = coreg_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx, &imptrux); 
251   init_interp(imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz);
252   length = imptruz-imptrlz;
253  
254   p = perfusion_alloc();
255   p->n1 = imptrlz;
256   p->n2 = imptruz;
257   p->te = *TE;
258   p->st_1d = s_2d_mask_mean(imptrs, fmask, p);
259   p->bolus_time = find_bolus_time(p);
260   p->s0_av = ave_s0_pre_bolus(p);
261   p->tscale = time[2]-time[1];
262   r2 = conv_st_to_r2(p);
263 
264   tv_erase(tv);
265   plot(PL_INIT, PL_TV, tv,
266        PL_AXIS_COLOR, black,
267        PL_X_TEXT, true,
268        PL_Y_TEXT, true,
269        PL_TITLE, "Concentration v Time",
270        PL_STYLE, PL_CROSS,
271        PL_COLOR, red,
272        PL_GRAPH_DATA, length, time, r2,
273        PL_PLOT,
274        NULL);
275 
276   stack_push((void *)mask, IMRECT, im_free);
277   perfusion_free(p);
278 
279   plot_data->mask = fmask;
280   plot_data->y1 = &r2[imptrlz];
281   plot_data->x = fvector_copy((void *)time, 0, length);
282   plot_data->y1_n2 = plot_data->x_n2 = length;
283   plot_data->y1_n1 = plot_data->x_n1 = 0;
284   set_pl_flow(plot_data);
285 
286 
287   return;
288 }
289 
290 void    gfit_measure_image(double err_thresh)
291 {
292   int        gamma_fit_pixel(Perfusion *p, Vec2 pos, void ***imptrs);
293   Perfusion *perfusion_alloc(void);
294   Perfusion *p = NULL;
295   List    *store = (List *)get_start_el();
296   Imrect    *im = NULL, *t0_im = NULL;
297   Imrect    *mask = NULL;
298   Imregion  *roi;
299   float      rcbv(Perfusion *p), emap(Perfusion *p), recirc(Perfusion *p);
300   float      bolus_m0(Perfusion *p), bolus_m1(Perfusion *p);
301   float      cbv, rcc, mtt, err, ttp;
302   Vec2       pos;
303   int        i, j, type;
304   void    ***imptrs=NULL;
305   
306   im = (Imrect *)(store->to);
307   roi = im->region;
308 
309   if (perf_TTP!=NULL) im_free(perf_TTP);
310   perf_TTP = im_alloc(im->height, im->width, roi, float_v);
311   if (perf_CBV!=NULL) im_free(perf_CBV);
312   perf_CBV = im_alloc(im->height, im->width, roi, float_v);
313   if (perf_MTT!=NULL) im_free(perf_MTT);
314   perf_MTT = im_alloc(im->height, im->width, roi, float_v);
315   if (perf_ERR!=NULL) im_free(perf_ERR);
316   perf_ERR = im_alloc(im->height, im->width, roi, float_v);
317   if (perf_RCC!=NULL) im_free(perf_RCC);
318   perf_RCC = im_alloc(im->height, im->width, roi, float_v);
319 
320   if (stack_check_types(IMRECT, NULL) != false)
321   {
322     mask = (Imrect *) stack_pop(&type);
323     roi = mask->region; 
324   }
325   imptrs = coreg_slice_init((Vec3 *)NULL);
326 
327   p = perfusion_alloc();
328   for (i = roi->ly; i < roi->uy; i++)
329   {
330     p->refit = 0;
331     for (j = roi->lx; j < roi->ux; j++)
332     {
333        if (mask == NULL || IM_CHAR(mask, i, j) > 0)
334        {
335           vec2_y(pos) = i + 0.5;
336           vec2_x(pos) = j + 0.5;
337           if (gamma_fit_pixel(p, pos,imptrs) != -1)
338           {
339              err = emap(p); 
340              cbv = rcbv(p);
341 /*
342              if (cbv > 0.0 && sqrt(err/cbv) < err_thresh)
343                 p->refit = 1;
344              else
345              {
346                 p->refit = 0;
347                 gamma_fit_pixel(p, pos,imptrs);
348              }
349 */
350              err = emap(p); 
351              rcc = recirc(p); 
352              cbv = rcbv(p);
353              if (cbv > 0.0 && sqrt(err/cbv) < err_thresh)
354                 p->refit = 1;
355 
356              if(cbv > 0.0)
357              {
358                 ttp = bolus_m0(p)/cbv;
359                 mtt = sqrt(ttp*ttp - 2.0*ttp*bolus_m0(p)/cbv + bolus_m1(p)/cbv);
360              }
361              else
362              {
363                 ttp = 0.0;
364                 mtt = 0.0;
365              }
366              IM_PIX_SET(perf_CBV, i, j, cbv);
367              IM_PIX_SET(perf_TTP, i, j, ttp - p->n1*p->tscale);
368              IM_PIX_SET(perf_MTT, i, j, mtt);
369              IM_PIX_SET(perf_ERR, i, j, err);
370              IM_PIX_SET(perf_RCC, i, j, rcc);
371           }
372        }
373      }
374      if (i%10 == 0) format("raster %d done \n",i);
375    }
376    perfusion_free(p);
377 
378   return;
379 }
380 
381 Imrect *gfit_get_image(int perf_type)
382 {
383   Perfusion *p = NULL;
384   Imrect    *im = NULL, *par_im = NULL;
385   Imregion  *roi;
386   int        i, j, type;
387 
388   if (perf_type ==0)
389       return(im_copy(perf_TTP));
390   else if (perf_type ==1)
391       return(im_copy(perf_CBV));
392   else if (perf_type ==2)
393       return(im_copy(perf_MTT));
394   else if (perf_type ==3)
395       return(im_copy(perf_ERR));
396   else if (perf_type ==4)
397       return(im_copy(perf_RCC));
398 }
399 
400 

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