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

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedCort_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    :  $Source: /home/tina/cvs/tina-tools/tinatool/tlmedical/tlmedCort_tool.c,v $
 23  * Date    :  $Date: 2007/02/15 01:55:50 $
 24  * Version :  $Revision: 1.3 $
 25  * CVS Id  :  $Id: tlmedCort_tool.c,v 1.3 2007/02/15 01:55:50 paul Exp $
 26  *
 27  * Author  :  mjs
 28  *
 29  * Notes : Tool and functionality to perform cerebral grey matter cortical 
 30  *         thickness estimation.
 31  *
 32  *********
 33 */
 34 
 35 #include "tlmedCort_tool.h"
 36 
 37 #if HAVE_CONFIG_H
 38   #include <config.h>
 39 #endif
 40 
 41 #include <stdio.h>
 42 #include <sys/param.h>
 43 #include <string.h>
 44 #include <math.h>
 45 
 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 <tinatool/draw/drawDef.h>
 51 #include <tinatool/draw/drawPro.h>
 52 #include <tina/image/imgDef.h>
 53 #include <tina/image/imgPro.h>
 54 #include <tina/geometry/geomDef.h>
 55 #include <tina/geometry/geomPro.h>
 56 #include <tina/medical/medDef.h>
 57 #include <tina/medical/medPro.h>
 58 #include <tinatool/draw/drawDef.h>
 59 #include <tinatool/draw/drawPro.h>
 60 #include <tinatool/wdgts/wdgtsDef.h>
 61 #include <tinatool/wdgts/wdgtsPro.h>
 62 #include <tinatool/tlvision/tlvisEdge_epi_mouse.h>
 63 #include <tinatool/tlbase/tlbaseDef.h>
 64 #include <tinatool/tlbase/tlbasePro.h>
 65 #include <tinatool/tlmedical/tlmedDef.h>
 66 #include <tinatool/tlmedical/tlmed_TalDef.h>
 67 #include <tinatool/tlmedical/tlmed_TalPro.h>
 68 #include <tinatool/tlmedical/tlmed_CoregPro.h>
 69 #include <tinatool/tlmedical/tlmedCort_funcs.h>
 70 
 71 static void *psi, *ppr, *plt, *put, *plen, *pmp, *pgr, *pcsf, *pext, *pfrac, *pmind, *pgm, *pedl;
 72 static void *ptalno, *pssp;
 73 static double sigma = 1.0; 
 74 static double precision = 0.01;
 75 static double low_edge_thres = 0.0;
 76 static double up_edge_thres = 0.9;
 77 static int len_thres = 5;
 78 static double grey_mid = -300;
 79 static double grey_mean = -500;
 80 static int extent = 20;
 81 static double edge_limit = 4.0;
 82 static double fraction = 0.5;
 83 static double skull_seg = 3.05;
 84 static int hist_no = 1;
 85 
 86 static List *gm_er_list = NULL;
 87 static List *bound_im_list = NULL;
 88 
 89 static Sequence *gm_seq = NULL; /*grey matter prob sequence */
 90 static Sequence *seq_gauss = NULL;
 91 
 92 static int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
 93 static void ***imptrs=NULL, ***imptrs_gm = NULL, ***imptrs_gm_prob = NULL; 
 94 static int colour_choice = 258; /*blue */
 95 static shistogram **shist_vec = NULL, **shist_vec_left = NULL, **shist_vec_right = NULL;
 96 static int hist_choice = 0;
 97 static List *surf_list= NULL;
 98 
 99 static int csf_edge = 0, csf_dip = 0, wm_edge = 0, wm_dip = 0, no_edge_dip = 0;
100 
101 static int count_line = 0, count_line_rm = 0;
102 
103 /* Returns the skull segmentation threshold parameter*/
104 double get_skull_seg_param(void)
105 {
106   return(skull_seg);
107 }
108 
109 /* Wrapper function for taking the grey matter probability images and calling 
110    "get_boundary_im" for finding the boundary of the skull. */
111 static void get_boundary_proc(void)
112 {
113   
114   List *gm_prob_list = NULL;
115   List *lptr = NULL;
116   Tv *imc_tv = NULL;
117 
118   if (gm_seq == NULL)
119     {
120       format("No GM prob sequence loaded\n");
121       return;
122     }
123 
124   if ((imc_tv = imcalc_tv_get()) == NULL)
125     {
126       format("Open imcalc tv first!\n");
127       return;
128     }
129   
130 
131   if (bound_im_list != NULL)
132     {
133       list_rm(bound_im_list, im_free);
134       bound_im_list = NULL;
135     }
136 
137   gm_prob_list = gm_seq->start;
138   
139   for (lptr = gm_prob_list; lptr != NULL; lptr = lptr->next)
140     {
141       /* go through GM prob list and obtain the boundary image */
142       Imrect *bound_im = NULL, *gm_im;
143       List *bound_el = NULL;
144      
145       gm_im = lptr->to; 
146       bound_im = get_boundary_im(gm_im);
147 
148       bound_el = link_alloc(bound_im, IMRECT);
149       bound_im_list = list_addtoend(bound_im_list, bound_el);
150                   
151     }
152 
153    return;
154 
155 }
156 
157 /* Frees all the stored histograms from memory*/
158 static void rm_hist_proc(void)
159 {
160   int i;
161 
162   if (shist_vec != NULL)
163     {
164       for (i=0; i<153; i++)
165         {
166           hfree(shist_vec[i]);
167         }
168       pvector_free(shist_vec, 0);
169       
170     }
171   
172 }
173 
174 /* Allocates and initialises all the required histograms for the bootstrap 
175    offset analysis ("check_canny_proc")*/
176 static void tal_hist_init_temp(void)
177 {
178 
179   int i;
180   List *lptr = NULL;
181   double sx, sy, sz;
182 
183   if (shist_vec != NULL)
184     {
185       for (i=0; i<153; i++)
186         {
187           hfree(shist_vec[i]);
188         }
189       pvector_free(shist_vec, 0);
190       
191     }
192 
193   if (shist_vec_left != NULL)
194     {
195       for (i=8; i<153; i++)
196         {
197           hfree(shist_vec_left[i]);
198         }
199       pvector_free(shist_vec_left, 8);
200       
201     }
202 
203   if (shist_vec_right != NULL)
204     {
205       for (i=8; i<153; i++)
206         {
207           hfree(shist_vec_right[i]);
208         }
209       pvector_free(shist_vec_right, 8);
210       
211     }
212   shist_vec = (shistogram **)pvector_alloc(0, 153);
213   shist_vec_left = (shistogram **)pvector_alloc(8, 153);
214   shist_vec_right = (shistogram **)pvector_alloc(8, 153);
215 
216   for (i=0; i<153; i++)
217     {
218       shistogram *tal_hist = NULL;
219 
220       tal_hist = hbook1(i, pLexRecords[i], -200, 200, 30);
221       shist_vec[i] = tal_hist;
222     }
223 
224   for (i=8; i<153; i++)
225     {
226       shistogram *tal_hist_left = NULL;
227       shistogram *tal_hist_right = NULL;
228 
229       tal_hist_left = hbook1(i, pLexRecords[i], -200, 200, 30);
230       tal_hist_right = hbook1(i, pLexRecords[i], -200, 200, 30);
231       shist_vec_left[i] = tal_hist_left;
232       shist_vec_right[i] = tal_hist_right;
233     }
234 
235 
236   return;
237 }
238 
239 /* Allocates and initialises all the histograms for the regional cortical 
240  thickness estimation*/
241 static void tal_hist_init(void)
242 {
243 
244   int i;
245   List *lptr = NULL;
246   double sx, sy, sz;
247 
248   if (shist_vec != NULL)
249     {
250       for (i=0; i<153; i++)
251         {
252           hfree(shist_vec[i]);
253         }
254       pvector_free(shist_vec, 0);
255       
256     }
257 
258   if (shist_vec_left != NULL)
259     {
260       for (i=8; i<153; i++)
261         {
262           hfree(shist_vec_left[i]);
263         }
264       pvector_free(shist_vec_left, 8);
265       
266     }
267 
268   if (shist_vec_right != NULL)
269     {
270       for (i=8; i<153; i++)
271         {
272           hfree(shist_vec_right[i]);
273         }
274       pvector_free(shist_vec_right, 8);
275       
276     }
277   shist_vec = (shistogram **)pvector_alloc(0, 153);
278   shist_vec_left = (shistogram **)pvector_alloc(8, 153);
279   shist_vec_right = (shistogram **)pvector_alloc(8, 153);
280 
281   for (i=0; i<153; i++)
282     {
283       shistogram *tal_hist = NULL;
284 
285       tal_hist = hbook1(i, pLexRecords[i], 0, extent/2.0, 50);
286       shist_vec[i] = tal_hist;
287     }
288 
289   for (i=8; i<153; i++)
290     {
291       shistogram *tal_hist_left = NULL;
292       shistogram *tal_hist_right = NULL;
293 
294       tal_hist_left = hbook1(i, pLexRecords[i], 0, extent/2.0, 50);
295       tal_hist_right = hbook1(i, pLexRecords[i], 0, extent/2.0, 50);
296       shist_vec_left[i] = tal_hist_left;
297       shist_vec_right[i] = tal_hist_right;
298     }
299 
300 
301   return;
302 }
303 
304 /* Takes the line structure describing the GM/WM boundary start point and the 
305    associated thickness (length), determines the Talairach region the start 
306    point is in and inserts the thickness into the respective histograms
307 */
308 static void tal_hist_fill(Line3 *line, Sequence *seq)
309 {
310 
311   float len = 0.0;
312   Vec3 talpos, start_pos, tal3;
313   ListEntry *list_entry = NULL;
314   double sx = 0.0, sy = 0.0, sz = 0.0;
315     
316   sx = seq->dim[0];
317   sy = seq->dim[1];
318   sz = seq->dim[2];
319   
320   len = line->length;
321   start_pos = line->p1;
322   
323   talpos = vec3(start_pos.el[0]/sx, start_pos.el[1]/sy, start_pos.el[2]/sz);
324   /* I think because we've got the actual position that we don't need to add 0.5? */
325   tal3 = coreg_proj(talpos.el[0], talpos.el[1], talpos.el[2]+0.5);
326   tal3 = tal_convert(tal3);
327   list_entry = hist_fill_search(tal3);
328   
329   if (list_entry != NULL)
330     {
331       int gross = 0, lobe = 0, sublobar = 0, gm_wm = 0, area = 0;
332       gross = list_entry->field1;
333       hfill1(shist_vec[gross], len, 1.0);
334       
335       lobe = list_entry->field2;
336       hfill1(shist_vec[lobe], len, 1.0);
337       
338       sublobar = list_entry->field3;
339       hfill1(shist_vec[sublobar], len, 1.0);
340       
341       
342       gm_wm = list_entry->field4;
343       hfill1(shist_vec[gm_wm], len, 1.0);
344       
345       area = list_entry->field5;
346       hfill1(shist_vec[area], len, 1.0);
347 
348       if (gross == 2) /* ie, in left cerebrum */
349         {
350           hfill1(shist_vec_left[lobe], len, 1.0);
351           hfill1(shist_vec_left[sublobar], len, 1.0);
352           hfill1(shist_vec_left[gm_wm], len, 1.0);
353           hfill1(shist_vec_left[area], len, 1.0);
354         }
355       if (gross == 3) /*ie, in right cerebrum */
356         {
357           hfill1(shist_vec_right[lobe], len, 1.0);
358           hfill1(shist_vec_right[sublobar], len, 1.0);
359           hfill1(shist_vec_right[gm_wm], len, 1.0);
360           hfill1(shist_vec_right[area], len, 1.0);
361         }
362         
363            
364     }
365   else
366     {}
367   
368 return;
369 
370 }
371 
372 
373 /* Takes the voxel position (i,j) of the (more precise) current GM/WM 
374    boundary position and a pointer to the list of the Gaussian smoothed 
375    T1 images and determines the 3D surface normal to the GM boundary at 
376    that voxel based on the adjacent 8 voxels. 
377  */
378 static Vec3 get_3D_orient(List *seq_list_ptr, int i/*column*/, int j/*row*/)
379 {
380 
381   Vec3 orient, or_res;
382   Sequence *seq = NULL;
383   Imrect *im_cen = NULL, *im_bel = NULL, *im_abv = NULL;
384   double x1=0.0, x2=0.0, y1=0.0, y2=0.0, z1=0.0, z2=0.0, cen = 0.0;
385   double sx=0.0, sy=0.0, sz=0.0;
386 
387   if ((seq = seq_get_current()) == NULL)
388     return;                     /* fixme: needs value */
389   
390   sx = seq->dim[0];
391   sy = seq->dim[1];
392   sz = seq->dim[2];
393   
394   im_cen = seq_list_ptr->to;
395   cen = im_get_pixf(im_cen, j, i);   /*image, row, column*/
396   x1 = im_get_pixf(im_cen, j, i-1); 
397   x2 = im_get_pixf(im_cen, j, i+1);
398   y1 = im_get_pixf(im_cen, j-1, i);
399   y2 = im_get_pixf(im_cen, j+1, i); 
400 
401   if (seq_list_ptr->last !=NULL) im_bel = seq_list_ptr->last->to; 
402   else im_bel = NULL;
403 
404   if (seq_list_ptr->next !=NULL) im_abv = seq_list_ptr->next->to;
405   else im_abv = NULL;
406 
407   if ((im_bel != NULL) && (im_abv != NULL))
408     {
409 
410       z1 = im_get_pixf(im_bel, j, i);
411       z2 = im_get_pixf(im_abv, j, i);
412       orient = vec3((x2-x1)/sx, (y2-y1)/sy, (z2-z1)/sz);
413       /*format("using above and below\n");*/
414     }
415 
416   else
417     {
418       if ((im_bel == NULL) && (im_abv != NULL))
419         {
420           z2 = im_get_pixf(im_abv, j, i);
421           orient = vec3((x2-x1)/sx, (y2-y1)/sy, ((z2-cen)*2)/sz);
422         }
423       else
424         {
425           if ((im_bel != NULL) && (im_abv == NULL))
426             {
427               z1 = im_get_pixf(im_bel, j, i);
428               orient = vec3((x2-x1)/sx, (y2-y1)/sy, ((cen-z1)*2)/sz);
429             }
430           else
431             {
432               /* at this point, z1 and z2 should still be zero */
433               orient = vec3((x2-x1)/sx, (y2-y1)/sy, 0.0);
434             }
435         }
436     }
437   /* This has already implicitly taken account of the scale factors */
438   /* so we have a unit vector pointing in the correct direction. */
439   or_res = vec3_unit(vec3_minus(orient));
440   return(or_res);
441 }
442 
443 /* For a given GM/WM boundary starting point (start_pos) and associated 
444    3D surface normal through the GM (unit_vec), get_frac determines the 
445    position of the opposing GM boundary. 
446  */
447 static Vec3 get_frac(Vec3 start_pos, Vec3 unit_vec, double frac)
448 {
449   /* start_pos is real world */
450   Vec3 end_pos_edge, end_pos_dip;
451   int i=0, j=0, n=0;
452   Sequence *seq = NULL;
453   double scalex=0.0, scaley=0.0, scalez=0.0;
454   Vec3 new_pos;
455   float new_val = 0.0, old_val = 0.0, gm_val=0.0;
456   double loc_frac = 0.0;
457   int new_slice;
458   int seq_start = 0, seq_end = 0;
459   List *lptr = NULL, *lptr_gm = NULL;
460   Imrect *grey_im = NULL, *prob_im = NULL;
461   int temp_pos = 0;
462   double final_pos = 0.0;
463   int edge_type = 0, dip_type = 0;
464   /*nb, for edges and dips, 0=no edge/dip, 1 = csf edge/dip, 2 = wm edge/dip */
465 
466   float *gm_vec = NULL, *gm_prob_vec = NULL;
467   float loc_avg = 0.0;
468   float loc_dist = 0.0;
469   float max_dist = 0.0;
470   int max_dist_pos = 0;
471   float max_dist_white = 0.0, max_dist_csf = 0.0;
472   int max_dist_pos_white = 0, max_dist_pos_csf = 0;
473   Vec3 unscaled_pos, unscaled_old_pos;
474   int dip_found = 0;
475 
476 
477   if ((seq = seq_get_current()) == NULL)
478     return;                     /* fixme: needs value */
479   
480   grey_im = get_seq_start_el(seq)->to;
481   prob_im = get_seq_start_el(gm_seq)->to;
482 
483   scalex = seq->dim[0];
484   scaley = seq->dim[1];
485   scalez = seq->dim[2];
486   
487   seq_start = seq->offset;
488   seq_end = get_seq_end();
489 
490    
491   unscaled_old_pos = vec3(start_pos.el[0]/scalex, start_pos.el[1]/scaley, start_pos.el[2]/scalez);
492   /* ie, in the image voxel co-ords from the actual real world brain co-ords*/
493   seq_voxel_vtype(prob_im->vtype);
494   /* sets the GM probability value at the start*/ 
495  old_val = seq_interp(imptrs_gm_prob, unscaled_old_pos);
496     
497   gm_vec = fvector_alloc(0, extent+1);
498   gm_prob_vec = fvector_alloc(0, extent+1);
499   
500   temp_pos = extent; /* only want end point to change if stopping criteria met */
501 
502   /* from the start point to the search extent, in increments of 1mm 
503      (unit vectors), arrays of the T1 image intensity values (gm_val) and 
504      the GM probability values (gm_prob_vec) are filled. 
505    */
506   for (i=0;i<extent+1 ; i++)
507     {
508       new_pos = vec3_sum(start_pos, vec3_times((double)i, unit_vec));
509       unscaled_pos = vec3(new_pos.el[0]/scalex, new_pos.el[1]/scaley, new_pos.el[2]/scalez);
510   
511       seq_voxel_vtype(grey_im->vtype);     
512       gm_val = seq_interp(imptrs_gm, unscaled_pos);
513       seq_voxel_vtype(prob_im->vtype);
514       new_val = seq_interp(imptrs_gm_prob, unscaled_pos);
515        
516       
517       gm_vec[i] = gm_val;
518       gm_prob_vec[i] = new_val;
519            
520     }
521   /* this for loop searches the GM probability array for the first occurrence
522      of a true edge, where the GM probability drops below the value defined by 
523      frac (default is 0.5)
524    */
525   for (i=0;i<extent+1; i++)
526     {
527       /* If both the image intensities and probabilities drop to 0.0, then
528          the edge of the brain has been reached and there is no further data.*/
529       if ((gm_vec[i] == 0.0) && (gm_prob_vec[i] == 0.0))
530         {
531           break;
532         }
533       
534       new_val = gm_prob_vec[i];
535       /* if the current GM prob value is less than the user defined fraction, 
536          and the previous value is greater then an edge has been found.*/
537       if ((new_val<frac) && (old_val>frac) && (new_val != 0.0))
538         {
539           float pos1, pos2;
540           temp_pos = i+1;
541           
542           pos1 = gm_prob_vec[temp_pos-2];
543           pos2 = gm_prob_vec[temp_pos-1];
544           loc_frac = (double)(pos1-frac)/(pos1-pos2);
545           
546           /* GM/CSF boundary*/
547           if (gm_vec[temp_pos-1] < grey_mean)
548             {
549               float temp =0.0;
550 
551               end_pos_edge = vec3_sum(start_pos, vec3_times((double)(temp_pos-2)+loc_frac, unit_vec));
552               edge_type = 1; /* found CSF boundary*/
553               temp = (float)i;
554               
555               if (temp>=extent/2.0)
556                 edge_type = 0; /* CSF boundary not within extent/2*/
557             }
558           /* GM/WM boundary*/
559           else
560             {
561               end_pos_edge = vec3_sum(start_pos, vec3_times((double)(temp_pos-2+loc_frac)/2.0, unit_vec));
562               edge_type = 2;
563             }
564 
565           break; /* the first occurrence of an edge has been found, so stop */
566           /* but this doesn't account for the case where there might be a WM 
567              edge further on if a CSF edge hasn't been found */
568         }
569       old_val = new_val;
570     }
571   
572   /* Now search the GM prob. vector for possible (partial volume) dips in value
573    */
574    for (j=1; j<extent; j++) 
575     {
576       /* Estimate the difference between the average of the previous and next
577        values with the current value*/
578       loc_avg = (gm_prob_vec[j-1] + gm_prob_vec[j+1])/2.0;
579       loc_dist = loc_avg - gm_prob_vec[j];
580       
581       /* if the current GM prob value is lower than those either side*/
582       if ((gm_prob_vec[j] < gm_prob_vec[j-1]) && (gm_prob_vec[j] < gm_prob_vec[j+1]))
583       {
584         float temp = (float)j;
585 
586         /* if it is a CSF dip, it's the largest CSF dip found and it's within
587          extent/2 then set this to the current CSF dip*/
588         if ((gm_vec[j] < grey_mean) && (loc_dist > max_dist_csf) && (temp < extent/2.0))
589           {
590             max_dist_csf = loc_dist;
591             max_dist_pos_csf = j;
592             dip_found = 1;
593           }
594         /* if it is a WM dip, it's the largest WM dip found then set this 
595            to the current CSF dip*/
596         if ((gm_vec[j] > grey_mean) && (loc_dist > max_dist_white))
597           {
598             max_dist_white = loc_dist;
599             max_dist_pos_white = j;
600             dip_found = 1;
601           }
602       }
603     }
604 
605    if (dip_found == 1)
606      {
607        /* If a dip has been found, then determine whether the CSF or WM dip 
608           is the largest*/
609        if (max_dist_csf > max_dist_white)
610          max_dist_pos = max_dist_pos_csf;
611        if (max_dist_white > max_dist_csf)
612          max_dist_pos = max_dist_pos_white;
613        
614        /* find the precise position of the dip using a linear approximation 
615           to quadratic that it is assumed the dip forms, offset because you 
616           know that the centre of the dip occurs some fraction (loc_dist/2) 
617           after the edge of the dip (see Tina Memo 2004-007).*/
618        j = max_dist_pos;
619        if (gm_prob_vec[j+1] > gm_prob_vec[j-1])
620          {
621            loc_frac = (gm_prob_vec[j+1] - gm_prob_vec[j-1])/(gm_prob_vec[j+1] - gm_prob_vec[j]);
622            final_pos = (double)j-1+loc_frac - (loc_dist/2);
623          }
624        else
625          {
626            loc_frac = (gm_prob_vec[j-1] - gm_prob_vec[j+1])/(gm_prob_vec[j-1] - gm_prob_vec[j]);
627            final_pos = (double)j+1 - loc_frac - (loc_dist/2); 
628          }
629        /* set dip type as CSF or WM, although this is duplication of the
630           above :-( */
631        if (gm_vec[j] < grey_mean)
632          {
633            end_pos_dip = vec3_sum(start_pos, vec3_times(final_pos, unit_vec));
634            dip_type = 1;
635          }
636        else 
637          {
638            end_pos_dip = vec3_sum(start_pos, vec3_times(final_pos/2, unit_vec));
639            dip_type = 2;
640          }
641      }
642 
643   fvector_free(gm_vec, 0);
644   fvector_free(gm_prob_vec, 0);
645    
646   /* set the colours for the overlay display based on edge/dip type*/
647   if (edge_type == 1)
648     {
649       colour_choice = orange;
650       csf_edge++;
651       return(end_pos_edge);
652     }
653   if ((edge_type == 2) && ((dip_type == 2) || (dip_type == 0)))
654     {
655       colour_choice = green;
656       wm_edge++;
657       return(end_pos_edge);
658     }
659   
660   if ((edge_type == 2) && (dip_type == 1))
661     {
662       colour_choice = blue;
663       csf_dip++;
664       return(end_pos_dip);
665     }
666   if ((edge_type == 0) && (dip_type != 0))
667     {
668       if (dip_type == 1)
669         {
670           colour_choice = blue;
671           csf_dip++;
672         }
673       else 
674         {
675           colour_choice = red;
676           wm_dip++;
677         }
678       return(end_pos_dip);
679     }
680   
681 
682   end_pos_edge = vec3(start_pos.el[0], start_pos.el[1], start_pos.el[2]);
683   colour_choice = black;
684   no_edge_dip++;
685   return(end_pos_edge);
686 
687 }
688 
689 /* For the current slice, draw the overlay of the start points and lines 
690    through the grey matter, as a 2D representation of the 3D case. */
691 static void overlay_proc(void)
692 {
693 
694   List *lptr_er = NULL, *string_list = NULL, *lptr2 = NULL;
695   Tv * tv = tv_get_next();
696   Line3 * line;
697   Vec2 start, end;
698   Sequence *seq = NULL;
699   double scalex =0.0, scaley =0.0, scalez = 0.0;
700   Graphic *gr = NULL;
701   int slice_no = 0, temp_slice=0;
702   Imrect *er_im = NULL;
703 
704   if ((seq = seq_get_current()) == NULL)
705     return;
706   
707   if (gm_er_list == NULL)
708     {
709       return;
710     }
711   
712   scalex = seq->dim[0];
713   scaley = seq->dim[1];
714   scalez = seq->dim[2];
715 
716   slice_no = get_slice_no();
717 
718  
719   for (lptr_er = gm_er_list,temp_slice=0 ; lptr_er != NULL, temp_slice<slice_no; lptr_er = lptr_er->next, temp_slice++);
720     
721   er_im = lptr_er->to; 
722   string_list = (List*) prop_get(er_im->props, STRING);
723       
724   for(lptr2 = string_list; lptr2 != NULL; lptr2 = lptr2->next)
725     {
726       List *tptr = NULL;
727       Tstring *edge_string = lptr2->to;
728       Edgel *edge_temp = NULL;
729           
730       
731       for (tptr = edge_string->start; /*tptr != NULL*/; tptr=tptr->next)
732         {
733           edge_temp = tptr->to;
734          
735           if ((line = prop_get(edge_temp->props, LINE3))!=NULL)
736             {
737          
738               int colour;
739               
740               if ((gr = prop_get(line->props,GRAPHIC))==NULL)
741                 colour = black;
742               else 
743                 colour = gr->colour;
744               
745               start = vec2(line->p1.el[0]/scalex, line->p1.el[1]/scaley);
746               end = vec2(line->p2.el[0]/scalex, line->p2.el[1]/scaley);
747               tv->cx = 128;
748               tv->cy = 128;
749               
750               tv_color_set(tv, colour);
751               tv_bigdot2(tv, start, 1);
752               tv_line2(tv, start, end);
753               tv_bigdot2(tv, end, 1);
754                   
755             }
756           else
757             format("null line detected\n");
758           
759           if (tptr == edge_string->end)
760             break;
761         }
762     }
763      
764 }
765 
766 /* calls the overlay_proc function*/
767 static void overlay_line_proc(void)
768 {
769   overlay_proc();
770 
771 }
772 
773 /* Takes a position (edgel) on the GM/WM boundary, calls get_frac to 
774    determine the opposing GM edge for that starting position, and puts 
775    the results in a line structure (which contains start, end and cortical 
776    thickness length) which is then returned.*/
777 static Line3 *find_min_dist(Edgel *edge, int i, int j, int slice_no)
778 {
779 
780   Sequence *seq = NULL;
781   List *seq_list_ptr = NULL;
782   double scalex, scaley, scalez;
783   Graphic *gr = NULL;
784   int m, n;
785 
786   if ((seq = seq_get_current()) == NULL)
787     return(NULL);
788   
789   if ((seq_list_ptr = get_seq_current_el(seq_gauss)) == NULL) /* this current is going to stay constant isnt it */
790     return(NULL);
791       
792 
793   if (edge == NULL)
794     return(NULL);
795   
796   
797   scalex = seq->dim[0];
798   scaley = seq->dim[1];
799   scalez = seq->dim[2];
800   
801   if (bound_im_list != NULL)
802     {
803       int bound_count=0, a=0;
804       List *lptr = NULL;
805       double res = 0.0;
806       Imrect *bound_im = NULL;
807     
808       
809       lptr = bound_im_list;
810       bound_count=seq->offset;
811       while (lptr != NULL && bound_count<slice_no)
812         {
813           lptr = lptr->next; 
814           bound_count += seq->stride;
815         }
816            
817       bound_im = lptr->to;
818       if (((double)im_get_pixf(bound_im, j, i)) == 1.0)
819         return(NULL); /* ie, if it's exterior to the brain.*/
820     }
821 
822   
823   if (((edge->type & EDGE_GET_CONN_MASK) == EDGE_CONN)|| ((edge->type & EDGE_GET_CONN_MASK) == EDGE_TERMIN))
824     {
825            
826       Line3 *conn_line = NULL;
827       List *line_el = NULL;
828       Vec3  e3_pos, unit_vec;
829       Vec3 e1_pos;
830       
831       e1_pos = vec3(edge->pos.el[0]*scalex, edge->pos.el[1]*scaley, (double)slice_no*scalez);
832       
833       unit_vec = get_3D_orient(seq_list_ptr, i, j); /*ie, ptr, column(x), row(y)*/
834       e3_pos = get_frac(e1_pos, unit_vec, fraction);
835       /*e1_pos is in the realworld brain co-ords */
836       conn_line = line3_make(e1_pos, e3_pos, LINE3); 
837       if (vec3_dist(e1_pos, e3_pos)==0.0) 
838         {
839           conn_line->v.el[0] = unit_vec.el[0];
840           conn_line->v.el[1] = unit_vec.el[1];
841           conn_line->v.el[2] = unit_vec.el[2];
842         }
843     
844       if ((gr = prop_get(conn_line->props,GRAPHIC))==NULL) 
845         {
846           gr=(Graphic *)ralloc(sizeof(Graphic));
847           conn_line->props = proplist_add(conn_line->props, (void *)gr, GRAPHIC, rfree);
848         }
849       gr->colour = colour_choice;
850       
851       return(conn_line);
852            
853     }
854  
855 }
856 
857 /* For every edge voxel in the image (ie every voxel on the GM/WM boundary), 
858    calls find_min_dist to get the line structure associated with the edgel  
859    structure start point. If the length of the line is not zero, tal hist fill
860    is called to insert the length into the appropriate regional histogram */
861 static void find_min_dist_proc(List *gm_list, int slice_no) /* this is not a list, it is a pointer in list*/
862 {
863   
864   int i, j;
865   int ux, lx, uy, ly;
866   Imrect *im = NULL;
867   Tv *tv  = NULL;
868   List * lptr = NULL;
869   Imrect *gm_er = NULL;
870   Imregion *imreg = NULL;
871   Sequence *seq = NULL;
872   double scalex = 0.0, scaley = 0.0, scalez = 0.0;
873   List *string_list = NULL;
874 
875   if ((seq = seq_get_current()) == NULL)
876     return;
877   /*gm_er = im_copy(gm_list->to);*/
878   gm_er = gm_list->to;
879 
880   tv = seq_tv_get();
881 
882 
883   if (tv != NULL && tv->tv_screen != NULL)    /* get roi from tv */
884     imreg = tv_get_im_roi(tv);
885 
886   ux = imreg->ux;
887   lx = imreg->lx;
888   uy = imreg->uy;
889   ly = imreg->ly;
890   
891  
892   string_list = (List*) prop_get(gm_er->props, STRING);
893       
894   for(lptr = string_list; lptr != NULL; lptr = lptr->next)
895     {
896       List *tptr = NULL;
897       Tstring *edge_string = lptr->to;
898 
899       for (tptr = edge_string->start; /*tptr != NULL*/; tptr=tptr->next)
900         {
901           Edgel *edge = NULL;
902           Line3 *line = NULL;
903           List *line_el = NULL;
904 
905           edge = tptr->to;
906           i = tina_int(edge->pos.el[0]);
907           j = tina_int(edge->pos.el[1]); 
908 
909           line = find_min_dist(edge, i, j, slice_no);
910 
911         
912           if (line != NULL)
913             {
914               
915               line_el = link_alloc(line, LINE3);
916               edge->props = proplist_add(edge->props, (void *) line, LINE3, line3_free);
917               if (line->length != 0.0)
918                 tal_hist_fill(line, seq);
919             }
920           if (tptr == edge_string->end)
921             break;
922         }
923     }
924   
925   rfree(imreg);
926   
927   return;  
928 
929 }
930 
931 /* takes the image in the mono tool and applies the iso-canny algorithm 
932    according to the parameters in the cortical thickness dialogue box.*/
933 static void test_canny_proc(void)
934 {
935 
936  Tv     *tv = NULL;   
937  Imrect *im = NULL;
938  Imrect *gm_er = NULL;
939  Imrect *csf_er = NULL;
940 
941  tv = mono_tv_get();
942  im = mono_image_get();
943  if (tv != NULL && tv->tv_screen != NULL)    /* get roi from tv */
944    im = im_subim(im, tv_get_im_roi(tv));
945 
946  mono_edges_set(iso_canny(im, sigma, precision, low_edge_thres, up_edge_thres, len_thres, -grey_mid));
947  if (tv != NULL && tv->tv_screen != NULL)
948    im_free(im);
949  tv_edges_conn(tv, mono_edges_get());
950  tv_flush(tv);
951  mono_geom_set((List *) NULL);       /* now out of date */
952  
953                                           
954 }
955 
956 /* Applies the iso-canny algorithm to all images in the T1 image sequence
957    and builds a list of all the edge strings found.*/
958 static void test_3D_mins_proc(void)
959 {
960 
961  Tv     *tv = NULL;   
962  Sequence *seq = NULL;
963  List *lptr = NULL;
964 
965 
966 
967  if ((seq = seq_get_current()) == NULL)
968    {
969      format("No sequence available\n");
970      return;
971    }
972 
973  if (gm_er_list != NULL)
974    {
975      list_rm(gm_er_list, im_free);
976      gm_er_list = NULL;
977    }
978    
979 
980  for (lptr = seq->start; lptr != NULL; lptr = lptr->next)
981    {
982      Imrect *im = NULL;
983      Imrect *gm_er = NULL;
984    
985      im = lptr->to;
986 
987      gm_er = iso_canny(im, sigma, precision, low_edge_thres, up_edge_thres, len_thres, -grey_mid);
988      gm_er_list = list_addtoend(gm_er_list, link_alloc(gm_er, EDGERECT));
989    }
990 }
991 
992 /* For every GM/WM boundary edge position in the image calculates the 
993    difference in grey level between the assumed GM/WM
994   midpoint value and a bootstrapped estimate of the true boundary
995   value by assuming that the voxel values either side of the
996   boundary (using the 3D surface normal direction to the boundary) by
997   some distance defined using the "extent" parameter in the
998   dialogue box are in pure tissue. The average of the intensity values
999   at these positions, minus the GM/WM midpoint values gives the offset 
1000   values. The value of the offset is inserted into a histogram according to 
1001   the region the start position is in. */
1002 static void check_canny_proc(void)
1003 {
1004   Sequence *seq = NULL;
1005   List *lptr = NULL, *lptr_er = NULL, *lptr_gauss = NULL, *lptr_bound = NULL;
1006   int slice_no = 0;
1007   double sx=0, sy=0, sz=0;
1008   int old_interp_choice;
1009 
1010  if ((seq = seq_get_current()) == NULL)
1011    {
1012      format("No sequence available\n");
1013      return;
1014    }
1015 
1016  if (gm_er_list == NULL)
1017    {
1018      return;
1019    }
1020   
1021  if (bound_im_list == NULL)
1022    {
1023      return;
1024 
1025    }
1026  
1027  
1028  if (seq_gauss != NULL)
1029     {
1030       seq_rm(seq_gauss);
1031       seq_gauss = NULL;
1032       seq_set_current(seq);
1033     }
1034  seq_gauss = seq_copy(seq);
1035  seq_gauss = seq_gauss_3d(seq_gauss, 5.0, 0.5);
1036   
1037 
1038  tal_hist_init_temp();
1039 
1040  slice_no = seq->offset;
1041  sx = seq->dim[0];
1042  sy = seq->dim[1];
1043  sz = seq->dim[2];
1044 
1045  seq_slice_init(seq);
1046  old_interp_choice = seq_interp_choice(1);
1047  imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx, &imptrux);
1048  seq_init_interp(imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz);
1049 
1050  lptr_bound = bound_im_list;
1051 
1052  for (lptr = seq->start, lptr_er = gm_er_list, lptr_gauss = seq_gauss->start; lptr != NULL, lptr_er != NULL, lptr_gauss != NULL; lptr = lptr->next, lptr_er = lptr_er->next, lptr_gauss = lptr_gauss->next)
1053    {
1054      Imrect *im = NULL;
1055      Imrect *er_im = NULL;
1056      Imrect *gauss_im = NULL;
1057      Imrect *bound_im = NULL;
1058      List *string_list = NULL, *lptr2 = NULL;
1059 
1060      im = lptr->to;
1061      er_im = lptr_er->to; 
1062      bound_im = lptr_bound->to;
1063 
1064      string_list = (List*) prop_get(er_im->props, STRING);
1065       
1066       for(lptr2 = string_list; lptr2 != NULL; lptr2 = lptr2->next)
1067         {
1068           List *tptr = NULL;
1069           Tstring *edge_string = lptr2->to;
1070           /*List *edge_list = edge_string->start;*/
1071           /* doesnt this miss off the last edge on the string?*/
1072           for(tptr = edge_string->start; /*tptr != edge_string->end*/; tptr = tptr->next)
1073             {
1074               Vec3 orient, normal, centre, pos1, pos2, pos1b, pos2b; 
1075               Line3 *line = NULL;
1076               Edgel *edge = tptr->to;
1077               double mean, val1, val2;
1078               int i, j;
1079 
1080               
1081               centre = vec3(edge->pos.el[0]*sx, edge->pos.el[1]*sy, slice_no*sz); /*maybe*/
1082               i = tina_int(edge->pos.el[0]);
1083               j = tina_int(edge->pos.el[1]);
1084               if (((double)im_get_pixf(bound_im, j, i)) == 0.0)
1085                 {
1086                   orient = get_3D_orient(lptr_gauss, i, j); /*column, row*/
1087                   /*normal = vec3_perp(orient); YOU MORON */
1088                   pos1 = vec3_sum(centre, vec3_times(extent,orient));
1089                   pos1b = vec3(pos1.el[0]/sx, pos1.el[1]/sy, pos1.el[2]/sz);
1090                   pos2 = vec3_sum(centre, vec3_minus(vec3_times(extent, orient)));
1091                   pos2b = vec3(pos2.el[0]/sx, pos2.el[1]/sy, pos2.el[2]/sz);
1092                   
1093                   val1 = seq_interp(imptrs, pos1b);
1094                   val2 = seq_interp(imptrs, pos2b);
1095                   
1096                   mean = ((val1+val2)/2)-grey_mid;
1097                   line = line3_alloc(LINE3);
1098                   line->p1 = pos1;
1099                   line->length = mean;
1100               
1101                   tal_hist_fill(line, seq);
1102                 }
1103               if (tptr == edge_string->end)
1104                 break;
1105             }
1106         }
1107    
1108       lptr_bound = lptr_bound->next;
1109       slice_no++;
1110    }
1111  seq_interp_choice(old_interp_choice);
1112 
1113 }
1114 
1115 /* Plots the offset values for the given histogram number. Writes the 
1116    offset peak position to the tinaTool panel.*/
1117 static void canny_hist_plot_proc()
1118 {
1119  Tv *graph_tv = NULL;
1120   float median = 0.0, temp = 0.0;
1121   shistogram **shist_choice = NULL;
1122 
1123   graph_tv = imcalc_graph_tv_get();
1124 
1125   if (hist_no >=8 && hist_no <53)
1126     {
1127       switch (hist_choice)
1128         {
1129         case 0: 
1130           {
1131             shist_choice = shist_vec;
1132             break;
1133           }
1134         case 1:
1135           {
1136             shist_choice = shist_vec_left;
1137             break;
1138           }
1139         case 2:
1140           {
1141             shist_choice = shist_vec_right;
1142             break;
1143           }
1144         default:
1145           shist_choice = shist_vec;
1146         }
1147     }
1148 
1149   else
1150     {
1151       shist_choice = shist_vec;
1152     }
1153   graph_hfit(graph_tv, shist_choice[hist_no]);
1154   temp = hmax1(shist_choice[hist_no], &median);
1155   format("Histogram of: %s ", pLexRecords[hist_no]);
1156   format("median: %f\n", median);
1157   
1158   return;
1159 
1160 }
1161 
1162 /* Loops over all regions and plots the histograms of offset values per region.
1163    Also prints the peak offset positions to the tinaTool panel.*/
1164 static void canny_macro_proc(void)
1165 {
1166   
1167   Tv *graph_tv = NULL;
1168   float median = 0.0, temp = 0.0;
1169   int i, entries;
1170   shistogram *ph = NULL;
1171 
1172   graph_tv = imcalc_graph_tv_get();
1173 
1174   for (i=1; i<53; i++)
1175     {
1176       ph = shist_vec[i];
1177       graph_hfit(graph_tv, ph);
1178       temp = hmax1(ph, &median);
1179       entries = ph->entries;
1180       format("%d\t%s\t%f\t%d\n", i, pLexRecords[i], median, entries);
1181     }
1182 
1183   for (i=8; i<53; i++)
1184     {
1185       ph = shist_vec_left[i];
1186       graph_hfit(graph_tv, ph);
1187       temp = hmax1(ph, &median);
1188       entries = ph->entries;
1189       format("%d\tleft  %s\t%f\t%d\n", i, pLexRecords[i], median, entries);
1190     }
1191 
1192   for (i=8; i<53; i++)
1193     {
1194       ph = shist_vec_right[i];
1195       graph_hfit(graph_tv, ph);
1196       temp = hmax1(ph, &median);
1197       entries = ph->entries;
1198       format("%d\tright %s\t%f\t%d\n", i, pLexRecords[i], median, entries);
1199     }
1200 
1201   return;
1202  
1203 }
1204 
1205 
1206 /* Initialises the histograms and image volumes (GM prob and T1 images) and 
1207    runs the cortical thickness algorithm*/
1208 static void get_3d_min_proc()
1209 {
1210   List *lptr = NULL;
1211   List *lptr_gm = NULL;
1212   List *lptr_gauss = NULL;
1213   Sequence *seq = NULL;
1214   int curr, end, count;
1215   int old_interp_choice = 1;
1216 
1217   if ((seq = seq_get_current()) == NULL)
1218     {
1219       format("No sequence available\n");
1220       return;
1221     }
1222  
1223   if (seq_gauss != NULL)
1224     {
1225       seq_rm(seq_gauss);
1226       seq_gauss = NULL;
1227       seq_set_current(seq);
1228     }
1229   seq_gauss = (Sequence *)seq_copy(seq);
1230   seq_gauss = (Sequence *)seq_gauss_3d(seq_gauss, 5.0, 0.5);
1231 
1232   lptr_gm = gm_er_list;
1233   lptr_gauss = seq_gauss->start;
1234  
1235   end = get_seq_end();
1236   curr = get_slice_no() + seq->offset;
1237 
1238   if ((curr<seq->offset) || (curr>end))
1239     {
1240       format("out of sequence bound %d %d %d\n", seq->offset, curr, end);
1241       return;
1242     }
1243 
1244   if (gm_seq == NULL)
1245     return;
1246 
1247   tal_hist_init();
1248 
1249   seq_slice_init(gm_seq);
1250   old_interp_choice = seq_interp_choice(1);
1251   imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx, &imptrux);
1252   imptrs_gm_prob = parray_copy((char **)imptrs, imptrlz, imptrly, imptruz, imptruy);
1253   imptrs_gm = seq_slice_init(seq);
1254   imptrs = NULL; 
1255   seq_init_interp(imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz);
1256 
1257   for (lptr = seq->start, count = seq->offset; lptr != NULL; lptr = lptr->next)
1258     {
1259       /* need to set the current sequence here i think */
1260       set_seq_current_el(seq, lptr);
1261       set_seq_current_el(seq_gauss, lptr_gauss);
1262       find_min_dist_proc(lptr_gm, count); 
1263       lptr_gm = lptr_gm->next;
1264       lptr_gauss = lptr_gauss->next;
1265       count++;
1266     }
1267   format("count is: %d\n", count);
1268   seq_interp_choice(old_interp_choice);
1269 
1270     
1271   seq_rm(seq_gauss);
1272   seq_gauss = NULL;
1273   seq_set_current(seq);
1274   
1275   format("csf edge %d dip %d, wm edge %d dip %d\n", csf_edge, csf_dip, wm_edge, wm_dip);
1276   csf_edge = csf_dip = wm_edge = wm_dip = 0;
1277 }
1278 
1279 /* Stores the currently loaded sequence (ie, the GM probability volume) 
1280    to memory*/                            
1281 static void get_gm_vol_proc()
1282 {
1283   Sequence *seq = NULL;
1284   List *lptr = NULL;
1285 
1286   if ((seq = seq_get_current()) == NULL)
1287     {
1288       format("No sequence available\n");
1289       return;
1290     }
1291   
1292   gm_seq = seq_copy(seq);
1293 
1294 }
1295  
1296 
1297 /* For selecting colour according to thickness length*/
1298 static int choose_point_colour(float length)
1299 {
1300   int colour = 0;
1301   
1302   colour = tina_int(length*255/10);
1303   if (colour>255) colour = 255;
1304   
1305   return(colour);
1306 }
1307 
1308 /* function to create 3D representations of the inner and outer GM cortical
1309  surfaces . It doesn't work very well, and is excluded from the tool itself*/
1310 static void make_surf_proc(void)
1311 {
1312   List *lptr = NULL, *lptr2 = NULL, *tptr = NULL, *lptr_bi = NULL;
1313   int count = 0;
1314   double zscale = gm_seq->dim[2];
1315 
1316   for(lptr = gm_er_list, lptr_bi = bound_im_list; lptr->next != NULL, lptr_bi->next != NULL; lptr = lptr->next, lptr_bi = lptr_bi->next)
1317     {
1318       List *string_list = NULL;
1319       Imrect *im = lptr->to;
1320       Imrect *im2 = lptr->next->to;
1321       Imrect *im_bi = lptr_bi->to;
1322       int ux, uy, lx, ly;
1323       int temp_count = 0;
1324 
1325       ux = im->region->ux;
1326       uy = im->region->uy;
1327       lx = im->region->lx;
1328       ly = im->region->ly;
1329 
1330       string_list = (List*) prop_get(im->props, STRING);
1331       
1332       for(lptr2 = string_list; lptr2 != NULL; lptr2 = lptr2->next)
1333         {
1334           Tstring *edge_string = lptr2->to;
1335           /*List *edge_list = edge_string->start;*/
1336 
1337           for(tptr = edge_string->start; ; tptr = tptr->next)
1338             {
1339               /*tptr->to is an edgel*/
1340               Vec3 *ppt1 = NULL, *ppt2 = NULL, *ppt3 = NULL, *ppt4 = NULL; 
1341               Vec3 point1, point2, point3, point4;
1342               List *points = NULL;
1343               List *surf = NULL;
1344               int i, j;
1345               Vec2 temp_vec;
1346               Edgel *edge3b = NULL;
1347               Edgel *edge4 = NULL, *edge5 = NULL;
1348               Edge_conn *conn = NULL;
1349               Edgel *edge1 = tptr->to;
1350               Edgel *edge2 = tptr->next->to;
1351               double dist=300.0, temp_dist=0.0;
1352               Vec2 orient_vec, diff_vec; 
1353 
1354               if (((double)im_get_pixf(im_bi, tina_int(edge1->pos.el[1]), tina_int(edge1->pos.el[0]))) == 0.0)
1355                 {
1356                   point1 = vec3(edge1->pos.el[0], edge1->pos.el[1], count*zscale);
1357                   
1358                   diff_vec = vec2_diff(edge2->pos, edge1->pos);
1359                   orient_vec = vec2_of_polar(1.0, edge1->orient);
1360                   if (vec2_dot(diff_vec, orient_vec)>=0.0)
1361                     {
1362                       point2 = vec3(edge2->pos.el[0], edge2->pos.el[1], count*zscale);             
1363                     }
1364                   else
1365                     {
1366                       if (tptr->last != NULL)
1367                         {
1368                           edge2 = tptr->last->to;
1369                           point2 = vec3(edge2->pos.el[0], edge2->pos.el[1], count*zscale);      
1370                         }
1371                       else 
1372                         {
1373                           point2 = vec3(edge1->pos.el[0]+0.05, edge1->pos.el[1], count*zscale);
1374                         }
1375                     }
1376 
1377                   
1378 
1379                   for(i=lx; i<ux; i++)
1380                     {
1381                       for(j=ly; j<uy; j++)
1382                         {
1383                           Edgel *edge3 = er_get_edge(im2, j, i); 
1384                           if (edge3 != NULL)
1385                             {
1386                               if (((edge3->type & EDGE_GET_CONN_MASK) == EDGE_CONN)|| ((edge3->type & EDGE_GET_CONN_MASK) == EDGE_TERMIN))
1387                                 {
1388                                   if((temp_dist = vec2_dist(edge1->pos, edge3->pos)) < dist)
1389                                     {
1390                                       dist = temp_dist; 
1391                                       temp_vec = vec2(edge3->pos.el[0], edge3->pos.el[1]);
1392                                       edge3b = edge3;
1393                                     }
1394                                 }
1395                             }
1396                         }
1397                     }
1398                   
1399                   point3 = vec3(edge3b->pos.el[0], edge3b->pos.el[1], (count+1)*zscale);
1400                   conn = (Edge_conn *) prop_get(edge3b->props, EDGE_CONN);
1401                   edge4 = conn->c2;
1402                   edge5 = conn->c1;
1403                   if (edge4 == NULL || edge5 == NULL)
1404                     {
1405                       point4 = vec3(edge3b->pos.el[0]+0.05, edge3b->pos.el[1], (count+1)*zscale);
1406                     }
1407                   /* also need to have if edge4 xor edge5 are null */
1408                   else
1409                     {
1410                       double dist_edge4 = 0.0, dist_edge5 = 0.0;
1411                       dist_edge4 = vec2_dist(edge2->pos, edge4->pos);
1412                       dist_edge5 = vec2_dist(edge2->pos, edge5->pos);
1413                       if (dist_edge4 < dist_edge5)
1414                         {
1415                           point4 = vec3(edge4->pos.el[0], edge4->pos.el[1], (count+1)*zscale);
1416                         }
1417                       else
1418                         {
1419                           point4 = vec3(edge5->pos.el[0], edge5->pos.el[1], (count+1)*zscale);
1420                         }
1421                     }
1422                   
1423                   ppt1 = vec3_make(point1);
1424                   ppt2 = vec3_make(point2);
1425                   ppt3 = vec3_make(point3);
1426                   ppt4 = vec3_make(point4);
1427                   
1428                   points = list_addtoend(points, link_alloc(ppt1, VEC3));
1429                   points = list_addtoend(points, link_alloc(ppt2, VEC3));
1430                   points = list_addtoend(points, link_alloc(ppt3, VEC3));
1431                   points = list_addtoend(points, link_alloc(ppt4, VEC3));
1432                   
1433                   surf = link_alloc(points, LIST);
1434                   surf_list = list_addtoend(surf_list, surf);
1435                 }
1436               if(tptr == edge_string->end)
1437                 break;
1438             }
1439          
1440         }
1441       count++;
1442     }
1443 }
1444 
1445 /* Function to draw the inner and outer surfaces, determined using 
1446    make_surf_proc, in 3D. Excluded from main tool. */
1447 static void test_surf_proc(void)
1448 {
1449     List *lptr = NULL;
1450     Tv *tv = NULL;
1451     double z = 0.0;
1452    
1453     tv = threed_tv_get();
1454     z = tv_dist3(tv, vec3_zero()); /*already set in threed_fulldraw */
1455     tv_save_draw(tv);
1456     tv_set_zbuff(tv, z - 500.0, z + 500.0); /*already set in threed_fulldraw */
1457     tv_zbuff_shade_set(tv, spec_shade_inside_1);
1458     tv_zbuff_epsilon_set(0.02);
1459     tv_erase(tv);
1460 
1461     for (lptr=surf_list; lptr != NULL; lptr = lptr->next)
1462       {
1463         List *points = NULL;
1464         Vec3 *ppt1, *ppt2, *ppt3, *ppt4;
1465 
1466         points = lptr->to;
1467         ppt1 = points->to;
1468         ppt2 = points->next->to;
1469         ppt3 = points->next->next->to;
1470         ppt4 = points->next->next->next->to;
1471         
1472         tv_zbuff_ruled3(tv, *ppt1, *ppt2, *ppt3, *ppt4);
1473       }
1474     
1475   
1476 }
1477 
1478 /* For all GM/WM boundary edge voxels, plots a 3D representation, in the 
1479    threed tv, of hte line thicknesses starting at the GM/WM boundary.*/
1480 static void plot_3d_point_proc(void)
1481 {
1482   Tv *threed_tv = threed_tv_get();
1483   List * points_3d = NULL, *lptr = NULL, *pptr = NULL;
1484   List *lptr_er = NULL, *string_list = NULL;
1485   Graphic *gr = NULL;
1486   Imrect *er_im = NULL;
1487 
1488   for (lptr_er = gm_er_list; lptr_er != NULL; lptr_er = lptr_er->next)
1489     {
1490       er_im = lptr_er->to; 
1491       string_list = (List*) prop_get(er_im->props, STRING);
1492       
1493       for(lptr = string_list; lptr != NULL; lptr = lptr->next)
1494         {
1495           List *tptr = NULL;
1496           Tstring *edge_string = lptr->to;
1497           Edgel *edge_temp = NULL;
1498           
1499           
1500           for (tptr = edge_string->start; /*tptr != NULL*/; tptr=tptr->next)
1501             { 
1502               Line3 *line = NULL;
1503               Point3 *point = NULL;
1504               edge_temp = tptr->to;
1505               
1506               if ((line = prop_get(edge_temp->props, LINE3))!=NULL)
1507                 {
1508                   point = point3_make(line->p1, POINT3);
1509 
1510                   if ((gr = prop_get(point->props,GRAPHIC))==NULL) 
1511                     {
1512                       gr=(Graphic *)ralloc(sizeof(Graphic));
1513                       point->props = proplist_add(point->props, (void *)gr, GRAPHIC, rfree);
1514                     }
1515                   gr->colour = choose_point_colour(line->length);
1516                   
1517                   pptr = link_alloc(point, POINT3);
1518                   points_3d = list_addtoend(points_3d, pptr);
1519                                   
1520                 }
1521               if (tptr == edge_string->end)
1522                 break;
1523 
1524             }
1525         }
1526     }
1527   threed_geom_set(points_3d);
1528   reclist_list_draw(threed_tv, points_3d, POINT3, geom_col_draw, NULL);
1529   tv_repaint(threed_tv);
1530 
1531    
1532 }
1533 
1534 
1535 static void tv_choice_proc(int val)
1536 {
1537 tv_set_next(threed_tv());
1538 }
1539 
1540 
1541 
1542 static void tv_proc(void)
1543 {
1544   threed_tv_set(threed_tv_make());
1545   tv_set_next(threed_tv());
1546 }
1547 
1548 /* Plots the GM probability profile along the 3D surface normal through the GM
1549    starting at the chosen start point (defined using Line)*/
1550 static void plot_profile(Sequence *seq, Line3 *line)
1551 {
1552   Tv *tv_graph = NULL;
1553   Vec3 start_pos, end_pos, current_pos;
1554   double xscale, yscale, zscale;
1555   Vec3 line_start, line_end;
1556   Vec3 diff, unit;
1557   int length = 0;
1558   float *xdata = NULL, *ydata = NULL;
1559   int i=0;
1560   Imrect *prob_im = NULL;
1561 
1562   if ((tv_graph = imcalc_graph_tv_get()) == NULL)
1563     return;
1564 
1565   if (tv_graph->tv_screen == NULL)
1566     return;
1567   
1568   prob_im = get_seq_start_el(gm_seq)->to;
1569 
1570   xscale = seq->dim[0];
1571   yscale = seq->dim[1];
1572   zscale = seq->dim[2];
1573   line_start = line->p1;
1574   line_end = line->p2;
1575 
1576   start_pos = vec3(line_start.el[0]/xscale, line_start.el[1]/yscale, line_start.el[2]/zscale);
1577   end_pos = vec3(line_end.el[0]/xscale, line_end.el[1]/yscale, line_end.el[2]/zscale);
1578 
1579   /*right. start_pos and end_pos here are the voxel start and end, not the absolute values */
1580 
1581   diff = vec3_diff(end_pos, start_pos);
1582   unit = vec3_unit(diff);
1583 
1584 
1585   length = tina_int(vec3_mod(diff)) + 2;
1586 
1587   xdata = fvector_alloc(0, length);
1588   ydata = fvector_alloc(0, length);
1589 
1590   for (i=0; i<length; i++)
1591     {
1592       current_pos = vec3_sum(start_pos, vec3_times((double)i, unit));
1593       
1594       xdata[i] = (float)i;
1595       seq_voxel_vtype(prob_im->vtype);
1596       ydata[i] = seq_interp(imptrs_gm_prob, current_pos);
1597     }
1598   
1599   tv_repaint(tv_graph);
1600 
1601   plot(PL_INIT, PL_TV, tv_graph,
1602        PL_AXIS_COLOR, black,
1603        PL_TITLE, "Pixel Profile Plot",
1604        PL_COLOR, blue,
1605        PL_GRAPH_DATA, length, xdata, ydata,
1606        PL_PLOT,
1607        NULL);
1608   
1609   /*tv_repaint(tv_graph);*/
1610 
1611   fvector_free(xdata, 0);
1612   fvector_free(ydata, 0);
1613   return;
1614 
1615 }
1616 
1617 /* draws cross-sections of image imtensity values along a given line through 
1618  the GM.*/
1619 static void draw_profile_region(Sequence *seq, Line3 *line)
1620 {
1621   Vec3 line_start, line_end;
1622   Vec3 start_pos, end_pos, unit, norm, diff, temp_start, current_pos;
1623   double xscale, yscale, zscale;
1624   Imrect *im = NULL, *im2 = NULL;
1625   int length = 0;
1626   int i, j;
1627   float value = 0.0;
1628   Tv *imc_tv = NULL, *imc_tv2 = NULL;
1629   Imrect *grey_im = NULL;
1630 
1631   xscale = seq->dim[0];
1632   yscale = seq->dim[1];
1633   zscale = seq->dim[2];
1634   line_start = line->p1;
1635   line_end = line->p2;
1636   
1637   grey_im = get_seq_start_el(seq)->to;
1638 
1639   start_pos = vec3(line_start.el[0]/xscale, line_start.el[1]/yscale, line_start.el[2]/zscale);
1640   end_pos = vec3(line_end.el[0]/xscale, line_end.el[1]/yscale, line_end.el[2]/zscale);
1641   /* start_pos and end_pos are the voxel positions, not the absolute scaled pos */
1642   /*length = tina_int(vec3_dist(start_pos, end_pos))+3;*/
1643   length = extent;
1644 
1645   diff = vec3_diff(end_pos, start_pos);
1646   unit = vec3_unit(diff);
1647   /*norm = vec3_perp(unit);*/
1648   norm = vec3_unitcross(vec3_ex(), unit);
1649 
1650   im = im_alloc(6, length, NULL, float_v);
1651   im2 = im_alloc(6, length, NULL, float_v);
1652  
1653   for (i=-2; i<3; i++)
1654     {
1655       temp_start = vec3_sum(start_pos, vec3_times((double)i, norm));
1656       for(j=0; j<length; j++)
1657         {
1658           current_pos = vec3_sum(temp_start, vec3_times((double)j, unit));
1659           seq_voxel_vtype(grey_im->vtype);
1660           value = seq_interp(imptrs_gm, current_pos);
1661           im_put_pixf((double)value, im, i+2, j);
1662         }
1663 
1664     }
1665   /* unit increments of the unscaled unit vector */
1666 
1667   im_put_pixf(-2000, im, 5, 0);
1668   im_put_pixf(200, im, 5, 1);
1669 
1670   stack_push(im, IMRECT, im_free);
1671   
1672   format("unit vector: ");
1673   vec3_format(unit);
1674   format("normal in yz plane: ");
1675   vec3_format(norm);
1676   norm = vec3_unitcross(vec3_ey(), unit);
1677   format("normal in xz plane: ");
1678   vec3_format(norm);
1679   format("length: %f, %f \n", vec3_mod(diff), line->length); /*ist unscaled, 2nd scaled */
1680   format("start pos (voxel): ");
1681   vec3_format(start_pos);
1682 
1683   for (i=-2; i<3; i++)
1684     {
1685       temp_start = vec3_sum(start_pos, vec3_times((double)i, norm));
1686       for(j=0; j<length; j++)
1687         {
1688           current_pos = vec3_sum(temp_start, vec3_times((double)j, unit));
1689           seq_voxel_vtype(grey_im->vtype);
1690           value = seq_interp(imptrs_gm, current_pos);
1691           im_put_pixf((double)value, im2, i+2, j);
1692         }
1693       
1694     }
1695 
1696   im_put_pixf(-2000, im2, 5, 0);
1697   im_put_pixf(200, im2, 5, 1);
1698 
1699   stack_push(im2, IMRECT, im_free);
1700   imc_tv = imcalc_tv_get();
1701   imc_tv2 = imcal2_tv_get();
1702   tv_init(imc_tv);
1703   tv_init(imc_tv2);
1704   
1705   tv_color_set(imc_tv, red);
1706   tv_color_set(imc_tv2, red);
1707   tv_bigdot2(imc_tv, vec2(vec3_mod(diff), 2.5), 1);
1708   tv_bigdot2(imc_tv2, vec2(vec3_mod(diff), 2.5), 1);
1709 
1710 
1711   return;
1712 }
1713 
1714 /* calls code to plot graph profile and image cross-sections */
1715 static void profile_proc(Tv *tv, Ipos pos)
1716 {
1717   List *lptr = NULL;
1718   List *lptr_curr = NULL;
1719   List *lptr_gm = NULL;
1720   Sequence *seq = NULL;
1721   int count=0;
1722   Edgel *edge = NULL;
1723   Line3 *line = NULL;
1724   int x, y;
1725   Vec2    v = {Vec2_id};
1726 
1727   if((seq = seq_get_current()) == NULL)
1728     {
1729       format("no sequence\n");
1730       return;
1731     }
1732 
1733   if ((lptr_gm = gm_er_list) == NULL)
1734     {
1735       format("Calculate GM/WM Boundary first\n");
1736       return;
1737     }
1738 
1739   lptr_curr = get_seq_current_el(seq);
1740 
1741   for (lptr = seq->start, count = seq->offset/**/; lptr != lptr_curr; lptr = lptr->next)
1742     {
1743       lptr_gm = lptr_gm->next;
1744       count++;
1745     }
1746   
1747   format("count is: %d\n", count);
1748 
1749   v = tv_backproj2(tv, pos);
1750   x = tina_int(vec2_x(v));
1751   y = tina_int(vec2_y(v));
1752 
1753 
1754   if ((edge = er_get_edge(lptr_gm->to, y, x)) == NULL)
1755     {
1756       format("no edge here, returning\n");
1757       return;
1758     }
1759   
1760   if ((line = prop_get(edge->props, LINE3)) == NULL)
1761     {
1762       format("line is null\n");
1763       return;
1764     }
1765 
1766   plot_profile(seq, line);
1767   draw_profile_region(seq, line);
1768 
1769   return;
1770 
1771 }
1772 
1773 /* Plots a single regional histogram*/
1774 static void tal_hist_plot_proc(void)
1775 {
1776   Tv *graph_tv = NULL;
1777   float median = 0.0;
1778   shistogram **shist_choice = NULL;
1779   float upper_quart = 0.0;
1780   float lower_quart = 0.0;
1781   int entries = 0;
1782 
1783   graph_tv = imcalc_graph_tv_get();
1784 
1785   if (hist_no >=8 && hist_no <153)
1786     {
1787       switch (hist_choice)
1788         {
1789         case 0: 
1790           {
1791             shist_choice = shist_vec;
1792             break;
1793           }
1794         case 1:
1795           {
1796             shist_choice = shist_vec_left;
1797             break;
1798           }
1799         case 2:
1800           {
1801             shist_choice = shist_vec_right;
1802             break;
1803           }
1804         default:
1805           shist_choice = shist_vec;
1806         }
1807     }
1808 
1809   else
1810     {
1811       shist_choice = shist_vec;
1812     }
1813   graph_hfit(graph_tv, shist_choice[hist_no]);
1814   median = hmedian(shist_choice[hist_no]);
1815   lower_quart = hquartile(shist_choice[hist_no], 0.25);
1816   upper_quart = hquartile(shist_choice[hist_no], 0.75);
1817   entries = shist_choice[hist_no]->entries;
1818   format("%d\t%s\t%f\t%d\t%f\t%f\n", hist_no, pLexRecords[hist_no], median, entries, lower_quart, upper_quart);  
1819   
1820   return;
1821 }
1822 
1823 /* Plots histograms for all regions.*/
1824 static void median_macro_proc(void)
1825 {
1826   
1827   Tv *graph_tv = NULL;
1828   float median = 0.0;
1829   float upper_quart = 0.0;
1830   float lower_quart = 0.0;
1831   int i, entries;
1832   shistogram *ph = NULL;
1833 
1834   graph_tv = imcalc_graph_tv_get();
1835 
1836   for (i=1; i<53; i++)
1837     {
1838       ph = shist_vec[i];
1839       graph_hfit(graph_tv, ph);
1840       median = hmedian(ph);
1841       lower_quart = hquartile(ph, 0.25);
1842       upper_quart = hquartile(ph, 0.75);
1843       entries = ph->entries;
1844       format("%d\t%s\t%f\t%d\t%f\t%f\n", i, pLexRecords[i], median, entries, lower_quart, upper_quart);
1845     }
1846   i=90;
1847   ph = shist_vec[i]; /* Hard code brodmann area 10 */
1848   graph_hfit(graph_tv, ph);
1849   median = hmedian(ph);
1850   lower_quart = hquartile(ph, 0.25);
1851   upper_quart = hquartile(ph, 0.75);
1852   entries = ph->entries;
1853   format("%d\t%s\t%f\t%d\t%f\t%f\n", i, pLexRecords[i], median, entries, lower_quart, upper_quart);
1854 
1855   for (i=8; i<53; i++)
1856     {
1857       ph = shist_vec_left[i];
1858       graph_hfit(graph_tv, ph);
1859       median = hmedian(ph);
1860       lower_quart = hquartile(ph, 0.25);
1861       upper_quart = hquartile(ph, 0.75);
1862       entries = ph->entries;
1863       format("%d\tleft  %s\t%f\t%d\t%f\t%f\n", i, pLexRecords[i], median, entries, lower_quart, upper_quart);
1864     }
1865   i=90;
1866   ph = shist_vec_left[i];
1867   graph_hfit(graph_tv, ph);
1868   median = hmedian(ph);
1869   lower_quart = hquartile(ph, 0.25);
1870   upper_quart = hquartile(ph, 0.75);
1871   entries = ph->entries;
1872   format("%d\tleft  %s\t%f\t%d\t%f\t%f\n", i, pLexRecords[i], median, entries, lower_quart, upper_quart);
1873       
1874 
1875   for (i=8; i<53; i++)
1876     {
1877       ph = shist_vec_right[i];
1878       graph_hfit(graph_tv, ph);
1879       median = hmedian(ph);
1880       lower_quart = hquartile(ph, 0.25);
1881       upper_quart = hquartile(ph, 0.75);
1882       entries = ph->entries;
1883       format("%d\tright %s\t%f\t%d\t%f\t%f\n", i, pLexRecords[i], median, entries, lower_quart, upper_quart);
1884     }
1885   i=90;
1886   ph = shist_vec_right[i];
1887   graph_hfit(graph_tv, ph);
1888   median = hmedian(ph);
1889   lower_quart = hquartile(ph, 0.25);
1890   upper_quart = hquartile(ph, 0.75);
1891   entries = ph->entries;
1892   format("%d\tright %s\t%f\t%d\t%f\t%f\n", i, pLexRecords[i], median, entries, lower_quart, upper_quart);
1893   return;
1894   
1895 
1896 }
1897 
1898 
1899 /*  takes a given position on the tv and prints the grey level value.*/
1900 void    seq_grey(Tv * tv, Ipos pos)
1901 {
1902   Vec2    v = {Vec2_id};
1903   Vec2    tv_backproj2();
1904   Imrect *im = NULL;
1905   List *seq_curr = NULL;
1906   int     x, y, val;
1907  
1908   seq_curr = get_current_seq_current_el();
1909   im = seq_curr->to;
1910 
1911   v = tv_backproj2(tv, pos);     
1912   x = vec2_x(v);
1913   y = vec2_y(v);
1914   val = im_get_pix(im, y, x);
1915  
1916   format("pos %10d %10d:     grey level = %10d\n", x, y, val);
1917 }
1918 
1919 
1920 /* returns the median of 5 values*/
1921 static float get_median_of_five(float *len)
1922 {
1923   float order[5];
1924   float median = 0.0;
1925   int i,j,k;
1926 
1927   order[0] = len[0];
1928 
1929  /* go through the temp list */
1930   for (i=1; i<5; i++)
1931     {
1932 
1933       /* if next is greater than previous */
1934       if (len[i] >= order[i-1])
1935         {
1936           order[i] = len[i];
1937         }
1938       else
1939         {
1940           /* if next is less than zeroth */
1941           if (len[i] < order[0])
1942             {
1943               /*shuffle all up one and set to zeroth */
1944               
1945               for (j = 4; j>0; j--)
1946                 {
1947               order[j] = order[j-1];
1948                 }
1949               order[0] = len[i];
1950             }
1951           else
1952             {
1953               /* if next is greater than zeroth and less than previous */
1954               if ((len[i] >= order[0]) && (len[i] < order[i-1]))
1955                 {
1956                   /*search though to find the point at which it should be inserted, 
1957                     then shuffle all above up and insert */
1958                   for (k=0; k<i; k++)
1959                     {
1960                       if ((len[i] >= order[k]) && (len[i] < order[k+1]))
1961                         break;
1962                     }
1963                   
1964                   for (j = 4; j>(k+1); j--)
1965                     {
1966                       order[j] = order[j-1];
1967                     }
1968                   order[k+1] = len[i];
1969                   /* nice to see none of this actually worked in the first place? */
1970                 } 
1971             }
1972         }
1973     }
1974   median = order[2];
1975   return(median);
1976 
1977 }
1978 
1979 /* replaces the data for one edge structure according to the new length
1980    (median val) that it has been assigned.*/
1981 static void change_edge_line(Edgel *edge, float median_val)
1982 {
1983   
1984   Line3 *line = NULL;
1985   Graphic *gr = NULL;
1986   int zero_line = 0;
1987 
1988   if ((line = prop_get(edge->props, LINE3)) == NULL)
1989       return;
1990 
1991   if (line->length == 0.0)
1992     zero_line = 1;
1993 
1994   line->length = median_val;
1995   line->p2 = vec3_sum(line->p1, vec3_times((double)median_val, line->v));
1996   line->p = vec3_times(0.5, vec3_sum(line->p1, line->p2));
1997   
1998   if ((gr = prop_get(line->props,GRAPHIC))==NULL)
1999     {  
2000       gr=(Graphic *)ralloc(sizeof(Graphic));
2001       line->props = proplist_add(line->props, (void *)gr, GRAPHIC, rfree);
2002     }
2003   
2004   if (zero_line == 0) /*ie, if it was a line before */
2005     gr->colour = yellow;
2006   else 
2007     gr->colour = white; /* ie, if start equalled end never seems to fill these in. perhaps there aren't any*/
2008   /* that seems rather ominous doesn't it? yes ...*/
2009 
2010   if ((zero_line == 0) && (median_val == 0.0))
2011     {
2012       count_line_rm++;
2013     }
2014   if ((zero_line == 1) && (median_val != 0.0))
2015     {
2016       count_line++;
2017     }
2018   
2019 }
2020 
2021 /* Function to determine whether the thickness at a given voxel is spuriously 
2022  long and needs replacing with the median of the window of 5 edge voxels 
2023  centred on the voxel of interest. If so, replaces the length, recalculates 
2024  the regional histograms and colours the changed lengths yellow.*/
2025 static void local_support_proc(void)
2026 {
2027   List *lptr_er = NULL;
2028   Sequence *seq = NULL;
2029   
2030   if ((seq = seq_get_current()) == NULL)
2031     return;
2032   
2033   tal_hist_init();
2034 
2035   for (lptr_er = gm_er_list; lptr_er != NULL; lptr_er = lptr_er->next)
2036     {
2037       Imrect *er_im = NULL;
2038       List *string_list = NULL, *lptr2 = NULL;
2039       
2040       er_im = lptr_er->to; 
2041       
2042       string_list = (List*) prop_get(er_im->props, STRING);
2043       
2044       for(lptr2 = string_list; lptr2 != NULL; lptr2 = lptr2->next)
2045         {
2046           List *tptr = NULL, *list_start = NULL, *list_end = NULL, *tptr2 = NULL;
2047           Tstring *edge_string = lptr2->to;
2048           int string_length = str_length(edge_string);
2049           Edgel *edge_temp = NULL, *edge_temp2 = NULL;
2050           Line3 *line_temp = NULL;
2051           float len[5];
2052           float median_val=0.0, diff_val = 0.0;
2053           int x;
2054 
2055           list_start = edge_string->start->next->next->next; /*want to start on 4th*/
2056           list_end = edge_string->end->last;
2057          
2058           
2059           for (tptr = edge_string->start, x=0; x<5; tptr=tptr->next, x++)
2060             {
2061               edge_temp = tptr->to;
2062               if ((line_temp = prop_get(edge_temp->props, LINE3)) != NULL)
2063                 {
2064                   len[x] = line_temp->length;
2065                 }
2066               else 
2067                 len[x] = 0.0;
2068             }
2069           
2070           tptr2 = edge_string->start;
2071           for (x = 0; x<3; x++)
2072             {
2073               median_val = get_median_of_five(len);
2074               diff_val = len[x]-median_val;
2075               /*hfill1(mean_diff_hist, diff_val, 1.0);*/
2076               if (fabsf(diff_val)>(float)edge_limit)
2077                 {
2078                   edge_temp2 = tptr2->to;
2079                   change_edge_line(edge_temp2, median_val);
2080                 }
2081               tptr2 = tptr2->next;
2082             }
2083           for (tptr = list_start; tptr != list_end; tptr = tptr->next)
2084             {
2085               
2086               double len1=0.0, len2=0.0, len3=0.0, len4=0.0, len5=0.0;
2087               edge_temp = tptr->next->next->to;
2088               line_temp = prop_get(edge_temp->props, LINE3);
2089               
2090               len[0] = len[1];
2091               len[1] = len[2];
2092               len[2] = len[3];
2093               len[3] = len[4];
2094               if (line_temp != NULL)
2095                 len[4] = line_temp->length;
2096               else 
2097                 len[4] = 0.0;
2098               
2099               
2100               median_val = get_median_of_five(len);
2101               diff_val = len[2]-median_val;
2102               if (fabsf(diff_val)>(float)edge_limit)
2103                 {
2104                   edge_temp2 = tptr->to;
2105                   change_edge_line(edge_temp2, median_val);
2106                 }
2107               /*hfill1(mean_diff_hist, diff_val, 1.0);      */
2108             }
2109           for (x = 3; x<5; x++)
2110             {
2111               tptr2 = list_end;
2112               median_val = get_median_of_five(len);
2113               diff_val = len[x]-median_val;
2114               /*hfill1(mean_diff_hist, diff_val, 1.0);*/
2115               if (fabsf(diff_val)>(float)edge_limit)
2116                 {
2117                   edge_temp2 = tptr2->to;
2118                   change_edge_line(edge_temp2, median_val);
2119                 }
2120             }
2121           
2122         }
2123     }
2124   /*graph_hfit(graph_tv, mean_diff_hist);*/
2125   for (lptr_er = gm_er_list; lptr_er != NULL; lptr_er = lptr_er->next)
2126     {
2127       Imrect *er_im = NULL;
2128       List *string_list = NULL, *lptr2 = NULL;
2129       
2130       er_im = lptr_er->to; 
2131       
2132       string_list = (List*) prop_get(er_im->props, STRING);
2133       
2134       for(lptr2 = string_list; lptr2 != NULL; lptr2 = lptr2->next)
2135         {
2136           List *tptr = NULL;
2137           Tstring *edge_string = lptr2->to;
2138 
2139           for (tptr = edge_string->start; ; tptr=tptr->next)
2140             {
2141               Edgel *edge_temp = NULL;
2142               Line3 *line = NULL;
2143 
2144               edge_temp = tptr->to;
2145               
2146               if ((line = prop_get(edge_temp->props, LINE3)) != NULL)
2147                 {
2148                   if (line->length != 0.0)
2149                     tal_hist_fill(line, seq);
2150                 }
2151 
2152 
2153               if (tptr == edge_string->end)
2154                 break;
2155             }
2156 
2157         }
2158     }  
2159   format("line_rm: %d, line_add: %d\n", count_line_rm, count_line);
2160 }
2161 
2162 
2163 Tv_mouse profile_mouse(void)
2164 {
2165   return (mouse_define(MOUSE_NAME, "profile mouse",
2166                        LEFT_NAME, "profile",
2167                        LEFT_DOWN, profile_proc, 
2168                        RIGHT_NAME, "grey", 
2169                        RIGHT_DOWN, seq_grey, NULL));
2170 }
2171 
2172 static void seq_mouse_proc(Tv_mouse(*func) ( /* ??? */ ))
2173 {
2174   tv_set_mouse(seq_tv(), (*func) ());
2175   (void) tv_set_activity(seq_tv(), MOUSE);
2176 }
2177 
2178 
2179 static void hist_choice_proc(int val)
2180 {
2181   hist_choice = val;
2182 
2183 }
2184 
2185 static void cort_params()
2186 {
2187    static void *dialog = NULL;
2188 
2189     if (dialog)
2190     {
2191         tw_show_dialog(dialog);
2192         return;
2193     }
2194 
2195     dialog = tw_dialog("Cortical Thickness Parameters");
2196 
2197     tw_label("Canny Parameters");
2198     tw_newrow();
2199     psi = (void *)tw_fglobal("Sigma: ", &sigma, 10);
2200     ppr = (void *)tw_fglobal("Precision: ", &precision, 10);
2201     tw_newrow();
2202     plt = (void *)tw_fglobal("Lower Thresh: ", &low_edge_thres, 10);
2203     put = (void *)tw_fglobal("Upper Thresh: ", &up_edge_thres, 10);
2204     tw_newrow();
2205     plen = (void *)tw_iglobal("Length Thresh: ", &len_thres, 10);
2206     tw_newrow();
2207     tw_label("Skull Segmentation Param: (3.0 -> 4.0)");
2208     tw_newrow();
2209     pssp = (void *)tw_fglobal("Skull Seg:", &skull_seg, 10);
2210     tw_newrow();
2211     tw_label("Cortical Thickness Parameters");
2212     tw_newrow();
2213     pgr = (void *)tw_fglobal("Grey/white midpoint: ", &grey_mid, 10);
2214     pgm = (void *)tw_fglobal("Grey Mean: ", &grey_mean, 10);
2215     tw_newrow();
2216     pext = (void *)tw_iglobal("Search extent: ", &extent, 10);
2217     pedl = (void *)tw_fglobal("Edge Limit: ", &edge_limit, 10);
2218     tw_newrow();
2219     pfrac = (void *)tw_fglobal("Fraction: ", &fraction, 10);
2220     /*pmind = (void *)tw_fglobal("Min diff scale: ", &min_diff_scale, 10);*/
2221     
2222     tw_end_dialog();
2223 }
2224 
2225 
2226 void    cortical_thickness_tool(int x, int y)
2227 {
2228  
2229    static void *tool = NULL;
2230 
2231     if (tool)
2232     {
2233         tw_show_tool(tool);
2234         return;
2235     }
2236     
2237     tv_proc();
2238 
2239     tool = (void *)tw_tool("Cortical Thickness Tool", x, y);
2240 
2241     
2242     tw_choice("Tv : ", tv_choice_proc, 0, "3D rep", NULL);
2243     tw_menubar("Mouse ", "seq", "profile", seq_mouse_proc, profile_mouse, NULL, NULL);
2244     tw_button("Cort params", cort_params, NULL);
2245     tw_label("         ");
2246     tw_help_button("Cortical");
2247     tw_newrow();
2248     tw_button("Test Boundary", test_canny_proc, NULL); 
2249     tw_newrow();
2250     tw_label("Cortical Thickness Algorithm");
2251     tw_newrow();
2252     tw_button("Get GM vol", get_gm_vol_proc, NULL);
2253     tw_label(" -> ");
2254     tw_button("Remove Skull", get_boundary_proc, NULL);
2255     tw_label(" -> ");
2256     tw_button("Get GM/WM", test_3D_mins_proc, NULL);
2257     tw_label(" -> ");
2258     tw_button("Get Thickness", get_3d_min_proc, NULL);
2259     tw_newrow();
2260     tw_button("Apply Median Filter", local_support_proc, NULL);
2261     tw_newrow();
2262     tw_choice("Hist: ", hist_choice_proc, 0, "Whole", "Left", "Right", NULL);
2263     ptalno = (void *)tw_iglobal("Hist no: ", &hist_no, 10);
2264     tw_newrow();
2265     tw_label("Plot Regional Thickness Histograms");
2266     tw_newrow();
2267     tw_button("Single Hist.", tal_hist_plot_proc, NULL);
2268     tw_button("All Hists.", median_macro_proc, NULL);
2269     tw_button("Remove Hists.", rm_hist_proc, NULL);
2270     tw_newrow();
2271     tw_label("Graphics Results");
2272     tw_newrow();
2273     tw_button("2D Overlay", overlay_line_proc, NULL);
2274     tw_button("3D Lines", plot_3d_point_proc, NULL);
2275     /*tw_button("Make 3D surface", make_surf_proc, NULL);
2276       tw_button("Draw 3D surface", test_surf_proc, NULL);*/
2277     tw_newrow();
2278     tw_label("GM/WM Boundary Placement Checks");
2279     tw_newrow();
2280     tw_button("Calculate Offset", check_canny_proc, NULL);
2281     tw_button("Plot Single Offset", canny_hist_plot_proc, NULL); 
2282     tw_button("Plot All Offsets", canny_macro_proc, NULL);
2283   
2284     tw_end_tool();
2285 
2286 
2287 }
2288 

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