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

Linux Cross Reference
Tina6/tina-tools/tinatool/tlmedical/tlmedFlow_tool.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    : 
 27  * Date    : 
 28  * Version :  
 29  * CVS Id  :  
 30  *
 31  * Author  : mjs
 32  *
 33  * Notes :
 34  *
 35  *********
 36 */
 37 
 38 /**
 39  * @file
 40  * @brief Functions to perform Cerebral Blood Flow calculations
 41  *
 42  * The functions in this file perform Cerebral Blood Flow calculations (as detailed in
 43  * Tina Memo 2002-001) for a volume of images, to produce normalised square-rooted cerebral blood 
 44  * volume (sqrt RCBV), 3D median filtered time of arrival (TTM) , net mean transit time (NMTT),  
 45  * and net cerebral blood flow (NCBF).
 46  * In addition, there is also the capability to produce regional estimates of these parameters
 47  * based on the Talairach stereotaxic atlas.
 48  */
 49 
 50 
 51 #include "tlmedFlow_tool.h"
 52 
 53 #if HAVE_CONFIG_H
 54   #include <config.h>
 55 #endif
 56 
 57 #include <stdio.h>
 58 #include <sys/param.h>
 59 #include <string.h>
 60 #include <math.h>
 61 
 62 #include <tina/sys/sysDef.h>
 63 #include <tina/sys/sysPro.h>
 64 #include <tina/math/mathDef.h>
 65 #include <tina/math/mathPro.h>
 66 #include <tina/image/imgDef.h>
 67 #include <tina/image/imgPro.h>
 68 #include <tina/geometry/geomDef.h>
 69 #include <tina/geometry/geomPro.h>
 70 #include <tina/medical/medDef.h>
 71 #include <tina/medical/medPro.h>
 72 #include <tinatool/draw/drawDef.h>
 73 #include <tinatool/draw/drawPro.h>
 74 #include <tinatool/wdgts/wdgtsDef.h>
 75 #include <tinatool/wdgts/wdgtsPro.h>
 76 #include <tinatool/tlvision/tlvisEdge_epi_mouse.h>
 77 #include <tinatool/tlbase/tlbaseDef.h>
 78 #include <tinatool/tlbase/tlbasePro.h>
 79 #include <tinatool/tlmedical/tlmedPro.h>
 80 #include <tinatool/tlmedical/tlmedDef.h>
 81 
 82 static Sequence *cbv_seq = NULL;
 83 static Sequence *ttm_seq = NULL;
 84 static Sequence *flow_seq = NULL;
 85 static Sequence *nmtt_seq = NULL;
 86 static Sequence *smoothed_seq = NULL;
 87 static Sequence *nmtt_mask_seq = NULL;
 88 
 89 static shistogram **shist_vec = NULL, **shist_vec_left = NULL, **shist_vec_right = NULL;
 90 static int hist_no = 1;
 91 static float *imptrs = NULL;
 92 static void *ptalno = NULL;
 93 static int hist_choice = 0;
 94 
 95 
 96 /** @brief Assigns the current sequence as the CBV volume
 97  *  @param void
 98  *  @return void
 99  *
100  *  Copies the current sequence, as shown in the sequence tool to the 
101  *  CBV sequence. If a CBV sequence already exists, this is deleted first.
102  *  The images in the current sequence should be the "raw" CBV images, as would be obtained
103  *  using the TINA NMR analysis tool.
104  */
105 
106 static void get_cbv_vol_proc(void)
107 {
108   Sequence *seq = NULL;
109   List *lptr = NULL;
110 
111   if ((seq = seq_get_current()) == NULL)
112     {
113       format("No sequence available\n");
114       return;
115     }
116  
117   if (cbv_seq != NULL)
118     {
119       seq_rm(cbv_seq);
120       cbv_seq = NULL;
121     }
122   cbv_seq = seq_copy(seq);
123 }
124 
125 
126 /** @brief Assigns the current sequence to be the TTM sequence 
127  *  @param void
128  *  @return void
129  *
130  *  Copies the current sequence, as shown in the sequence tool to the 
131  *  TTM sequence. If a TTM sequence already exists, the TTM, NMTT, and 3D median smoothed TTM
132  *  sequences are all deleted first.
133  *  The images in the current sequence should be the raw TTM images as would be obtained using the 
134  *  TINA NMR analysis software
135  */
136 static void get_ttm_vol_proc(void)
137 {
138   Sequence *seq = NULL;
139   List *lptr = NULL;
140   
141   if ((seq = seq_get_current()) == NULL)
142     {
143       format("No sequence available\n");
144       return;
145     }
146 
147   if (ttm_seq != NULL)
148     {
149       seq_rm(ttm_seq);
150       seq_rm(nmtt_seq);    
151       seq_rm(smoothed_seq);
152       ttm_seq = NULL;
153       nmtt_seq = NULL;
154       smoothed_seq = NULL;
155     }
156 
157 
158   ttm_seq = seq_copy(seq);
159 }
160 
161 /** @brief Processes the CBV sequence to obtain the normalised square-rooted RCBV maps
162  *  @param void
163  *  @return void
164  *
165  *  The function initially replaces each image in the sequence with its square-root, and then 
166  *  produces a maximum intensity projection (MIP) of the images. The MIP undergoes tangential
167  *  smoothing, and the maximum voxel value in the image found. This voxel is assumed to be 
168  *  100% blood, so the square-root CBV volume is normalised to this voxel value to give % maps of
169  *  blood volume.
170  *  The MIP is pushed onto the stack (visible on Imcalc TV) and the current sequence is set to the 
171  *  sqrt RCBV maps (visible on the Sequence Tv).
172  */
173 static void process_cbv_proc(void)
174 {
175   List *lptr = NULL;
176   Imrect *mip_im = NULL;
177   Imrect *mip_im2 = NULL;
178   double max_val = 0.0;
179   int x, y;
180 
181   if (cbv_seq == NULL)
182     {
183       format("No CBV sequence\n");
184       return;
185     }
186 
187   for (lptr = cbv_seq->start; lptr != NULL; lptr = lptr->next)
188     {
189       Imrect *seq_im = NULL, *seq_im2 = NULL;
190      
191       seq_im = lptr->to;
192       seq_im2 = im_sqrt(seq_im);
193       lptr->to = seq_im2;
194       im_free(seq_im);
195 
196       if (lptr == cbv_seq->start)
197         {
198           mip_im = im_maxsel(seq_im2, seq_im2);
199         }
200       else
201         {
202           mip_im2 = im_maxsel(mip_im, seq_im2);
203           im_free(mip_im);
204           mip_im = im_copy(mip_im2);
205           im_free(mip_im2);
206           mip_im2 = NULL;
207         }
208     }
209   
210   mip_im2 = im_tsmooth(mip_im);
211   im_free(mip_im);
212   max_val = im_locate_max(mip_im2, &y, &x);
213   stack_push((void *)mip_im2, IMRECT, im_free);
214   format("max value %5.3f \n",max_val);
215 
216   if (max_val > 2000.0) 
217     {
218       format("Warning: max sqrt(CBV) value is greater than 2000,\n check there are no outliers in the data\n");
219     }
220   
221   
222   for (lptr = cbv_seq->start; lptr != NULL; lptr = lptr->next)
223     {
224       Imrect *seq_im = NULL, *seq_im2 = NULL;
225       
226       seq_im = lptr->to;
227       seq_im2 = im_times((1.0/max_val), seq_im);
228       lptr->to = seq_im2;
229       im_free(seq_im);
230       
231     }
232   seq_set_current(cbv_seq);
233 
234  
235 }
236 
237 /** @brief Processes the TTM maps to produce the NMTT maps
238  *  @param void
239  *  @return void
240  *
241  *  Initially the TTM images are thresholded at 0.0 to remove spurious negative values.
242  *  The TTM maps are then 3d median smoothed, using the slice of interest and the adjacent two slices.
243  *  The smoothed images are maintained in a separate sequence (smoothed_seq).
244  *  The NMTT maps are then produced using the spatial differentiation in 3D of the smoothed 
245  *  TTM maps. Note that the NMTT volume is 2 slices shorter than the TTM volumes.
246  */
247 static void process_ttm_proc(void)
248 {
249   /* 3d median smooths entire volume 
250      nb, need to maintain two sequences here, and delete original
251      at end of conversion
252      calculate 3d spatial derivative.nb, two less (? or is it 4) images
253      in sequence than before. 
254   */
255         
256   List *lptr = NULL;
257   List *smoothed_images = NULL, *nmtt_images = NULL;
258   double mtt_scale = 0.0;
259 
260   if (ttm_seq == NULL)
261     {
262       format("No TTM sequence\n");
263       return;
264     }
265 
266   /* threshold images first */
267 
268   for (lptr = ttm_seq->start; lptr != NULL; lptr = lptr->next)
269     {
270       Imrect *thresh_im = NULL, *temp = NULL;
271 
272       temp = lptr->to;
273       thresh_im = im_thresh(0.0, temp);
274       lptr->to = thresh_im;
275       im_free(temp); 
276 
277     }
278     
279 
280   /* 3d median smooth - produces a new sequence */
281   
282   for (lptr = ttm_seq->start; lptr != NULL; lptr = lptr->next)
283     {
284       Imrect *im_new = NULL, *im_bef = NULL, *im_aft = NULL, *im_cen = NULL;
285       List *smoothed_link = NULL;
286 
287       im_cen = lptr->to;
288       /* The following if statement deals with the cases at the start and end of the sequence
289          where there are no below and above slices respectively. In these cases the slice of
290          interest replaces the missing slice */
291       if ((lptr->last != NULL) && (lptr->last->to != NULL))
292         {
293           im_bef = lptr->last->to;
294         }
295       else 
296         {
297           im_bef = lptr->to;
298         }
299       if ((lptr->next != NULL) && (lptr->next->to != NULL))
300         {
301           im_aft = lptr->next->to;
302         }
303       else
304         {
305           im_aft = lptr->to;
306         }
307       im_new = im_median_3D(im_bef, im_cen, im_aft);
308       
309 
310       smoothed_link = link_alloc(im_new, IMRECT);
311       smoothed_images = list_addtoend(smoothed_images, smoothed_link);
312     }
313   /* can be deleted later, just for debugging */
314   smoothed_seq = seq_alloc();
315   smoothed_seq->start = smoothed_images;
316   smoothed_seq->current = smoothed_images;
317   smoothed_seq->end = list_get_end(smoothed_images);
318     
319   /* seq_set_current(smoothed_seq);*/
320 
321   /* calculate 3d spatial derivative (creates a new sequence of NMTT maps)*/
322   /* should also create the relevant masks */
323   /* 1.79688/(3*2) */
324 
325   mtt_scale = (ttm_seq->dim[0])/(2*ttm_seq->dim[2]);
326 
327   if (smoothed_images->next->to == NULL)
328     {
329       format("2nd image in seq is null");
330       return;
331     }
332 
333   for (lptr = smoothed_images->next; lptr->next != NULL; lptr = lptr->next)
334     {
335           /* bear in mind the sequence will now be two images shorter (2-24)*/
336       Imrect * im_cen = NULL, *im_bef = NULL, *im_aft = NULL;
337       Imrect *im_grad_h = NULL, *im_grad_v = NULL, *im_diff_z = NULL, *im_z = NULL;
338       Imrect *im_gradsq = NULL, *im_grad_z_sq = NULL, *im_gradsq_tot = NULL;
339       Imrect *im_nmtt = NULL;
340       List *nmtt_link = NULL;
341       
342       im_cen = lptr->to;
343       im_bef = lptr->last->to; 
344       im_aft = lptr->next->to;
345       
346       im_grad_h = imf_grad_h(im_cen); 
347       im_grad_v = imf_grad_v(im_cen);
348       im_gradsq = imf_sumsq(im_grad_h, im_grad_v);
349       
350       im_diff_z = im_diff(im_aft, im_bef);
351       im_z = im_times(mtt_scale, im_diff_z);
352       im_grad_z_sq = imf_prod(im_z, im_z);
353       
354       im_gradsq_tot = im_sum(im_gradsq, im_grad_z_sq);
355       im_nmtt = im_sqrt(im_gradsq_tot);
356       nmtt_link = link_alloc(im_nmtt, IMRECT);
357       nmtt_images = list_addtoend(nmtt_images, nmtt_link);
358       
359       im_free(im_grad_h);
360       im_free(im_grad_v);
361       im_free(im_gradsq);
362       im_free(im_diff_z);
363       im_free(im_z);
364       im_free(im_grad_z_sq);
365       im_free(im_gradsq_tot);
366 
367 
368     }
369       
370   nmtt_seq = seq_alloc();
371   nmtt_seq->start = nmtt_images;
372   nmtt_seq->current = nmtt_images;
373   nmtt_seq->end = list_get_end(nmtt_images);
374   nmtt_seq->dim[0] = ttm_seq->dim[0];
375   nmtt_seq->dim[1] = ttm_seq->dim[1];
376   nmtt_seq->dim[2] = ttm_seq->dim[2];
377   nmtt_seq->dim[3] = ttm_seq->dim[3];
378   nmtt_seq->offset = ttm_seq->offset+1;
379   nmtt_seq->stride = 1;
380   seq_set_current(nmtt_seq);
381     
382  
383 }
384 
385 /** @brief  Produces masks for applying to the NMTT maps to remove edge artefacts
386  *  @param  void
387  *  @return void
388  *
389  *  The NMTT maps have edge artefacts particularly at the boundary between the brain
390  *  and background. Masks to remove these artefacts can be created using a convolution of 
391  *  the binary masks of the smoothed TTM maps used to create the NMTT maps. 
392  *  More specifically, the TTM map at the slice of interest is binarised and barrel shifted 
393  *  by one pixel in all directions to produce 4 masks. These are convolved, along with masks 
394  *  slice above and below the slice of interest.
395  *  The NMTT masks are set to be the current sequence (visible in the Sequence Tool Tv).
396  *  NB, the volume of NMTT masks is 2 images shorter than the volume of TTM images used in 
397  *  their creation.
398  */
399 static void get_nmtt_masks_proc(void)
400 {
401   List *lptr = NULL, *nmtt_mask_list = NULL;
402 
403 
404   if (nmtt_mask_seq != NULL)
405     {
406       seq_rm(nmtt_mask_seq);
407       nmtt_mask_seq = NULL;
408     }
409 
410   if (smoothed_seq == NULL)
411     return;
412 
413   if (smoothed_seq->start->next == NULL)
414     return;
415 
416   for (lptr=smoothed_seq->start->next; lptr->next != NULL; lptr = lptr->next)
417     {
418       List *nmtt_mask_link = NULL;
419       Imrect *befa = NULL, *befb = NULL, *afta = NULL, *aftb = NULL;
420       Imrect *left = NULL, *right = NULL;
421       Imrect *up = NULL, *down = NULL;
422       Imrect *cena = NULL, *cenb = NULL;
423       Imrect *temp1=NULL, *temp2=NULL, *nmtt_mask = NULL; 
424 
425       befa = im_bthresh(0.0, lptr->last->to);
426       afta = im_bthresh(0.0, lptr->next->to);
427       cena = im_bthresh(0.0, lptr->to);
428 
429       befb = im_cast(befa, uchar_v);
430       aftb = im_cast(afta, uchar_v);
431       cenb = im_cast(cena, uchar_v);
432 
433       left = im_bshift(cenb, -1, 0);
434       right = im_bshift(cenb, 1, 0);
435       up = im_bshift(cenb, 0, -1);
436       down = im_bshift(cenb, 0, 1);
437 
438       temp1 = im_prod(befb, aftb);
439       temp2 = im_prod(temp1, left); im_free(temp1);
440       temp1 = im_prod(temp2, right); im_free(temp2);
441       temp2 = im_prod(temp1, up); im_free(temp1);
442       nmtt_mask = im_prod(temp2, down); im_free(temp2);
443       
444       im_free(befa); im_free(befb);
445       im_free(afta); im_free(aftb);
446       im_free(cena); im_free(cenb);
447       im_free(left); im_free(right);
448       im_free(up); im_free(down);
449 
450       nmtt_mask_link = link_alloc(nmtt_mask, IMRECT);
451       nmtt_mask_list = list_addtoend(nmtt_mask_list, nmtt_mask_link);
452 
453     }
454 
455   nmtt_mask_seq = seq_alloc();
456   nmtt_mask_seq->start = nmtt_mask_list;
457   nmtt_mask_seq->current = nmtt_mask_list;
458   nmtt_mask_seq->end = list_get_end(nmtt_mask_list);
459   nmtt_mask_seq->dim[0] = ttm_seq->dim[0];
460   nmtt_mask_seq->dim[1] = ttm_seq->dim[1];
461   nmtt_mask_seq->dim[2] = ttm_seq->dim[2];
462   nmtt_mask_seq->dim[3] = ttm_seq->dim[3];
463   nmtt_mask_seq->offset = ttm_seq->offset+1;
464   nmtt_mask_seq->stride = 1;
465   seq_set_current(nmtt_mask_seq);
466   
467 }
468 
469 /** @brief  Produces the log Net Cerebral Blood Flow images
470  *  @param  void
471  *  @return void
472  *
473  *  Calculates the log(NCBF) as described in Tina Memo 2002-001. 
474  *  The log(NCBF) sequence is set to the current sequence as visible in
475  *  the Sequence tool Tv.
476  */
477 static void get_flow_vol_proc(void)
478 {
479 
480   List *nmtt_lptr = NULL, *cbv_lptr = NULL, *mask_lptr = NULL;
481   double scale_factor = 0.0;
482   double voxel_length = 0.0;
483   List *flow_list = NULL;
484 
485   
486   if (flow_seq != NULL)
487     {
488       seq_rm(flow_seq);
489       flow_seq = NULL;
490     }
491   if ((nmtt_seq == NULL) || (cbv_seq == NULL) || (nmtt_mask_seq == NULL))
492     return;
493 
494   /* scale factor is 60secs*voxel length* (M/rho)^2/3 */
495   voxel_length = ttm_seq->dim[0]/10.0;
496   scale_factor = 60*voxel_length*20.4684;
497 
498    for (nmtt_lptr = nmtt_seq->start, mask_lptr = nmtt_mask_seq->start, cbv_lptr = cbv_seq->start->next; nmtt_lptr != NULL; nmtt_lptr = nmtt_lptr->next, mask_lptr = mask_lptr->next, cbv_lptr = cbv_lptr->next)
499     {
500       List *flow_link = NULL;
501       Imrect *flow_im = NULL, *nmtt_im = NULL, *cbv_im = NULL, *mask_im = NULL;
502       Imrect *temp1 = NULL, *temp2 = NULL, *cbv_sq = NULL, *cbv_sq_mask = NULL;
503       Imrect *log_flow_im = NULL, *thresh_log_flow_im = NULL;
504       double noise = 0.0, sigma_x = 0.0, sigma_y=0.0;
505 
506       nmtt_im = nmtt_lptr->to;
507       cbv_im = cbv_lptr->to;
508 
509       /* misnomer in nomenclature, the cbv_im is actually the sqrt cbv image */
510       mask_im = mask_lptr->to;
511       cbv_sq = im_prod(cbv_im, cbv_im);
512       
513       cbv_sq_mask = im_prod(cbv_sq, mask_im);
514       temp1 = im_prod(nmtt_im, mask_im);
515 
516       sigma_x = imf_diffx_noise(temp1, nmtt_im->region);
517       sigma_y = imf_diffy_noise(temp1, nmtt_im->region);
518       noise = (sigma_x + sigma_y)/2;
519 
520       temp2 = im_div(cbv_sq_mask, temp1, noise, noise);
521       flow_im = im_times(scale_factor, temp2);
522       log_flow_im = im_log(flow_im);
523       thresh_log_flow_im = im_thresh(-10.0, log_flow_im);
524 
525       im_free(temp1);
526       im_free(temp2);
527       im_free(cbv_sq);
528       im_free(cbv_sq_mask);
529       im_free(flow_im);
530       im_free(log_flow_im);
531       flow_link = link_alloc(thresh_log_flow_im, IMRECT);
532       flow_list = list_addtoend(flow_list, flow_link);
533      
534     }
535   flow_seq = seq_alloc();
536   flow_seq->start = flow_list;
537   flow_seq->current = flow_list;
538   flow_seq->end = list_get_end(flow_list);
539   flow_seq->dim[0] = ttm_seq->dim[0];
540   flow_seq->dim[1] = ttm_seq->dim[1];
541   flow_seq->dim[2] = ttm_seq->dim[2];
542   flow_seq->dim[3] = ttm_seq->dim[3];
543   flow_seq->offset = ttm_seq->offset+1;
544   flow_seq->stride = 1;
545   seq_set_current(flow_seq);
546  
547 }
548 
549 /** @brief  Creates (initialises) histograms for each Talairach region
550  *  @param  void
551  *  @return void
552  *
553  *  Initialises the histograms that will contain the regional (Taliarach) values 
554  *  of the flow maps. There are three sets of histograms corresponding to the whole brain 
555  *  and the left and right hemispheres. 
556  *  
557  */
558 static void init_tal_proc(void)
559 {
560  
561 
562   int i;
563 
564   if (shist_vec != NULL)
565     {
566       for (i=0; i<153; i++)
567         {
568           hfree(shist_vec[i]);
569         }
570      
571     }
572 
573   if (shist_vec_left != NULL)
574     {
575       for (i=0; i<153; i++)
576         {
577           hfree(shist_vec_left[i]); 
578         }
579          
580     }
581 
582   if (shist_vec_right != NULL)
583     {
584       for (i=0; i<153; i++)
585         {
586           hfree(shist_vec_right[i]);
587         }
588        
589     }
590   
591   shist_vec = (shistogram **)pvector_alloc(0, 153);
592   shist_vec_left = (shistogram **)pvector_alloc(0, 153);
593   shist_vec_right = (shistogram **)pvector_alloc(0, 153);
594 
595   for (i=0; i<153; i++)
596     {
597       shistogram *tal_hist = NULL;
598 
599       tal_hist = hbook1(i, pLexRecords[i], -10, 10, 50);
600       shist_vec[i] = tal_hist;
601       /*format("%d %s\n", i, pLexRecords[i]);*/
602     }
603 
604   for (i=0; i<153; i++)
605     {
606       shistogram *tal_hist_left = NULL;
607       shistogram *tal_hist_right = NULL;
608 
609       tal_hist_left = hbook1(i, pLexRecords[i], -10, 10, 50);
610       tal_hist_right = hbook1(i, pLexRecords[i], -10, 10, 50);
611       shist_vec_left[i] = tal_hist_left;
612       shist_vec_right[i] = tal_hist_right;
613       
614       }
615 
616   /* ideally also want something here to initialise flow to tal atlas in coreg tool if not already performed. */
617 
618   return;
619 
620 
621 
622 }
623 
624 /** @brief  Fills previously created regional histograms with log(NCBF) values
625  *  @param  void
626  *  @return void
627  *
628  *  Assumes that the volume of interest has already been registered to the Talairach 
629  *  template, or that the correct AIR file has been loaded. Requires the Coreg tool to 
630  *  be opened.
631  *  The function iterates over all the voxels in the log(NCBF) volume and projects the positions 
632  *  onto the Talairach atlas to determine which region each voxel falls in. If the voxel is
633  *  within the brain (not background) its value is added to the appropriate histograms. There 
634  *  are histograms covering a hierarchy of subdivisions, from lobar down through to the level of 
635  *  Brodmann area, and both left and right hemisphere and whole-brain histograms exist for all areas.
636  */
637 static void hist_flow_proc(void)
638 {
639   /* scan entire volumes
640      take flow value and input into correct histogram*/
641   float vox_val = 0.0;
642   Vec3 talpos, start_pos, tal3;
643   ListEntry *list_entry = NULL;
644   int i, j, k;
645   int seq_end = 0;
646   Imrect *im = NULL;
647   Imregion *imreg = NULL;
648   List *lptr = NULL;
649   int lx, ux, ly, uy;
650 
651   if (flow_seq == NULL)
652     {
653       format("There's no flow sequence \n");
654       return;
655     }
656 
657   format("Please ensure you have registered the flow volume to the Taliarach atlas or loaded in the appropriate air file.\n");
658 
659   seq_end = get_end_frame(flow_seq); 
660   lptr = flow_seq->start;
661   im = lptr->to;
662   imreg = im->region;
663   lx = imreg->lx; ly = imreg->ly; ux = imreg->ux; uy = imreg->uy;
664 
665   for (k = flow_seq->offset, lptr = flow_seq->start; k < seq_end+1, lptr != NULL; k++, lptr = lptr->next)
666     {
667       im = lptr->to;
668 
669       for (i = lx; i < ux; i++)
670         {
671  
672           for (j = ly; j < uy; j++)
673             {
674 
675               talpos = vec3(i, j, k); /*where ij and k are the voxel positions might need to add 0.5 to all */
676               vox_val = im_get_pixf(im, j, i); /* I HOPE :) */
677   /* I think because we've got the actual position that we don't need to add 0.5? */
678               tal3 = coreg_proj(talpos.el[0], talpos.el[1], talpos.el[2]+0.5);
679               /*!! this might prove problematic. need to do coreg_init or something, but I'm not sure which data set it needs to be done on */
680               tal3 = tal_convert(tal3);
681               list_entry = hist_fill_search(tal3);
682               
683               if ((list_entry != NULL) && (vox_val != 0.0000))
684                 {
685                   int gross = 0, lobe = 0, sublobar = 0, record = 0;
686                   int sublob1 = 0, sublob2 = 0;
687                   gross = list_entry->field1;
688                   hfill1(shist_vec[gross], vox_val, 1.0);
689                   
690                   lobe = list_entry->field2;
691                   hfill1(shist_vec[lobe], vox_val, 1.0);
692                   
693                   sublobar = list_entry->field3;
694                   hfill1(shist_vec[sublobar], vox_val, 1.0);
695                   
696                   sublob1 = list_entry->field4;
697                   hfill1(shist_vec[sublob1], vox_val, 1.0);
698                   
699                   sublob2 = list_entry->field5;
700                   hfill1(shist_vec[sublob2], vox_val, 1.0);
701 
702                   if (lobe >7 && lobe <53 && sublobar > 7 && sublobar <53)
703                     {
704                       if (gross == 2) /* ie, in left cerebrum */
705                         {
706                           hfill1(shist_vec_left[lobe], vox_val, 1.0);
707                           hfill1(shist_vec_left[sublobar], vox_val, 1.0);
708                           hfill1(shist_vec_left[sublob1], vox_val, 1.0);
709                           hfill1(shist_vec_left[sublob2], vox_val, 1.0);
710                         }
711                       if (gross == 3) /*ie, in right cerebrum */
712                         {
713                           hfill1(shist_vec_right[lobe], vox_val, 1.0);
714                           hfill1(shist_vec_right[sublobar], vox_val, 1.0);
715                           hfill1(shist_vec_right[sublob1], vox_val, 1.0);
716                           hfill1(shist_vec_right[sublob2], vox_val, 1.0);
717                         }
718                     }
719 
720                 }
721               else
722                 /*format("no desc here\n");*/
723                 {}
724             }
725         }
726     }
727   return;
728 }
729 
730 /** @brief Sets the value of the global hist_choice variable
731  *  @param int val The value determining histogram type 
732  *  @return void
733  *
734  *  Sets the value (0, 1, or 2) of the global variable determining histogram type 
735  *  (whole-brain, left or right hemisphere).
736  */
737 static void hist_choice_proc(int val)
738 {
739   hist_choice = val;
740 
741 }
742 
743 /** @brief  Plots the user-determined regional histogram of log(NCBF) values
744  *  @param  void
745  *  @return void
746  *
747  *  Takes the histopgram type (as specified by "Hist: Whole/Left/Right" on the user
748  *  interface) and histogram number (as specified by "Hist no :" on the user interface)
749  *  and plots the corresponding histogram on the Imcalc Graph Tv. Also provides the 
750  *  region name and histogram mean in the text box of the TinaTool.
751  */
752 static void hist_plot_proc(void)
753 {
754   Tv *graph_tv = NULL;
755   float hist_mean = 0.0;
756   shistogram **shist_choice = NULL;
757 
758   if ((graph_tv = imcalc_graph_tv_get()) == NULL)
759     {
760       format("Open and install Graph Tv first\n");
761       return;
762     }
763   
764   if ((shist_vec == NULL) || (shist_vec_left == NULL) || (shist_vec_right == NULL))
765     {
766       format("No histograms exist\n");
767       return;
768     }
769 
770   if ((hist_no >152) ||(hist_no < 1))
771     {
772       format("No such histogram %d\n", hist_no);
773       return;
774     }
775 
776   if (hist_no >=8 && hist_no <153)
777     {
778       switch (hist_choice)
779         {
780         case 0: 
781           {
782             shist_choice = shist_vec;
783             break;
784           }
785         case 1:
786           {
787             shist_choice = shist_vec_left;
788             break;
789           }
790         case 2:
791           {
792             shist_choice = shist_vec_right;
793             break;
794           }
795         default:
796           shist_choice = shist_vec;
797         }
798     }
799 
800   else
801     {
802       shist_choice = shist_vec;
803     }
804   graph_hfit(graph_tv, shist_choice[hist_no]);
805   hist_mean = shist_choice[hist_no]->mean;
806   format("Histogram of: %s ", pLexRecords[hist_no]);
807   format("mean: %f\n", hist_mean);
808   
809   return;
810 }
811 
812 
813 /** @brief Code for the flow tool interface
814  *  @param int x  tool x-dimension
815  *  @param int y  tool y-dimension
816  *  @return void
817  *
818  *
819  */
820 void flow_tool(int x, int y)
821 {
822   static void *tool = NULL;
823 
824   if (tool)
825     {
826       tw_show_tool(tool);
827       return;
828     }
829   
830   /*tv_proc();*/
831     
832     tool = (void *)tw_tool("Flow Calculation Tool", x, y);
833     
834     tw_label("Initialise Image Volumes");
835     tw_newrow();
836     tw_button("CBV volume", get_cbv_vol_proc, NULL);
837     tw_button("TTM volume", get_ttm_vol_proc, NULL);
838     tw_newrow();
839     tw_label("Perform Flow Calculations");
840     tw_newrow();
841     tw_button("Process CBV", process_cbv_proc, NULL);
842     tw_label(" -> ");
843     tw_button("process TTM", process_ttm_proc, NULL);
844     tw_label(" -> ");
845     tw_button("NMTT Masks", get_nmtt_masks_proc, NULL);
846     tw_label(" -> ");
847     tw_button("Get Flow Vol", get_flow_vol_proc, NULL);
848     tw_newrow();
849     tw_label("Initialise Regional Histograms");
850     tw_newrow();
851     tw_button("Init Hist.", init_tal_proc, NULL);
852     tw_button("Calc Flow Hists", hist_flow_proc, NULL);
853     tw_newrow();
854     tw_choice("Hist: ", hist_choice_proc, 0, "Whole", "Left", "Right", NULL);
855     ptalno = (void *)tw_iglobal("   Hist no: ", &hist_no, 10);
856     tw_button("Plot Hist", hist_plot_proc, NULL);
857 
858     tw_end_tool();
859 }
860 
861 
862 
863 
864 
865 
866 

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