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

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

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