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

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

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