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

Linux Cross Reference
Tina6/tina-tools/tinatool/tlmedical/tlmedOctants_tool.c

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

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  *
  6  * Redistribution and use in source and binary forms, with or without modification,
  7  * are permitted provided that the following conditions are met:
  8  *
  9  *   . Redistributions of source code must retain the above copyright notice,
 10  *     this list of conditions and the following disclaimer.
 11  *
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation
 14  *     and/or other materials provided with the distribution.
 15  *
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this
 18  *     software without specific prior written permission.
 19  *
 20  *
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :
 37  * Date    :
 38  * Version :
 39  * CVS Id  :
 40  *
 41  * Notes :
 42  *
 43  *********
 44 */
 45 
 46 #include "tlmedOctants_tool.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #    include <config.h>
 50 #endif
 51 
 52 
 53 #include <math.h>
 54 #include <stdlib.h>
 55 #include <stdio.h>
 56 #include <sys/param.h>
 57 #include <tina/sys/sysDef.h>
 58 #include <tina/sys/sysPro.h>
 59 #include <tina/math/mathDef.h>
 60 #include <tina/math/mathPro.h>
 61 #include <tina/file/filePro.h>
 62 #include <tina/image/imgPro.h>
 63 #include <tina/geometry/geomPro.h>
 64 
 65 #include <tinatool/gphx/gphxDef.h>
 66 #include <tinatool/gphx/gphxPro.h>
 67 #include <tinatool/wdgts/wdgtsDef.h>
 68 #include <tinatool/wdgts/wdgtsPro.h>
 69 #include <tinatool/draw/drawPro.h>
 70 #include <tinatool/tlbase/tlbasePro.h>
 71 #include <tinatool/tlmedical/tlmed_SegPro.h>
 72 #include <tinatool/tlmedical/tlmed_CoregPro.h>
 73 #include <tinatool/tlmedical/tlmedOctants_alloc.h>
 74 #include <tinatool/tlmedical/tlmedOctants_csfcount12.h>
 75 
 76 
 77 /* Files */
 78 
 79 static void *p_dir=NULL, *p_data_file = NULL, *p_thresh_file = NULL, *p_tvs_file = NULL,
 80 *p_output_file = NULL;
 81 static char directory_name[MAXPATHLEN], data_list_name[MAXPATHLEN], mean_model_name[MAXPATHLEN],
 82 tvs_file_name[MAXPATHLEN],
 83 output_file_name[MAXPATHLEN];
 84 
 85 static FILE *p_file_global=NULL, *p_file_global2=NULL; 
 86 /* Global so that the file does not have to be passed through loop_on_data_list */
 87 
 88 /* Stuff for coreg checker */
 89 
 90 static void *g_data_list=NULL;
 91 static char g_norm_brain_path[512], g_rough_air_path[512];
 92 static int g_no_files=0, g_norm_no_frames=0;
 93 
 94 /* General parameters */
 95 
 96 static int knear=0;
 97 static int fixmeanage_flag = 0;
 98 static double global_mean_age=65.0;
 99 static int scale_opt_flag=0;
100 
101 /* Parameters for box-in-brain threshold determination */
102 
103 static int box_lx = 70, box_ly=80, box_lz=14, box_ux=180, box_uy=150, box_uz=20;        /* Tuned box params */
104 static int box2_lx = 70, box2_ly=80, box2_lz=27, box2_ux=180, box2_uy=150, box2_uz=32;  /* Tuned box params */
105 
106 /* Parameters for merging and re-splitting the sequence */
107 
108 static int seqtv_lx, seqtv_ux, seqtv_ly, seqtv_uy;
109 static int seq_lx, seq_ux, seq_ly, seq_uy, seq_lz, seq_uz;
110 
111 
112 /******************************** Handling functions for external tools ********************/
113 
114 
115 static void seq_one_proc(void)
116 {
117     Sequence *seq = seq_get_current();
118     Imrect *seq_im=NULL, *new_im=NULL;
119     Imregion *roi=NULL;
120     List *lptr=NULL, *nptr=NULL;
121     int no_frames=0, i, j, k, height, width, end_frame;
122     float pix;
123     
124     nptr=seq->start;
125     seq_im = (Imrect *)nptr->to;
126     roi = seq_im->region;
127     for(lptr=nptr; lptr!=NULL; lptr=lptr->next)
128     {
129         no_frames++;
130     }
131     seq_lx = roi->lx;
132     seq_ux = roi->ux;
133     seq_ly = roi->ly;
134     seq_uy = roi->uy;
135     seq_lz = 0;
136     seq_uz = no_frames;
137     
138     height = (seq_uy-seq_ly)*no_frames;
139     width = seq_ux-seq_lx;
140     
141     roi = roi_alloc(0, 0, width, height);
142     new_im = im_alloc(height, width, roi, float_v);
143     
144     k=0;
145     for(lptr=nptr; lptr!=NULL; lptr=lptr->next)
146     {
147         seq_im = lptr->to;
148         for(i=seq_ly; i<seq_uy; i++)
149         {
150             for(j=seq_lx; j<seq_ux; j++)
151             {
152                 pix = (float)im_get_pix(seq_im, i, j);
153                 im_put_pixf(pix, new_im, k, (j-seq_lx));
154             }
155             k++;
156         }
157     }
158 
159     dd_list_rm(seq->start, im_free);
160     seq->start = NULL;
161     seq->current = NULL;
162     seq->end = NULL;
163     seq_frame_insert(new_im, seq);
164     end_frame = get_end_frame(seq);
165     set_seq_end(end_frame);
166     tv_init(seq_tv_get());
167     seq_show_current(seq);
168     seq_param_update(seq);
169     seq_init_funcs_run(); 
170     tv_set_roi(seq_tv_get(), roi);
171     rfree(roi);
172 }
173 
174 
175 static void one_seq_proc(void)
176 {
177     Sequence *seq = seq_get_current();
178     List *lptr=NULL;
179     Imrect *seq_im=NULL, *new_im=NULL;
180     Imregion *seq_roi = NULL, *roi=NULL;
181     int i, j, k, down, across, end_frame;
182     float pix;
183     
184     lptr = seq->start;
185     seq_im = im_copy((Imrect *)lptr->to);
186     seq_roi = seq_im->region;
187     dd_list_rm(seq->start, im_free);
188     seq->start = NULL;
189     seq->current = NULL;
190     seq->end = NULL;
191     
192     down=0;
193     roi = roi_alloc(seq_lx, seq_ly, seq_ux, seq_uy);
194     for(k=seq_lz; k<seq_uz; k++)
195     {
196         new_im = im_alloc((seq_uy-seq_ly), (seq_ux-seq_lx), roi, float_v);
197         for(i=seq_ly; i<seq_uy; i++)
198         {
199             if(i<seqtv_ly||i>=seqtv_uy)
200             {
201                 for(j=seq_lx; j<seq_ux; j++)
202                 {
203                     im_put_pixf(0.0, new_im, i, j);
204                 }
205             }
206             else
207             {
208                 across=0;
209                 for(j=seq_lx; j<seq_ux; j++)
210                 {
211                     if(j<seqtv_lx||j>=seqtv_ux)
212                     {
213                         im_put_pixf(0.0, new_im, i, j);
214                     }
215                     else
216                     {
217                         pix = im_get_pixf(seq_im, down, across);
218                         im_put_pixf(pix, new_im, i, j);
219                         across++;
220                     }
221                 }
222                 down++;
223             }
224         }
225         stack_push(im_copy(new_im), IMRECT, im_free);
226         seq_frame_insert(new_im, seq);
227         end_frame = get_end_frame(seq);
228         set_seq_end(end_frame);
229     }
230 
231     tv_init(seq_tv_get());
232     seq_show_current(seq);
233     seq_param_update(seq);
234     seq_init_funcs_run(); 
235     tv_set_roi(seq_tv_get(), roi);
236     rfree(roi);
237 }
238 
239 
240 static void read_tset(FILE *fp, Octants *octants)
241 {
242         double *a=NULL;
243 
244         a = (double *)octants->a;
245 
246         fscanf(fp," %s %d %f", octants->name, &(octants->diag), &(octants->age));
247         fscanf(fp," --a %lf -- b %lf --c %lf", &(a[0]), &(a[1]), &(a[2]));
248         fscanf(fp," -+a %lf -+ b %lf -+c %lf", &(a[3]), &(a[4]), &(a[5]));
249         fscanf(fp," +-a %lf +- b %lf +-c %lf", &(a[6]), &(a[7]), &(a[8]));
250         fscanf(fp," ++a %lf ++ b %lf ++c %lf", &(a[9]), &(a[10]), &(a[11]));
251         fscanf(fp,"  x mm %lf, y mm %lf, z mm %lf", &(a[12]), &(a[13]), &(a[14]));
252 }
253 
254 
255 static double dist(double *scales, int i, int person, void **tsets)
256 {
257         double distance=1.0;
258         Octants *octants_i=NULL, *octants_p=NULL;
259 
260         octants_i = (Octants *)tsets[i];
261         octants_p = (Octants *)tsets[person];
262 
263         distance *= exp(-(octants_i->middle - octants_p->middle - octants_i->front + octants_p->front)*
264                         (octants_i->middle - octants_p->middle - octants_i->front + octants_p->front)
265                         /(4.0*scales[0]*scales[0]));
266 
267         distance *= exp(-(octants_i->middle - octants_p->middle - octants_i->back + octants_p->back)*
268                         (octants_i->middle - octants_p->middle - octants_i->back + octants_p->back)
269                         /(4.0*scales[1]*scales[1]));
270 
271         distance *= exp(-(octants_i->front - octants_p->front + octants_i->middle - octants_p->middle 
272         + octants_i->back - octants_p->back)* (octants_i->front - octants_p->front + octants_i->middle 
273         - octants_p->middle + octants_i->back - octants_p->back)
274                         /(6.0*scales[2]*scales[2]));
275 
276         distance *= exp(-(octants_i->left - octants_i->right - octants_p->left + octants_p->right)
277                         *(octants_i->left - octants_i->right - octants_p->left + octants_p->right)
278                         /(4.0*scales[3]*scales[3]));
279 
280         distance *= exp(-(octants_i->top - octants_i->bot - octants_p->top + octants_p->bot)
281                         *(octants_i->top - octants_i->bot - octants_p->top + octants_p->bot)
282                         /(4.0*scales[4]*scales[4]));
283 
284         return(distance);
285 }
286 
287 
288 static int near_prob(double *scales, int person, int num_dat, void  **tsets)
289 {
290         int i, best_gfac;
291         double gfac[12];
292         double total=0.0;
293         Octants *octants=NULL;
294 
295         for (i=0;i<12;i++) gfac[i] = 0.0;
296         for (i=0;i<num_dat;i++)
297         {
298                 if(i!=person)
299                 {
300                         octants = (Octants *)tsets[i];  
301                         gfac[((int)(octants->diag))] += dist(scales, i, person, tsets); 
302                 }
303         }
304 
305         best_gfac = 11;
306         total = gfac[1] + gfac[2] + gfac[3] + gfac[4] + gfac[5] + gfac[6] + gfac[7] + gfac[8] + gfac[9] 
307         + gfac[10] + gfac[11];
308         for (i=1;i<12;i++) gfac[i] /= total;
309         for (i=1;i<12;i++)
310         {
311                 if (gfac[i]>gfac[best_gfac]) best_gfac = i; 
312         }
313         if (gfac[best_gfac]>0.00) return(best_gfac);
314         else return(0);
315 } 
316 
317 
318 static void near_prob_pab(double *scales, int person, int num_dat, void  **tsets, double
319 *best_gfac_prob, int *best_gfac_i)
320 {
321         int i, best_gfac;
322         double gfac[12];
323         double total=0.0;
324         Octants *octants=NULL;
325 
326         for (i=0;i<12;i++) gfac[i] = 0.0;
327         for (i=0;i<num_dat;i++)
328         {
329                 if(i!=person)
330                 {
331                         octants = (Octants *)tsets[i];  
332                         gfac[((int)(octants->diag))] += dist(scales, i, person, tsets); 
333                 }
334         }
335 
336         for(i=1; i<12; i++) total += gfac[i];
337         for (i=1;i<12;i++) gfac[i] /= total;
338         best_gfac = 11;
339 
340         for (i=1;i<12;i++)
341         {
342                 if (gfac[i]>gfac[best_gfac]) best_gfac = i;
343         }
344         if(gfac[best_gfac]>0.00)
345         {
346                 *best_gfac_prob = gfac[best_gfac];
347                 *best_gfac_i = best_gfac;
348                 return;
349         }
350         else
351         {
352                 *best_gfac_prob = 0.00;
353                 *best_gfac_i = 0;
354                 return;
355         }
356 }
357 
358 
359 static void xgobi_output(void **tsets, int knear)
360 {
361         FILE *f_output_dat_ptr = NULL, *f_output_col_ptr = NULL, *f_output_row_ptr = NULL, *f_output_v12_ptr = NULL;
362         char col_output_pathname[MAXPATHLEN], dat_output_pathname[MAXPATHLEN], row_output_pathname[MAXPATHLEN], v12_output_pathname[MAXPATHLEN];
363         int i, j, diag;
364         Octants *octants=NULL;
365         double *a=NULL, age, volume, total, tiv;
366 
367         (void) strip_spaces(output_file_name);
368         (void) strip_spaces(directory_name);
369         (void) string_append(col_output_pathname, output_file_name, ".colors", NULL);
370         (void) string_append(row_output_pathname, output_file_name, ".row", NULL);
371         (void) string_append(dat_output_pathname, output_file_name, ".dat", NULL);
372         (void) string_append(v12_output_pathname, output_file_name, ".v12_norm", NULL);
373 
374         f_output_col_ptr = open_file_for_writing(directory_name, col_output_pathname);
375         if(f_output_col_ptr==NULL) return;
376         f_output_row_ptr = open_file_for_writing(directory_name, row_output_pathname);
377         if(f_output_row_ptr==NULL) return;
378         f_output_dat_ptr = open_file_for_writing(directory_name, dat_output_pathname);
379         if(f_output_dat_ptr==NULL) return;
380         f_output_v12_ptr = open_file_for_writing(directory_name, v12_output_pathname);
381         if(f_output_v12_ptr==NULL) return;
382 
383         for(i=0; i<knear; i++)
384         {
385                 octants = (Octants *)tsets[i];
386                 a = (double *)octants->a;
387                 age = (double)(octants->age);
388                 diag = (int)(octants->diag);
389                 volume = a[12]*a[13]*a[14];
390                 tiv = 0.538*volume;
391 
392                 fprintf(f_output_v12_ptr, "%d \t", diag);
393                 fprintf(f_output_v12_ptr, "%f \t", age);
394                 total = 0.0;
395                 for(j=0; j<12; j++)
396                 {
397                         total = total+a[j];
398                         fprintf(f_output_v12_ptr, "%f \t", a[j]);
399                 }
400                 fprintf(f_output_v12_ptr, " %f ", (total*tiv));
401                 fprintf(f_output_v12_ptr, " %f ", volume);
402                 fprintf(f_output_v12_ptr, " %f ", tiv);
403                 fprintf(f_output_v12_ptr, " %f \n", total);
404 
405 
406                 fprintf(f_output_dat_ptr, " %d %f %f %f %f %f %f\n", (octants->diag), 
407                 ((octants->middle - octants->front)/sqrt(2.0)),
408                                                 ((octants->middle - octants->back)/sqrt(2.0)),
409                                                 ((octants->front + octants->middle + octants->back)/sqrt(3.0)), 
410                                                 ((octants->left - octants->right)/sqrt(2.0)), 
411                                                 ((octants->top - octants->bot)/sqrt(2.0)),
412                                                 (octants->age/100.0));
413 
414                 
415                 /* New diagnosis numbers for optimisation */
416 
417                 if (octants->diag == 1) fprintf(f_output_col_ptr, "SkyBlue\n");                 /* normal */
418                 if (octants->diag == 2) fprintf(f_output_col_ptr, "Orange\n");                  /* alz */
419                 if (octants->diag == 3) fprintf(f_output_col_ptr, "Red\n");                     /* vasc */
420                 if (octants->diag == 4) fprintf(f_output_col_ptr, "Green\n");                   /* ftd */
421                 if (octants->diag == 5) fprintf(f_output_col_ptr, "Yellow\n");                  /* Lewy */
422                 if (octants->diag == 6) fprintf(f_output_col_ptr, "Blue\n");                    /* Schiz */
423                 if (octants->diag == 7) fprintf(f_output_col_ptr, "White\n");                   /* Possible Alz */
424                 if (octants->diag == 8) fprintf(f_output_col_ptr, "White\n");                   /* Possible Vasc */
425                 if (octants->diag == 9) fprintf(f_output_col_ptr, "Grey\n");                    /* Undefined */
426 
427                 fprintf(f_output_row_ptr, "%s\n", octants->name);
428         }
429         fclose(f_output_dat_ptr);
430         fclose(f_output_col_ptr);
431         fclose(f_output_row_ptr);
432         fclose(f_output_v12_ptr);
433 }
434 
435 
436 static void pab_hprint(FILE *p_file, shistogram *confusion_matrix)
437 {
438         /* Print out the values in the confusion matrix */
439 
440         int xbins, ybins, i, j;
441         float **array;
442         
443         xbins = (int)confusion_matrix->xbins;
444         ybins = (int)confusion_matrix->ybins;
445         array = (float **)confusion_matrix->array;
446 
447         fprintf(p_file, "Confusion Matrix\n\n");
448 
449         for(i=0; i<ybins; i++)
450         {
451                 for(j=0; j<xbins; j++)
452                 {
453                         fprintf(p_file, "%f  ", ((float)array[i][j]));
454                 }
455                 fprintf(p_file, "\n");
456         }
457 }
458 
459 
460 static double scale_select_cost_function(int n, double *scales, void *tset)
461 {
462         Octants *octants=NULL;
463         double cost=0, best_gfac_prob;
464         int i, best_gfac_i;
465         void **tset_here=NULL;
466         
467         tset_here = (void **)tset;
468 
469         for (i=0;i<knear;i++)
470         {
471                 octants = (Octants *)tset_here[i];
472                 near_prob_pab(scales, i, knear, tset_here, &best_gfac_prob, &best_gfac_i);
473                 if (octants->diag != near_prob(scales,i,knear, tset_here))
474                 {
475                         cost += best_gfac_prob;
476                 }
477         }
478 
479         return cost;
480 }
481 
482 
483 static void octants_results_classifier(void)
484 {
485         int i, j, best_gfac_i;
486         void **tset=NULL;
487         double *a=NULL, volume = 0, fit_a, fit_b, fit_c, fit_d, total=0, *averages=NULL, tiv, fit_pat, fit_mean;
488         double *scales=NULL, *lambda=NULL, mean_age, best_gfac_prob;
489         float age=0;
490         FILE *p_output_file = NULL;
491         FILE *p_input_file = NULL;
492         Octants *octants=NULL;
493         shistogram *confusion_matrix=NULL, *new_confusion_matrix=NULL;
494 
495         p_output_file = open_file_for_writing(directory_name, output_file_name);
496         p_input_file = open_file_for_reading(directory_name, tvs_file_name);
497         if(p_input_file==NULL||p_output_file==NULL) return;
498 
499         scales = dvector_alloc(0, 5);
500 
501         scales[0] =0.010;
502         scales[1] =0.007;
503         scales[2] =0.028;
504         scales[3] =0.007;
505         scales[4] =0.017;
506 
507         averages = dvector_alloc(0, 12);
508 
509         for(j=0; j<12; j++)
510         {
511                 averages[j] = 0;
512         }
513         total=0;
514 
515         tset = pvector_alloc(0, knear);
516         confusion_matrix = hbook2(19,"confusion matrix\0", 0.0, 11.0, 11, 0.0, 11.0, 11);
517         new_confusion_matrix = hbook2(20,"new confusion matrix\0", 0.0, 11.0, 11, 0.0, 11.0, 11);
518 
519         /* Calculate Mean Age */
520 
521         mean_age = 0;
522         for(i=0; i<knear; i++)
523         {
524                 octants = octants_alloc();
525                 tset[i] = (void *)octants;
526                 read_tset(p_input_file, octants);
527                 mean_age += octants->age;
528         }
529         mean_age=mean_age/knear;
530         if(fixmeanage_flag==1) mean_age = global_mean_age;
531 
532         /* Calculate averages */
533 
534         for (i=0;i<knear;i++)
535         {
536                 octants = (Octants *)tset[i];
537                 a = (double *)octants->a;
538                 age = octants->age;
539                 volume = a[12]*a[13]*a[14];
540 
541                 for(j=0; j<12; j++)
542                 {
543                         total += a[j];
544                         averages[j] += a[j];
545                 }
546         }
547         for(j=0; j<12; j++)
548         {
549                 averages[j] = 100*averages[j]/total;
550         }
551 
552         format("The averages are:\n");
553         for(j=0; j<12; j++)
554         {
555                 format("b[%d]=%f\n", j, averages[j]);
556         }
557 
558         /* Weibull age correction: see Tina Memo no. 2004-002 for details */
559 
560 
561         fit_a = 0.0727;
562         fit_b = 1.4970;
563         fit_c = 3.0885;
564 
565         for (i=0;i<knear;i++)
566         {
567                 octants = (Octants *)tset[i];
568 
569                 age = (octants->age);
570                 a = (double *)octants->a;
571                 volume = a[12]*a[13]*a[14];
572                 tiv = 0.538*volume;
573 
574                 for(j=0; j<12; j++)
575                 {
576                         a[j] = a[j] / tiv;
577                         fit_d = 0.009872 * averages[j];
578                         fit_pat = fit_d*(fit_a + (1-fit_a)*(1-exp(-pow((age/(fit_b*100.0)), fit_c))));
579                         fit_mean = fit_d*(fit_a + (1-fit_a)*(1-exp(-pow((mean_age/(fit_b*100.0)),fit_c))));
580                         a[j] = a[j]*(fit_mean/fit_pat);
581                 }
582 
583                 octants->front = sqrt(a[0] + a[3]  +a[6] +a[9]);
584                 octants->middle = sqrt(a[4] + a[10] +a[1] +a[7]);
585                 octants->back  = sqrt(a[2] + a[5]  +a[8] +a[11]);
586                 octants->left =  sqrt(a[0] + a[1]  +a[2] +a[3] +a[4] +a[5]);
587                 octants->right = sqrt(a[6] + a[7]  +a[8] +a[9] +a[10] +a[11]);
588                 octants->top = sqrt(a[0] + a[1] + a[2] + a[6] +a[7] +a[8]);
589                 octants->bot = sqrt(a[3] +a[4] +a[5] +a[9] +a[10] + a[11]);
590         }
591 
592 
593         if(scale_opt_flag==1)
594         {
595                 lambda = dvector_alloc(0, 5);
596                 for(i=0; i<5; i++) lambda[i] = 0.005;
597                 simplexmin2(5, scales, lambda, scale_select_cost_function, tset, 1.0e-8, NULL);
598 
599                 for(i=0; i<5; i++) format("scales[%d] = %f \n", i, scales[i]);
600                 dvector_free(lambda, 0);
601         }
602 
603         for (i=0;i<knear;i++)
604         {
605                 octants = (Octants *)tset[i];
606                 near_prob_pab(scales, i, knear, tset, &best_gfac_prob, &best_gfac_i);
607                 hfill2(new_confusion_matrix, (float)(octants->diag+0.5), ((float)best_gfac_i + 0.5),
608                 ((float)best_gfac_prob));
609                 hfill2(confusion_matrix, (float)(octants->diag+0.5), 
610                 (float)(near_prob(scales, i, knear, tset) + 0.5), 1.0);
611                 if (octants->diag != near_prob(scales,i,knear, tset))
612                         fprintf(p_output_file,"%s %d %d \n", octants->name, octants->diag, 
613                         near_prob(scales, i, knear, tset));
614         }
615 
616         pab_hprint(p_output_file, confusion_matrix);
617         pab_hprint(p_output_file, new_confusion_matrix);
618 
619         hprint(p_output_file, confusion_matrix);
620         fclose(p_input_file);
621         fclose(p_output_file);
622         xgobi_output(tset, knear);
623 
624         for(i=0; i<knear; i++)
625         {
626                 octants = (Octants *)tset[i];
627                 octants = octants_free(octants);
628         }
629         pvector_free(tset, 0);
630         tset=NULL;
631 
632         hfree(confusion_matrix);
633         confusion_matrix=NULL;
634         hfree(new_confusion_matrix);
635         new_confusion_matrix=NULL;
636         dvector_free(averages, 0);
637         dvector_free(scales, 0);
638 }
639 
640 
641 static List *stack_to_list_backwards(void)
642 {
643         /* Remove images from stack and put them into a list, but backwards */
644 
645         Imrect *im=NULL;
646         List *lptr=NULL, *nptr=NULL;
647         int type;
648 
649         while((im=(Imrect *)stack_pop(&type))!=NULL)
650         {
651                 lptr = link_alloc(im, type);
652                 if(nptr!=NULL) nptr->last = lptr;
653                 lptr->next = nptr;
654                 nptr = lptr;
655         }
656 
657         tv_repaint(imcalc_tv_get());
658         tv_repaint(imcal2_tv_get());
659         return lptr;
660 }
661 
662 
663 /************************* Backend functions for loop_on_data_list ***************************************/
664 
665 
666 static void octants_coregref_loader(char *volume_path, int volume_frames)
667 {
668         double coreg_centrex=128.0, coreg_centrey=128.0, coreg_centrez=24.5;
669 
670         /* Load and latch the first coreg volume */
671 
672         seq_load_wrapper(volume_path, 0, 1, volume_frames, 1);
673         coreg_centrez = 0.5+(volume_frames/2);
674         reset_centre(coreg_centrex, coreg_centrey, coreg_centrez);
675         coreg_reset_threshbordscale(-2500, 5, 3);
676         coreg_zoom_proc();
677         coreg_par_choice_notify_proc(7);
678         coreg_check_proc(1);
679         coreg_reslice_choice_proc(1);
680         seq_str_rm_proc();
681 }
682 
683 static List *norm_mask_prepare(char *norm_brain_path, int norm_no_frames, char *image_path, int no_frames)
684 {
685         double coreg_centrex=128.0, coreg_centrey=128.0, coreg_centrez=24.5;
686         List *lptr=NULL;
687         char air_path[512], norm_mask_path[512];
688 
689         /* Load the normal brain masks and reslice them in the coordinate system of the volumes in the data list */
690 
691         (void) string_append(norm_mask_path, norm_brain_path, "_bincut", NULL);
692         seq_load_wrapper(norm_mask_path, 0, 1, norm_no_frames, 1);
693         coreg_centrez = 0.5;
694         coreg_reslice_choice_proc(0);
695         reset_centre(coreg_centrex, coreg_centrey, coreg_centrez);
696 
697         (void) string_append(air_path, image_path, ".tair_mi", NULL);
698         coreg_air_load(air_path);
699         coreg_batch_reslice();
700         lptr = stack_to_list_backwards();
701         stack_clear();
702         seq_str_rm_proc();
703         return(lptr);
704 }
705 
706 
707 static void bulk_coreg_backend(int no_frames, char *image_path, char *pat_name, float pat_age, int pat_diag,
708                         char *norm_brain_path, char *rough_air_path, int norm_no_frames, int i)
709 {
710         char output_path[512];
711         struct air16 air1;
712         int j;
713 
714         /* Coregister the normal brain to the volumes in the data list, and output tair files */
715 
716         octants_coregref_loader(image_path, no_frames);
717         seq_load_wrapper(norm_brain_path, 0, 1, norm_no_frames, 1);
718         coreg_air_load(rough_air_path); 
719         for(j=0; j<3; j++) mui_coreg_auto_proc();
720         get_coreg_trans(&air1);
721         (void) string_append(output_path, image_path, ".tair_mi", NULL);
722         write_air16(output_path, 0, &air1);
723         seq_str_rm_proc();
724         coreg_check_proc(0);
725 }
726 
727 
728 static Imrect *im_subim_pab(Imrect *im)
729 {
730         Imrect *im_new=NULL, *im_smooth=NULL;
731         int i, j;
732         float pix;
733         
734         im_new = im_alloc((box_uy-box_ly), (box_ux-box_lx), NULL, im->vtype);
735         for(i=0; i<box_uy; i++)
736         {
737                 for(j=0; j<box_ux; j++)
738                 {
739                         pix = im_get_pixf(im, (i+box_ly), (j+box_lx));
740                         im_put_pixf(pix, im_new, i, j);
741                 }
742         }
743         
744         im_smooth = im_tsmooth(im_new);
745         im_free(im_new);
746         return(im_smooth);
747 }
748 
749 
750 static Imrect *im_subim_pab2(Imrect *im)
751 {
752         Imrect *im_new=NULL, *im_smooth=NULL;
753         int i, j;
754         float pix;
755         
756         im_new = im_alloc((box2_uy-box2_ly), (box2_ux-box2_lx), NULL, im->vtype);
757         for(i=0; i<box2_uy; i++)
758         {
759                 for(j=0; j<box2_ux; j++)
760                 {
761                         pix = im_get_pixf(im, (i+box2_ly), (j+box2_lx));
762                         im_put_pixf(pix, im_new, i, j);
763                 }
764         }
765         
766         im_smooth = im_tsmooth(im_new);
767         im_free(im_new);
768         return(im_smooth);
769 }
770 
771 
772 static void bulk_em_seg_proc_prep_images(int no_frames, char *image_path, char *pat_name, float pat_age, int pat_diag,
773                         char *norm_brain_path, char *rough_air_path, int norm_no_frames, int i)
774 {
775         double coreg_centrex=128.0, coreg_centrey=128.0, coreg_centrez=24.5;
776         char air_path[512];
777         List *lptr=NULL, *nptr=NULL;
778         Imrect *im=NULL;
779         int j=0;
780 
781         /* Load the volumes from the air file in the coordinate frame of the normal brain,
782         and determine CSF thresholds */
783 
784         octants_coregref_loader(norm_brain_path, norm_no_frames);
785         seq_load_wrapper(image_path, 0, 1, no_frames, 1);
786         (void) string_append(air_path, image_path, ".tair_mi", NULL);
787         coreg_air_load(air_path);
788         coreg_air_invert_proc();
789         coreg_centrez = 0.5;
790         reset_centre(coreg_centrex, coreg_centrey, coreg_centrez);
791         coreg_batch_reslice();
792         seq_str_rm_proc();
793         coreg_check_proc(0);
794         octants_coregref_loader(norm_brain_path, norm_no_frames);
795         stack_seq(NULL);
796         coreg_redraw();
797 
798         lptr = get_current_seq_start_el();
799         for(nptr=lptr; nptr!=NULL; nptr=nptr->next)
800         {
801                 j++;
802                 im = (Imrect *)nptr->to;
803                 if(j>box_lz && j<box_uz)
804                 {
805                         stack_push(im_subim_pab(im), IMRECT, im_free);
806                 }
807                 else if(j>box2_lz && j<box2_uz)
808                 {
809                         stack_push(im_subim_pab2(im), IMRECT, im_free);
810                 }
811         }
812         seq_str_rm_proc();
813         coreg_check_proc(0);
814         stack_seq(NULL);
815 }
816 
817 
818 void bulk_em_seg_proc_backend(int no_frames, char *image_path, char *pat_name, float pat_age, int pat_diag,
819                         char *norm_brain_path, char *rough_air_path, int norm_no_frames, int i)
820 {
821         char model_input_path[512], model_output_path[512];
822         Mixmodel *model=NULL;
823 
824         /* Load the middle of the sequence, perform the segmentation, and write out the model */
825 
826         string_append(model_input_path, directory_name, "/", mean_model_name, NULL);
827         bulk_em_seg_proc_prep_images(no_frames, image_path, pat_name, pat_age, pat_diag, norm_brain_path, 
828                                         rough_air_path, norm_no_frames, i);
829         model = input_mixmodel(model_input_path); 
830         if(model==NULL)
831         {
832                 format("Model file for %s not found", pat_name);
833                 return;
834         }
835 
836         em_set_model(model);
837         seq_one_proc();
838         run_edgeseg_proc();
839         (void) string_append(model_output_path, image_path, "_EM_model", NULL);
840         model = em_get_model();
841         output_mixmodel(model, model_output_path);
842         seq_str_rm_proc();
843 }
844 
845 
846 void output_model_parameters_proc_backend(int no_frames, char *image_path, char *pat_name,
847                 float pat_age, int pat_diag, char *norm_brain_path, char *rough_air_path,
848                 int norm_no_frames, int i)
849 {
850         char model_input_path[512];
851         Mixmodel *model=NULL;
852         double tmp;
853         int j, k;
854 
855         (void) string_append(model_input_path, image_path, "_EM_model", NULL);
856         model = input_mixmodel(model_input_path);
857 
858         seq_load_wrapper(image_path, 25, 1, 30, 1);
859 
860         for (i = 0; i < model->nmix; i++)
861         {
862                 for (j = 0; j < model->ndim; j++)
863                 {
864                         fprintf(p_file_global, " %f ", model->vectors[i][j]);
865                 }
866                 fprintf(p_file_global, " %f ", model->volume[i]);
867         }
868 
869         for (k = 0; k < model->nmix; k++)
870         {
871                 for (i = 0; i < model->ndim; i++)
872                 {
873                         for (j = 0; j < model->ndim; j++)
874                         {
875                                 tmp = model->cov[k][i][j];
876                                 if (tmp != 0.0)
877                                         tmp /= sqrt(fabs(tmp));
878                                 fprintf(p_file_global, " %f ", tmp);
879                         }
880                 }
881         }
882 
883         for (i = 0; i < model->nmix; i++)
884         {
885                 for (j = 0; j < model->nmix; j++)
886                 {
887                         fprintf(p_file_global, " %f  ", model->par_dens[i][j]);
888                 }
889         }
890 
891         if (model->gradscale  != NULL)
892         {
893                 for (i = 0; i < model->nmix; i++)
894                 {
895                         for (j = 0; j < model->nmix; j++)
896                         {
897                                 fprintf(p_file_global, "  %f  ", model->gradscale[i][j]);
898                         }
899                 } 
900         }
901 
902         fprintf(p_file_global, "\n");
903         fprintf(p_file_global2, "%s \n", pat_name);
904 
905         seq_str_rm_proc();
906 }
907 
908 
909 void bulk_em_seg2_proc_backend(int no_frames, char *image_path, char *pat_name, float pat_age, int pat_diag,
910                         char *norm_brain_path, char *rough_air_path, int norm_no_frames, int i)
911 {
912         Sequence *seq=NULL;
913         Imrect *im=NULL, *prod=NULL, *bin_prod=NULL, *mask=NULL;
914         Mixmodel *model=NULL;
915         char model_input_path[512], *tissue_type=NULL, seq_output_path[512];
916         List *lptr=NULL, *nptr=NULL, *mptr=NULL, *optr=NULL;
917         double xscale, yscale, zscale;
918         int j, filename_length;
919 
920         tissue_type = em_get_str();
921         strcpy(tissue_type, "CSF\0");
922 
923         octants_coregref_loader(image_path, no_frames);
924         lptr = norm_mask_prepare(norm_brain_path, norm_no_frames, image_path, no_frames);
925 
926         (void) string_append(model_input_path, image_path, "_EM_model", NULL);
927         model = input_mixmodel(model_input_path);
928         format("Name=%s", image_path);
929 
930         for(j=0; j<no_frames; j++)
931         {
932                 model = input_mixmodel(model_input_path);
933                 em_set_model(model);
934                 seq_load_wrapper(image_path, j, 1, j, 1);
935                 seq_get_iscales(&xscale, &yscale, &zscale);
936                 if((seq = seq_get_current())==NULL)
937                 {
938                         format("No frames wrong for %s\n", pat_name);
939                         break;
940                 }
941                 seq_type_choice_proc(3);
942 
943                 em_estep_grad_proc();
944                 em_prior_proc();
945                 em_estep_grad_proc();
946                 em_prob_class_plot_proc();
947                 seq_str_rm_proc();
948         }
949         
950         stack_seq(NULL);
951         seq=seq_get_current();
952         seq->dim[0]=xscale;
953         seq->dim[1]=yscale;
954         seq->dim[2]=zscale;
955 
956         (void) string_append(seq_output_path, image_path, "_maskedEM", NULL);
957         filename_length = strlen(seq_output_path);
958         seq->filename = cvector_copy(seq_output_path, 0, (filename_length+1));
959 
960         nptr=seq->start;
961         optr = lptr;
962 
963         for(mptr=nptr; mptr!=NULL; mptr=mptr->next)
964         {
965                 im = (Imrect *)mptr->to;
966                 mask = (Imrect *)optr->to;
967                 optr = optr->next;
968                 prod = im_prod(mask, im);
969                 bin_prod = im_bthresh(0.5, prod);
970                 mptr->to = (Imrect *)im_cast(bin_prod, uchar_v);
971                 im_free(bin_prod);
972                 im_free(prod);
973                 im_free(im);
974         }
975 
976         seq_save(seq, seq_output_path, NULL, NULL, NULL, NULL, 1);
977 
978         seq_str_rm_proc();
979         list_rm(lptr, im_free);
980         lptr=NULL;
981         coreg_check_proc(0);
982 }
983 
984 
985 static void bulk_classifier_proc_backend(int no_frames, char *image_path, char *pat_name,
986         float pat_age, int pat_diag,char *norm_brain_path, char *rough_air_path, int norm_no_frames, int i)
987 {
988         char masked_path[512], air_path[512];
989         List *center_point_list=NULL, *plist=NULL;
990         Match *nptr=NULL;
991         float *weight_vector=NULL, *brain_dim=NULL, *count=NULL;
992         int j=0;
993 
994         /* Read in the data file, calculate the octants, write the full tvs file for all data sets
995            and do the classification */
996 
997         brain_dim = fvector_alloc(0, 3);
998         count = fvector_alloc(0, 12);
999 
1000         (void) string_append(masked_path, image_path, "_maskedEM", NULL);
1001         (void) string_append(air_path, image_path, ".tair_mi", NULL);
1002         octants_coregref_loader(masked_path, no_frames);
1003         coreg_reslice_choice_proc(0);
1004         seq_load_wrapper(masked_path, 0, 1, no_frames, 1);
1005         coreg_air_load(air_path);
1006         coreg_air_invert_proc();
1007         center_point_list = csf_count12_octants_proc(brain_dim, NULL);
1008 
1009         for(plist = center_point_list; plist!=NULL; plist=plist->next)
1010         {
1011                 nptr = (Match *)plist->to;
1012                 weight_vector = nptr->to2;
1013                 count[j] = weight_vector[0];
1014                 count[(j+6)] = weight_vector[1];
1015                 j++;
1016         }
1017 
1018         fprintf(p_file_global, " %s %d %f\n", pat_name, pat_diag, pat_age);
1019         fprintf(p_file_global, " --a %f --b %f  --c %f \n",count[0],count[1],count[2]);
1020         fprintf(p_file_global, " -+a %f -+b %f  -+c %f \n",count[3],count[4],count[5]);
1021         fprintf(p_file_global, " +-a %f +-b %f  +-c %f \n",count[6],count[7],count[8]);
1022         fprintf(p_file_global, " ++a %f ++b %f  ++c %f \n",count[9],count[10],count[11]);
1023         fprintf(p_file_global, " x mm %f, y mm %f, z mm %f \n", brain_dim[0],brain_dim[1],brain_dim[2]);
1024 
1025         center_point_list = center_point_list_destroy(center_point_list);
1026         fvector_free(brain_dim, 0);
1027         fvector_free(count, 0);
1028         seq_str_rm_proc();
1029         coreg_check_proc(0);
1030 }
1031 
1032 
1033 static int read_data_list_header(FILE *p_data_list, int *no_files, char *rough_air_path, int
1034 *norm_no_frames, char *norm_brain_path)
1035 {
1036         char *error=NULL, *error2=NULL;
1037 
1038         fscanf(p_data_list,"%d \n", no_files);
1039         error = fgets(rough_air_path, 512, p_data_list);
1040         fscanf(p_data_list,"%d", norm_no_frames);
1041         error2 = fgets(norm_brain_path, 512, p_data_list);
1042         (void) strip_spaces(rough_air_path);
1043         (void) strip_spaces(norm_brain_path);
1044         
1045         if(error==NULL||error2==NULL) return 0;
1046         return 1;
1047 }
1048 
1049 
1050 static int read_data_list_loop(FILE *p_data_list, int *no_frames, char *image_path, char *pat_name,
1051 int *pat_diag, float *pat_age)
1052 {
1053         char *error=NULL;
1054 
1055         fscanf(p_data_list,"%d", no_frames);
1056         error = fgets(image_path, 512, p_data_list);
1057         (void) strip_spaces(image_path);
1058         fscanf(p_data_list," %s %d %f", pat_name, pat_diag, pat_age);
1059         
1060         if(error==NULL) return 0;
1061         return 1;
1062 }
1063 
1064 
1065 static void loop_on_data_list(void(*function_to_call)(int no_frames, char *image_path, char
1066                 *pat_name, float pat_age, int pat_diag, char *norm_brain_path, char *rough_air_path, 
1067                 int norm_no_frames, int counter))
1068 {
1069         /* Generic function to open data list, load its data, then loop over the image volumes listed,  */
1070         /* calling the lower level functions.  The format for the data list is:                         */
1071         /*                                                                                              */
1072         /* [No of records]                                                                              */
1073         /* [path to the rough tair file] (for initial approximate coregistration)                       */
1074         /* [No of slice in reference brain] [path to reference brain]                                   */
1075         /* [No of slice in image volume] [path to image volume]                                         */
1076         /* [patient name] [patient disease code] [patient age]                                          */
1077         /*                                                                                              */
1078         /* So, there should be {no_records} sets of the last two lines                                  */
1079 
1080         FILE *p_data_list=NULL;
1081         int no_files=0, no_frames = 0, norm_no_frames=0, pat_diag, i, code;
1082         char rough_air_path[512], norm_brain_path[512], image_path[512], pat_name[30];
1083         float pat_age;
1084 
1085         p_data_list = open_file_for_reading(directory_name, data_list_name);
1086         if(p_data_list==NULL) return;
1087 
1088         code = read_data_list_header(p_data_list, &(no_files), rough_air_path, 
1089                 &(norm_no_frames), norm_brain_path);
1090         if(code==0) return;
1091 
1092         for(i=0; i<no_files; i++)
1093         {
1094                 code = read_data_list_loop(p_data_list, &(no_frames), image_path, pat_name,
1095                 &(pat_diag), &(pat_age));
1096                 if(code==0) break;
1097                 function_to_call(no_frames, image_path, pat_name, pat_age, pat_diag, norm_brain_path, 
1098                 rough_air_path, norm_no_frames, i);
1099         }
1100 
1101         fclose(p_data_list);
1102 }
1103 
1104 
1105 /************************* Coreg checker *****************************************************************/
1106 
1107 
1108 static void coreg_checker_proc(void)
1109 {
1110         int code;
1111 
1112         g_data_list = open_file_for_reading(directory_name, data_list_name);
1113         if(g_data_list==NULL) return;
1114         code = read_data_list_header(g_data_list, &(g_no_files), g_rough_air_path, &(g_norm_no_frames), g_norm_brain_path);
1115         if(code==0) format("Error reading file header\n");
1116 }
1117 
1118 
1119 static void coreg_checker_next(void)
1120 {
1121         int no_frames, pat_diag, code;
1122         char image_path[512], pat_name[30], air_path[512];
1123         float pat_age;
1124         Sequence *seq=NULL;
1125 
1126         if(g_data_list==NULL)
1127         {
1128                 error("Data list not open\n", warning);
1129                 return;
1130         }
1131         seq=seq_get_current();
1132         if(seq!=NULL)
1133         {
1134                 seq_str_rm_proc();
1135                 coreg_check_proc(0);
1136         }
1137         code = read_data_list_loop(g_data_list, &(no_frames), image_path, pat_name, &(pat_diag), &(pat_age));
1138         if(code==0)
1139         {
1140                 format("Error reding file record\n");
1141                 return;
1142         }
1143 
1144         coreg_check_proc(0);
1145         (void) string_append(air_path, image_path, ".tair_mi", NULL);
1146         octants_coregref_loader(image_path, no_frames);
1147         seq_load_wrapper(g_norm_brain_path, 0, 1, g_norm_no_frames, 1);
1148         coreg_air_load(air_path);
1149         format(air_path);
1150 }
1151 
1152 
1153 static void coreg_checker_end(void)
1154 {
1155         Sequence *seq=NULL;
1156 
1157         seq=seq_get_current();
1158         if(seq!=NULL)
1159         {
1160                 seq_str_rm_proc();
1161                 coreg_check_proc(0);
1162         }
1163         if(g_data_list!=NULL) fclose(g_data_list);
1164 }
1165 
1166 
1167 /********************** Box-in-brain display *************************/
1168 
1169 
1170 static void display_boxmask_proc(void)
1171 {
1172         /* Write something crap so that you can see the box in the brain */
1173 
1174         Sequence *seq=seq_get_current();
1175         List *lptr = get_current_seq_start_el(), *nptr=NULL;
1176         int i=0, j=0, k=0, pix;
1177         Imrect *im=NULL, *im2=NULL;
1178         Imregion *roi2=NULL;
1179 
1180         for(nptr=lptr; nptr!=NULL; nptr=nptr->next)
1181         {
1182                 im = (Imrect *)nptr->to;
1183                 roi2 = im->region;
1184 
1185                 im2 = im_alloc((roi2->ux-roi2->lx), (roi2->uy-roi2->ly), NULL, float_v);
1186                 for(j = roi2->ly; j<roi2->uy; j++)
1187                 {
1188                         for(k=roi2->lx; k<roi2->ux; k++)
1189                         {
1190                                 pix = im_get_pixf(im, j, k);
1191                                 if(box_lx<k && k<box_ux && box_ly<j && j<box_uy && i>box_lz && i<box_uz)
1192                                 {
1193                                         im_put_pixf(pix, im2, j, k);
1194                                 }
1195                                 else if(box2_lx<k && k<box2_ux && box2_ly<j && j<box2_uy && i>box2_lz && i<box2_uz)
1196                                 {
1197                                         im_put_pixf(pix, im2, j, k);
1198                                 }
1199                                 else
1200                                 {
1201                                         im_put_pixf(0.0, im2, j, k);
1202                                 }
1203                         }
1204                 }
1205                 im_free(im);
1206                 nptr->to = im2;
1207 
1208                 i++;
1209         }
1210         coreg_slice_init(seq, seq_dim_to_vec3());
1211         coreg_coords_init();
1212 }
1213 
1214 
1215 /********************** Wrappers for calls to loop_on_data_list *************************/
1216 
1217 
1218 static void bulk_coreg_proc(void)
1219 {
1220         loop_on_data_list(bulk_coreg_backend);
1221 }
1222 
1223 
1224 static void image_prep_proc(void)
1225 {
1226         loop_on_data_list(bulk_em_seg_proc_prep_images);
1227 }
1228 
1229 
1230 static void bulk_em_seg_proc(void)
1231 {
1232         loop_on_data_list(bulk_em_seg_proc_backend);
1233 }
1234 
1235 
1236 static void bulk_em_seg2_proc(void)
1237 {
1238         loop_on_data_list(bulk_em_seg2_proc_backend);
1239 }
1240 
1241 
1242 static void bulk_classifier_proc(void)
1243 {
1244         p_file_global = open_file_for_writing(directory_name, tvs_file_name);
1245         if(p_file_global==NULL)
1246         {
1247                 format("Could not open TVS file\n");
1248                 return;
1249         }
1250         loop_on_data_list(bulk_classifier_proc_backend);
1251         fclose(p_file_global);
1252 }
1253 
1254 
1255 static void output_model_parameters_proc(void)
1256 {
1257         char output_file_name_dat[512], output_file_name_row[512];
1258 
1259         (void) string_append(output_file_name_dat, output_file_name, ".dat", NULL);
1260         (void) string_append(output_file_name_row, output_file_name, ".row", NULL);
1261 
1262         p_file_global2 = open_file_for_writing(directory_name, output_file_name_row);
1263         p_file_global = open_file_for_writing(directory_name, output_file_name_dat);
1264         if(p_file_global==NULL||p_file_global2==NULL)
1265         {
1266                 format("Could not open output file\n");
1267                 return;
1268         }
1269         loop_on_data_list(output_model_parameters_proc_backend);
1270         fclose(p_file_global);
1271         fclose(p_file_global2);
1272 }
1273 
1274 
1275 /*********************** Utility functions *********************************************************/
1276 
1277 
1278 static void tvs_to_v12_proc(void)
1279 {
1280         /* Input a tiv file and a tvs file and combine them */
1281 
1282         int i, j, diag;
1283         FILE *p_input_tvs=NULL, *p_output_file=NULL;
1284         Octants *octants=NULL;
1285         double *a=NULL, age, volume, total, tiv;
1286 
1287         p_input_tvs = open_file_for_reading(directory_name, tvs_file_name);
1288         p_output_file = open_file_for_writing(directory_name, output_file_name);
1289         if(p_input_tvs==NULL) return;
1290 
1291 
1292         octants = octants_alloc();
1293         for (i=0;i<knear;i++)
1294         {
1295                 read_tset(p_input_tvs, octants);
1296 
1297                 a = (double *)octants->a;
1298                 age = (double)(octants->age);
1299                 diag = (int)(octants->diag);
1300                 volume = a[12]*a[13]*a[14];
1301                 tiv = 0.538*volume;
1302 
1303                 fprintf(p_output_file, "%d \t", diag);
1304                 fprintf(p_output_file, "%f \t", age);
1305                 total = 0.0;
1306                 for(j=0; j<12; j++)
1307                 {
1308                         total = total+a[j];
1309                         fprintf(p_output_file, "%f \t", a[j]/tiv);
1310                 }
1311                 fprintf(p_output_file, " %f ", total);
1312                 fprintf(p_output_file, " %f ", volume);
1313                 fprintf(p_output_file, " %f ", tiv);
1314                 fprintf(p_output_file, " %f \n", (total/tiv));
1315         }
1316         octants = octants_free(octants);
1317         fclose(p_input_tvs);
1318         fclose(p_output_file);
1319 }
1320 
1321 
1322 /*********************** Choices *********************************************************/
1323 
1324 
1325 static void fixmeanage_choice_proc(int value)
1326 {
1327         fixmeanage_flag = value;
1328 }
1329 
1330 
1331 static void scale_opt_choice_proc(int value)
1332 {
1333         scale_opt_flag = value;
1334 }
1335 
1336 /*********************** Tool and dialog creation ****************************************/
1337 
1338 
1339 static void segmentation_dialog_proc(void)
1340 {
1341         static void *dialog=NULL;
1342 
1343         if(dialog)
1344         {
1345                 tw_show_dialog(dialog);
1346                 return;
1347         }
1348         dialog = tw_dialog("Segmentation Parameters");
1349 
1350         tw_iglobal("Box lx:", &box_lx, 10);
1351         tw_iglobal("    ux:", &box_ux, 10);
1352         tw_newrow();
1353         tw_iglobal("Box ly:", &box_ly, 10);
1354         tw_iglobal("    uy:", &box_uy, 10);
1355         tw_newrow();
1356         tw_iglobal("Box lz:", &box_lz, 10);
1357         tw_iglobal("    uz:", &box_uz, 10);
1358         tw_newrow();
1359         tw_newrow();
1360         tw_iglobal("Box2 lx:", &box2_lx, 10);
1361         tw_iglobal("    ux:", &box2_ux, 10);
1362         tw_newrow();
1363         tw_iglobal("Box2 ly:", &box2_ly, 10);
1364         tw_iglobal("    uy:", &box2_uy, 10);
1365         tw_newrow();
1366         tw_iglobal("Box2 lz:", &box2_lz, 10);
1367         tw_iglobal("    uz:", &box2_uz, 10);
1368         
1369         tw_end_dialog();
1370 }
1371 
1372 
1373 void  octants_tool(int x, int y)
1374 {
1375         static void *tool=NULL;
1376 
1377         if (tool)
1378         {
1379                 tw_show_tool(tool);
1380                 return;
1381         }
1382 
1383         tool = (void *)tw_tool("DODECANTS Tool", x, y);
1384 
1385         p_dir = (void *) tw_sglobal("Directory:", directory_name, 32);
1386         tw_newrow();
1387         p_data_file = (void *)tw_sglobal("Data list:", data_list_name, 32);
1388         tw_newrow();
1389         p_thresh_file = (void *)tw_sglobal("Initial model:", mean_model_name, 32);
1390         tw_newrow();
1391         p_tvs_file = (void *)tw_sglobal(".tvs file:", tvs_file_name, 32);
1392         tw_newrow();
1393         p_output_file = (void*) tw_sglobal("Output files:", output_file_name, 32);
1394         tw_newrow();
1395 
1396         tw_label("1. MI Coregistration:");
1397         tw_newrow();
1398         tw_button("Bulk coreg", bulk_coreg_proc, NULL);
1399         tw_button("CC start", coreg_checker_proc, NULL);
1400         tw_button("CC next", coreg_checker_next, NULL);
1401         tw_button("CC end", coreg_checker_end, NULL);
1402         tw_newrow();
1403 
1404         tw_label("2. EM Segmentation");
1405         tw_newrow();
1406         tw_button("Params", segmentation_dialog_proc, NULL);
1407         tw_button("Display box", display_boxmask_proc, NULL);
1408         tw_button("Image Prep", image_prep_proc, NULL);
1409         tw_button("Seq->one", seq_one_proc, NULL);
1410         tw_button("One->seq", one_seq_proc, NULL);
1411         tw_newrow();
1412         tw_button("EM: model", bulk_em_seg_proc, NULL);
1413         tw_button("EM: seg", bulk_em_seg2_proc, NULL);
1414         tw_button("Model output", output_model_parameters_proc, NULL);
1415         tw_newrow();
1416         tw_button("Bulk CSF count", bulk_classifier_proc, NULL);
1417         tw_button("TVS->V12", tvs_to_v12_proc, NULL);
1418         tw_newrow();
1419 
1420         tw_label("4. KNN Classifier:");
1421         tw_newrow();
1422         tw_choice("Fix mean age", fixmeanage_choice_proc, 0, "Off", "On", NULL);
1423         (void *)tw_fglobal("Mean age:  ", &global_mean_age, 5);
1424         tw_newrow();
1425         tw_choice("Optimise scales", scale_opt_choice_proc, 0, "Off", "On", NULL);
1426         tw_newrow();
1427         (void *)tw_iglobal("KNEAR: ", &knear, 5);
1428         tw_button("Run classifier", octants_results_classifier, NULL);
1429         tw_newrow();
1430 
1431         tw_end_tool();
1432 }
1433 

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