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

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedCoreg_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 "tlmedCoreg_tool.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #   include <config.h>
 50 #endif
 51 
 52 #include <math.h>
 53 #include <stdio.h>
 54 #include <tina/sys/sysDef.h>
 55 #include <tina/sys/sysPro.h>
 56 #include <tina/image/imgDef.h>
 57 #include <tina/image/imgPro.h>
 58 #include <tina/geometry/geomDef.h>
 59 #include <tina/math/mathDef.h>
 60 #include <tina/math/mathPro.h>
 61 #include <tina/file/fileDef.h>
 62 #include <tina/file/filePro.h>
 63 #include <tina/medical/medPro.h>
 64 #include <tinatool/draw/drawDef.h>
 65 #include <tinatool/draw/drawPro.h>
 66 #include <tinatool/wdgts/wdgtsPro.h>
 67 #include <tinatool/tlbase/tlbasePro.h>
 68 #include <tinatool/tlmedical/tlmedCoreg_view.h>
 69 #include <tinatool/tlmedical/tlmedCoreg_auto.h>
 70 #include <tinatool/tlmedical/tlmedCoreg_noise.h>
 71 
 72 
 73 static double xcentre=128.0, ycentre=128.0, zcentre=40.0;
 74 static double bscale = 0.75, border = 2.0;
 75 static int fnum=0;
 76 static struct air16 air1;
 77 static void *ape1=NULL, *ape2=NULL, *ape3=NULL;
 78 static void *ape4=NULL, *ape5=NULL, *ape6=NULL;
 79 static void *ape7=NULL, *ape8=NULL, *ape9=NULL;
 80 static void *ate1=NULL, *ate2=NULL, *ate3=NULL;
 81 static void *af1=NULL, *af2=NULL;
 82 static void *avs1=NULL, *avs2=NULL, *avs3=NULL;
 83 static void *avr1=NULL, *avr2=NULL, *avr3=NULL;
 84 static void *ads1=NULL,*ads2=NULL,*ads3=NULL;
 85 static void *adr1=NULL,*adr2=NULL,*adr3=NULL;
 86 static void *pfnum=NULL;
 87 static void *p_thresh=NULL, *p_border=NULL,*p_scale=NULL;
 88 static void *reslice_choice=NULL;
 89 static void *latch_choice=NULL;
 90 static void *p_air_file=NULL;
 91 
 92 extern double xscale;
 93 extern double yscale;
 94 extern double zscale;
 95 extern double nxscale;
 96 extern double nyscale;
 97 extern double nzscale;
 98 extern double ex[3], ey[3], ez[3], ec[3];
 99 static double zoomfac =1.0;
100 static void *psx,*psy,*psz;
101 static char datafile[128];
102 static char output_path[128];
103 static char cpfile[128];
104 
105 static int coreg_im;
106 static int mask = 3;
107 
108 /* Coreg Interface flags */
109 /* mi_switch : 0 = traditional mutual information
110                1 = unbiased mutual information
111    mod_switch : for edge-based coregistration
112                 0 = edges have same sign in both image volumes
113                 1 = edges in the two image volumes have opposite signs
114    a_switch : representation of angles in the transformation model
115               0 = degenerate quaternions, with bulk of rotation stored
116               1 = non-degenerate quaternions: full rotation used
117               2 = Euler angles (for MI)
118    bin_switch : 0 = no of bins in MI joint histogram fully automatic
119                 1 = no of bins in MI joint histogram limited
120    smooth_switch: 0 = gaussian smoothing of joint histogram
121                   1 = iterated tangential smoothing
122    smooth_power: s.d. of Gaussian, iterations of tangential
123    peak_switch : main peak remover
124    coreg_noise : noise to be added to imptrs
125 */
126 
127 int mi_switch=0;
128 int mod_switch=0;
129 int a_switch=0;
130 int bin_switch=0;
131 int bin_size_switch=0;
132 int maxig_switch=0;
133 double coreg_threshold=600;
134 int smooth_switch=0;
135 double smooth_power=4.0;
136 int peak_switch=1;
137 double coreg_noise=0.0;
138 double coreg_target_noise=0.0;
139 int step_switch = 0;
140 
141 
142 /********************* Air file handling functions ********************/
143 
144 
145 static void air_init()
146 {
147    air1.r.x_dim = air1.s.x_dim = 256;
148    air1.r.y_dim = air1.s.y_dim = 256;
149    air1.r.z_dim = air1.s.z_dim = 40;
150 
151    air1.r.x_size = air1.s.x_size = 1.0;
152    air1.r.y_size = air1.s.y_size = 1.0;
153    air1.r.z_size = air1.s.z_size = 1.0;
154 
155    air1.e[0][0] = 1.0;
156    air1.e[0][1] = 0.0;
157    air1.e[0][2] = 0.0;
158    air1.e[1][0] = 0.0;
159    air1.e[1][1] = 1.0;
160    air1.e[1][2] = 0.0;
161    air1.e[2][0] = 0.0;
162    air1.e[2][1] = 0.0;
163    air1.e[2][2] = 1.0;
164    air1.e[3][0] = 0.0;
165    air1.e[3][1] = 0.0;
166    air1.e[3][2] = 0.0;
167 }
168 
169 
170 void air_dump_proc(struct air16 *air1)
171 {
172     double yaw, pitch, roll, new_e[4][4];
173     double min_size = MIN(MIN(air1->s.x_size,air1->s.y_size),air1->s.z_size);
174 
175     new_e[0][0] = air1->e[0][0]*air1->s.x_size/min_size;
176     new_e[1][0] = air1->e[1][0]*air1->s.x_size/min_size;
177     new_e[2][0] = air1->e[2][0]*air1->s.x_size/min_size;
178     new_e[3][0] = air1->e[3][0]*air1->s.x_size/min_size;
179     new_e[0][1] = air1->e[0][1]*air1->s.y_size/min_size;
180     new_e[1][1] = air1->e[1][1]*air1->s.y_size/min_size;
181     new_e[2][1] = air1->e[2][1]*air1->s.y_size/min_size;
182     new_e[3][1] = air1->e[3][1]*air1->s.y_size/min_size;
183     new_e[0][2] = air1->e[0][2]*air1->s.z_size/min_size;
184     new_e[1][2] = air1->e[1][2]*air1->s.z_size/min_size;
185     new_e[2][2] = air1->e[2][2]*air1->s.z_size/min_size;
186     new_e[3][2] = air1->e[3][2]*air1->s.z_size/min_size;
187 
188     format("Standard %s \n",&air1->s_file);
189     format("Reslice %s \n",&air1->r_file);
190 
191     pitch = asin(new_e[1][0]);
192     pitch *= 180.0/3.141592;
193     roll = asin(new_e[2][1]);
194     roll *= 180.0/3.141592;
195     yaw =  asin(new_e[0][2]);
196     yaw *= 180.0/3.141592;
197 
198     format("pitch %f\n",pitch);
199     format("roll %f\n",roll);
200     format("yaw %f\n",yaw);
201     format("trans x %f\n",new_e[3][0]);
202     format("trans y %f\n",new_e[3][1]);
203     format("trans z %f\n",new_e[3][2]);
204 
205 }
206 
207 
208 void reset_air()
209 {
210     tw_fglobal_reset(ape1);
211     tw_fglobal_reset(ape2);
212     tw_fglobal_reset(ape3);
213     tw_fglobal_reset(ape4);
214     tw_fglobal_reset(ape5);
215     tw_fglobal_reset(ape6);
216     tw_fglobal_reset(ape7);
217     tw_fglobal_reset(ape8);
218     tw_fglobal_reset(ape9);
219     tw_fglobal_reset(ate1);
220     tw_fglobal_reset(ate2);
221     tw_fglobal_reset(ate3);
222     tw_sglobal_reset(af1);
223     tw_sglobal_reset(af2);
224     tw_fglobal_reset(avs1);
225     tw_fglobal_reset(avs2);
226     tw_fglobal_reset(avs3);
227     tw_fglobal_reset(avr1);
228     tw_fglobal_reset(avr2);
229     tw_fglobal_reset(avr3);
230     tw_iglobal_reset(ads1);
231     tw_iglobal_reset(ads2);
232     tw_iglobal_reset(ads3);
233     tw_iglobal_reset(adr1);
234     tw_iglobal_reset(adr2);
235     tw_iglobal_reset(adr3);
236 }
237 
238 
239 static void read_air_proc()
240 {
241   strcpy(cpfile,datafile);
242   parse_fname(NULL,cpfile,fnum);
243   read_air16(cpfile,&air1);
244   set_coreg_centre(xcentre,ycentre,zcentre);
245   set_coreg_trans(air1);
246   reset_air();
247 }
248 
249 
250 static void write_air_proc()
251 {
252   get_coreg_trans(&air1);
253   parse_fname(NULL,datafile,fnum);
254   write_air16(datafile,0,&air1);
255 }
256 
257 
258 void coreg_air_invert_proc()
259 {
260     char swapf[128];
261     struct key_info key;
262     double new_e[4][4];
263     double min_size;
264 
265     get_coreg_trans(&air1);
266     strcpy(swapf,air1.r_file);
267     strcpy(air1.r_file,air1.s_file);
268     strcpy(air1.s_file,swapf);
269     key = air1.r;
270     air1.r = air1.s;
271     air1.s = key;
272 /* rescale back to cubic coordinates */
273     min_size = MIN(MIN(air1.r.x_size,air1.r.y_size),air1.r.z_size);
274     new_e[0][0] = air1.e[0][0]*air1.r.x_size/min_size;
275     new_e[0][1] = air1.e[1][0]*air1.r.x_size/min_size;
276     new_e[0][2] = air1.e[2][0]*air1.r.x_size/min_size; 
277     new_e[0][3] = air1.e[3][0]*air1.r.x_size/min_size;
278     new_e[1][0] = air1.e[0][1]*air1.r.y_size/min_size; 
279     new_e[1][1] = air1.e[1][1]*air1.r.y_size/min_size;
280     new_e[1][2] = air1.e[2][1]*air1.r.y_size/min_size; 
281     new_e[1][3] = air1.e[3][1]*air1.r.y_size/min_size;
282     new_e[2][0] = air1.e[0][2]*air1.r.z_size/min_size;
283     new_e[2][1] = air1.e[1][2]*air1.r.z_size/min_size;
284     new_e[2][2] = air1.e[2][2]*air1.r.z_size/min_size;
285     new_e[2][3] = air1.e[3][2]*air1.r.z_size/min_size;
286 
287 /*  translation of coordinates */
288     new_e[3][0] = -(new_e[0][3]*new_e[0][0]
289                   + new_e[1][3]*new_e[1][0]
290                   + new_e[2][3]*new_e[2][0]);
291     new_e[3][1] = -(new_e[0][3]*new_e[0][1]
292                   + new_e[1][3]*new_e[1][1]
293                   + new_e[2][3]*new_e[2][1]);
294     new_e[3][2] = -(new_e[0][3]*new_e[0][2]
295                   + new_e[1][3]*new_e[1][2]
296                   + new_e[2][3]*new_e[2][2]);
297 
298     min_size = MIN(MIN(air1.s.x_size,air1.s.y_size),air1.s.z_size);
299 /* back to non cubic pixel coordinates */
300     air1.e[0][0] = new_e[0][0]*min_size/air1.s.x_size;
301     air1.e[1][0] = new_e[1][0]*min_size/air1.s.x_size;
302     air1.e[2][0] = new_e[2][0]*min_size/air1.s.x_size;
303     air1.e[3][0] = new_e[3][0]*min_size/air1.s.x_size;
304     air1.e[0][1] = new_e[0][1]*min_size/air1.s.y_size;
305     air1.e[1][1] = new_e[1][1]*min_size/air1.s.y_size;
306     air1.e[2][1] = new_e[2][1]*min_size/air1.s.y_size;
307     air1.e[3][1] = new_e[3][1]*min_size/air1.s.y_size;
308     air1.e[0][2] = new_e[0][2]*min_size/air1.s.z_size;
309     air1.e[1][2] = new_e[1][2]*min_size/air1.s.z_size;
310     air1.e[2][2] = new_e[2][2]*min_size/air1.s.z_size;
311     air1.e[3][2] = new_e[3][2]*min_size/air1.s.z_size;
312     set_coreg_trans(air1);
313     reset_air();
314 }
315 
316 
317 /*************** Coords and slice handling functions *******************/
318 
319 
320 static void coreg_seq_init_funcs(void)
321 {
322     /* Wrapper around coreg_slice_init and coreg_coords_init allowing  */ 
323     /* them to be called from a void function                          */
324 
325     Sequence *seq = seq_get_current();
326     
327     coreg_slice_init(seq, seq_dim_to_vec3());
328     coreg_coords_init();
329 }
330 
331 
332 static void coreg_seq_init_funcs_init(void)
333 {
334     /* Add seq_init_funcs to seq tools init funcs */
335 
336     seq_init_funcs_add(coreg_seq_init_funcs);
337 }
338 
339 
340 void reset_centre(double new_x, double new_y, double new_z)
341 {
342     xcentre = new_x;
343     ycentre = new_y;
344     zcentre = new_z;
345     tw_fglobal_reset(psx);
346     tw_fglobal_reset(psy);
347     tw_fglobal_reset(psz);
348 }
349 
350 
351 static void inc_fnum_proc(void)
352 {
353   fnum++;
354   tw_iglobal_reset(pfnum);
355 }
356 
357 
358 static void dec_fnum_proc(void)
359 {
360   fnum--;
361   tw_iglobal_reset(pfnum);
362 }
363 
364 
365 void coreg_zoom_proc()
366 {
367    set_coreg_centre(xcentre,ycentre,zcentre);
368    set_coreg_zoom(zoomfac,mask);
369    reset_tv_coords();
370    coreg_redraw();
371 }
372 
373 
374 static void coreg_push_proc()
375 {
376    Imregion *roi;
377    Imrect *im = NULL;
378 
379    if (coreg_im == 0)
380    {
381        roi = tv_get_im_roi(zcoreg_tv_get());
382        im = seq_slicez((float)zcentre,roi, coreg_bproj);
383        rfree(roi);
384    }
385    if (coreg_im == 1)
386    {
387        roi = tv_get_im_roi(ycoreg_tv_get());
388        im = seq_slicey((float)ycentre,roi, coreg_bproj);
389        rfree(roi);
390    }
391    if (coreg_im == 2)
392    {
393        roi = tv_get_im_roi(xcoreg_tv_get());
394        im = seq_slicex((float)xcentre,roi,coreg_bproj);
395        rfree(roi);
396    }
397    if (im!=NULL) stack_push((void *)im, IMRECT,im_free);
398    imcalc_draw(imcalc_tv_get());
399    imcalc_draw(imcal2_tv_get());
400    image_choice_reset();
401 }
402 
403 
404 void coreg_check_proc(int val)
405 {
406 /* Call to reset the choice immediately allows this function to be called from the octants tool */
407 
408     static int old_val=0;
409     Sequence *seq = seq_get_current();
410 
411     if (val==0) latch_clear();
412     if(val>0)
413     {
414        if (old_val <1) strcpy(air1.s_file,seq->filename);
415        latch_slices(val);
416     }
417     old_val = val;
418     tw_choice_reset(latch_choice, val);
419 }
420 
421 
422 static void last_slice_proc()
423 {
424     if(coreg_im == 0) reset_centre(xcentre,ycentre,zcentre-1.0);
425     if(coreg_im == 1) reset_centre(xcentre,ycentre-1.0,zcentre);
426     if(coreg_im == 2) reset_centre(xcentre-1.0,ycentre,zcentre);
427 }
428 
429 
430 static void next_slice_proc()
431 {
432     if(coreg_im == 0) reset_centre(xcentre,ycentre,zcentre+1.0);
433     if(coreg_im == 1) reset_centre(xcentre,ycentre+1.0,zcentre);
434     if(coreg_im == 2) reset_centre(xcentre+1.0,ycentre,zcentre);
435 }
436 
437 
438 void coreg_batch_reslice()
439 {
440         /* Reslice whole volume at once along the z axis. PAB 12/10/01 */
441 
442         int i, old_coreg_im, uy, ly;
443         Imrect *y_im;
444         Imregion *roi;
445         float old_zslice;
446 
447         y_im = get_ycoreg_im();
448         if(y_im==NULL) return;
449         roi = y_im->region;
450 
451         ly = roi->ly;
452         uy = roi->uy;
453 
454         old_zslice=zcentre;
455         old_coreg_im=coreg_im;
456         coreg_im=0;
457 
458         for(i=ly; i<uy; i++)
459         {
460                    zcentre = (float)i +0.5;
461                    reset_centre(xcentre, ycentre, zcentre);
462                    coreg_push_proc();
463         }
464 
465         zcentre = old_zslice;
466         reset_centre(xcentre, ycentre, zcentre);
467         coreg_im=old_coreg_im;
468 }
469 
470 
471 static void batch_seq_proc()
472 {
473         Vec3 *iscale=NULL;
474 
475         /* Delete and replace sequence with results of batch reslicing */
476 
477         iscale = vec3_alloc();
478         iscale->el[0] = nxscale;
479         iscale->el[1] = nyscale;
480         iscale->el[2] = nzscale;
481         stack_seq(iscale);
482 }
483 
484 
485 static void set_scales()
486 {
487     seq_scales_update(nxscale, nyscale, nzscale);
488 }
489 
490 
491 static void set_rfile_proc(struct air16 *rair)
492 {
493     set_seq_rfile(rair,fnum);
494 }
495 
496 
497 /************ Coregistration and coil correction functions ***************/
498 
499 
500 static void batch_norm_proc()
501 {
502    /* coil correct whole volume NAT 01/05/2002
503     *   Add mod for tool / lib separation - GAB 18 Feb 2003
504     *   Altered to receive image list and loop through images pushing to stack
505     *   and displaying.
506     */
507 
508    List   *imlist, *imptr;
509 
510    imlist = seq_norm(seq_get_current(), 15.0, 100.0);
511 
512    for (imptr = imlist; imptr != NULL; imptr = imptr->next)
513    {
514       Imrect *im = (Imrect *) imptr->to;
515 
516       stack_push((void*) im,IMRECT,im_free);
517       imcalc_draw(imcalc_tv_get());
518    }
519 }
520 
521 
522 static void coreg_noise_proc(float noise_multiplier)
523 {
524    if(noise_multiplier!=0) coreg_pab_noise(noise_multiplier);
525 }
526 
527 
528 static void coreg_target_noise_proc(float noise_multiplier)
529 {
530    if(noise_multiplier!=0) coreg_add_noise(noise_multiplier);
531 }
532 
533 
534 static void coreg_noise_caller_proc()
535 {
536    coreg_noise_proc(coreg_noise);
537 }
538 
539 
540 static void coreg_target_noise_caller_proc()
541 {
542    coreg_target_noise_proc(coreg_target_noise);
543 }
544 
545 
546 void coreg_auto_proc()
547 {
548    coreg_auto(coreg_threshold, bscale, border, mask);
549    reset_air();
550    reset_tv_coords();
551    coreg_redraw();
552 }
553 
554 
555 void mui_coreg_auto_proc()
556 {
557    mui_coreg_auto(coreg_threshold, bscale, border, mask);
558    reset_air();
559    reset_tv_coords();
560    coreg_redraw();
561    imptrs_reset();
562 }
563 
564 
565 static void coreg_covar_proc()
566 {
567         if(mi_switch!=1)
568         {
569                 error("Switch to unbiased MI first!\n", warning);
570                 return;
571         }
572         if(step_switch!=2&&a_switch!=2)
573         {
574                 error("Must use Euler angles with Thacker or Brainweb steps.\n", warning);
575         }
576 
577         coreg_covar(coreg_threshold,border,mask);
578 }
579 
580 
581 /************** Bulk processing functions *****************************/
582 /* Description of the functions:                                      */
583 /*                                                                    */
584 /* bulk_coreg_proc performs 1000 coregistrations of the latched       */
585 /* data, always starting from the current air file, adding a new      */
586 /* noise field at the s.d. given in the interface each time, and      */
587 /* outputs the coreg_vec to the filname specified in the interface.   */
588 /*                                                                    */
589 /* process_bulk_coreg_results takes the output from bulk_coreg_proc   */
590 /* and calculates a covariance matrix, outputting two files: the      */
591 /* matrix and the average of the coreg_vecs.                          */
592 /*                                                                    */
593 /* bulk_covar_proc calls bulk_covar in tlmedCoreg_auto.c to scan 100  */
594 /* points around the minimum of the coreg cost function, calculating  */
595 /* the covariance at each point adn writing all 100 matrices to a     */
596 /* file.                                                              */
597 /*                                                                    */
598 /**********************************************************************/
599 
600 
601 void process_bulk_coreg_results(void)
602 {
603         /* Take the results from the above function and process to produce a covariance matrix */
604 
605         FILE *input_file=NULL, *output_file=NULL, *output_file_av=NULL;
606         char output_name[128], output_name_av[128];
607         double **data_matrix=NULL, *means=NULL, **cov_matrix=NULL;
608         int i, j, k, length=9, rows=0;
609 
610         (void) strip_spaces(output_path);
611         if((input_file = fopen(output_path, "r"))==NULL)
612         {
613                 format("Could not open input file\n");
614                 return;
615         }
616         (void) string_append(output_name, output_path, ".matrix", NULL);
617         (void) string_append(output_name_av, output_path, ".av", NULL);
618         if((output_file = fopen(output_name, "w"))==NULL)
619         {
620                 format("Could not open output file\n");
621                 return;
622         }
623         if((output_file_av = fopen(output_name_av, "w"))==NULL)
624         {
625                 format("Could not open output av file\n");
626                 return;
627         }
628         data_matrix = darray_alloc(0, 0, 1000, length);
629         means = dvector_alloc(0, length);
630         cov_matrix = darray_alloc(0, 0, length,length);
631 
632         for(i=0; i<length; i++)
633         {
634                 means[i] = 0.0;
635                 for(j=0; j<length; j++)
636                 {
637                         cov_matrix[i][j] = 0.0;
638                 }
639         }
640 
641         for(i=0; i<1000; i++)
642         {
643                 for(j=0; j<length; j++)
644                 {
645                         fscanf(input_file, "%lf", &data_matrix[i][j]);
646                         if(feof(input_file)) break;
647                         data_matrix[i][j] = data_matrix[i][j]/100.0;
648                         means[j] += data_matrix[i][j];
649                 }
650                 if(feof(input_file)) break;
651                 rows++;
652         }
653 
654         for(i=0; i<length; i++) means[i] = means[i]/rows;
655         for(i=0; i<rows; i++)
656         {
657                 for(j=0; j<length; j++)
658                 {
659                         data_matrix[i][j] -= means[j];
660                 }
661         }
662 
663         for(k=0; k<rows; k++)
664         {
665                 for(i=0; i<length; i++)
666                 {
667                         for(j=0; j<length; j++)
668                         {
669                                 cov_matrix[i][j] += data_matrix[k][i] *
670                                         data_matrix[k][j];
671                         }
672                 }
673         }
674         for(i=0; i<length; i++)
675         {
676                 fprintf(output_file_av, "%f ", means[i]);
677                 for(j=0; j<length; j++)
678                 {
679                         cov_matrix[i][j] = cov_matrix[i][j] / (rows-1);
680                         fprintf(output_file, "%f ", cov_matrix[i][j]*1000000);
681                 }
682                 fprintf(output_file, "\n");
683         }
684 
685         darray_free(data_matrix, 0, 0, 1000, length);
686         dvector_free(means, 0);
687         darray_free(cov_matrix, 0, 0, length, length);
688         fclose(input_file);
689         fclose(output_file);
690         fclose(output_file_av);
691 }
692 
693 void bulk_coreg_proc(void)
694 {
695         /* Repeatedly add noise to the latched images, repeat the
696         coregistration, and output the coreg vector. Assume that the
697         images have already been loaded and latched */
698 
699         FILE *output_file=NULL;
700         int i, j, npar;
701         double *vec=NULL;
702 
703         vec = dvector_alloc(0, 10);
704 
705         (void) strip_spaces(output_path);
706         if((output_file = fopen(output_path, "w"))==NULL)
707         {
708                 format("Could not open output file\n");
709                 return;
710         }
711 
712         for(i=0; i<1000; i++)
713         {
714                 read_air_proc();
715                 coreg_target_noise_proc(1.0);
716                 coreg_noise_proc(coreg_noise);
717                 mui_coreg_auto_proc();
718                 stack_clear();
719                 npar = coreg_get_vec(vec, mask);
720 
721                 for(j=0; j<npar; j++)
722                 {
723                         fprintf(output_file, "%f ", vec[j]);
724                 }
725                 fprintf(output_file, "\n");
726                 imptrs_reset();
727         }
728         fclose(output_file);
729         dvector_free(vec, 0);
730         process_bulk_coreg_results();
731 }
732 
733 
734 void process_bulk_covar_proc(void)
735 {
736         FILE *input_file=NULL, *output_file=NULL;
737         char output_path_2[128];
738         double **mat, *step_vec=NULL, **master_array=NULL;
739         int i, j, k, no_records=100, length=9;
740 
741         (void) strip_spaces(output_path);
742         if((input_file = fopen(output_path, "r"))==NULL)
743         {
744                 format("Could not open input file\n");
745                 return;
746         }
747         (void) string_append(output_path_2, output_path, ".stddev", NULL);
748         if((output_file = fopen(output_path_2, "w"))==NULL)
749         {
750                 format("Could not open output file\n");
751                 return;
752         }
753         fscanf(input_file,"%d \n", &no_records);
754 
755         mat = darray_alloc(0, 0, length, length);
756         master_array = darray_alloc(0, 0, no_records, (length*2));
757 
758         for(i=0; i<no_records; i++)
759         {
760                 step_vec = dvector_alloc(0, length);
761                 for(j=0; j<length; j++)
762                 {
763                         fscanf(input_file, "%lf", &step_vec[j]);
764                 }
765                 for(j=0; j<length; j++)
766                 {
767                         for(k=0; k<length; k++)
768                         {
769                                 fscanf(input_file, "%lf", &mat[j][k]);
770                         }
771                 }
772                 for(j=0; j<length; j++)
773                 {
774                         master_array[i][j] = step_vec[j];
775                 }
776                 for(j=0; j<length; j++)
777                 {
778                         master_array[i][j+length] = sqrt(mat[j][j]);
779                 }
780 
781                 dvector_free(step_vec, 0);
782         }
783 
784         for(i=0; i<no_records; i++)
785         {
786                 for(j=0; j<(length*2); j++)
787                 {
788                         fprintf(output_file, "%f ", master_array[i][j]);
789                 }
790                 fprintf(output_file, "\n");
791         }
792 
793         darray_free(master_array, 0, 0, no_records, (length*2));
794         darray_free(mat, 0, 0, length, length);
795 
796         fclose(output_file);
797         fclose(input_file);
798 }
799 
800 
801 
802 static void bulk_covar_proc()
803 {
804         if(mi_switch!=1)
805         {
806                 error("Switch to unbiased MI first!\n", warning);
807                 return;
808         }
809         if(step_switch!=2&&a_switch!=2)
810         {
811                 error("Must use Euler angles with Thacker or Brainweb steps.\n", warning);
812         }
813 
814         (void) strip_spaces(output_path);
815         coreg_covar_bulk(coreg_threshold, border, mask, output_path);
816 }
817 
818 
819 /**************** External communication for octants tool ********************/
820 
821 
822 void coreg_reset_threshbordscale(double ext_threshold, double ext_border, double ext_scale)
823 {
824         coreg_threshold = ext_threshold;
825         border = ext_border;
826         bscale = ext_scale;
827         tw_fglobal_reset(p_thresh);
828         tw_fglobal_reset(p_border);
829         tw_fglobal_reset(p_scale);
830 }
831 
832 
833 void coreg_air_load(char *filename)
834 {
835         strcpy(datafile, filename);
836         tw_sglobal_reset(p_air_file);
837         read_air_proc();
838 }
839 
840 
841 /**************** Dialogues, choices etc... ********************************/
842 
843 
844 static void     tv_choice_proc(int val)
845 {
846     coreg_im = val;
847 
848     switch (val)
849     {
850         case 0:
851           if (!zcoreg_tv_get()) zcoreg_tv_make();
852           tv_set_next(zcoreg_tv_get());
853         break;
854         case 1:
855           if (!ycoreg_tv_get()) ycoreg_tv_make();
856           tv_set_next(ycoreg_tv_get());
857         break;
858         case 2:
859           if (!xcoreg_tv_get()) xcoreg_tv_make();
860           tv_set_next(xcoreg_tv_get());
861         break;
862     }
863 }
864 
865 
866 void coreg_reslice_choice_proc(int val)
867 {
868 /* Call to reset the choice immediately allows this function to be called from the octants tool */
869 
870    seq_interp_choice(val);
871    tw_choice_reset(reslice_choice, val);
872 }
873 
874 
875 void coreg_par_choice_notify_proc(int value)
876 {
877    mask = value;
878    set_coreg_zoom(zoomfac,mask);
879 }
880 
881 
882 static void mod_choice_proc(int val)
883 {
884    mod_switch = val;
885 }
886 
887 
888 static void mi_choice_proc(int val)
889 {
890         mi_switch=val;
891 }
892 
893 
894 static void angle_rep_choice_proc(int val)
895 {
896         a_switch = val;
897 }
898 
899 
900 static void bin_limiters_choice_proc(int val)
901 {
902         bin_switch = val;
903 }
904 
905 
906 static void bin_size_choice_proc(int val)
907 {
908         bin_size_switch=val;
909 }
910 
911 static void maxig_switch_proc(int val)
912 {
913         maxig_switch = val;
914 }
915 
916 
917 static void smooth_switch_proc(int val)
918 {
919         smooth_switch = val;
920 }
921 
922 
923 static void peak_switch_proc(int val)
924 {
925         peak_switch = val;
926 }
927 
928 
929 static void step_defs_proc(int val)
930 {
931         step_switch = val;
932 }
933 
934 
935 static void air_param_dialog_proc()
936 {
937     static void *dialog = NULL;
938 
939     if (dialog)
940     {
941         tw_show_dialog(dialog);
942         return;
943     }
944     dialog = (void *)tw_dialog("AIR Parameters");
945 
946     af1 = (void *) tw_sglobal("Standard:",air1.s_file,42);
947     tw_newrow();
948     ads1 = (void *) tw_iglobal("dimensions :",&air1.s.x_dim,10);
949     ads2 = (void *) tw_iglobal("",&air1.s.y_dim,10);
950     ads3 = (void *) tw_iglobal("",&air1.s.z_dim,10);
951     tw_button("set S", set_seq_sfile, &air1);
952     tw_newrow();
953     avs1 = (void *) tw_fglobal("Scale:", &nxscale, 12);
954     avs2 = (void *) tw_fglobal("", &nyscale, 12);
955     avs3 = (void *) tw_fglobal("", &nzscale, 12);
956     tw_newrow();
957     af2 = (void *) tw_sglobal("Reslice :",air1.r_file,42);
958     tw_newrow();
959     adr1 = (void *) tw_iglobal("dimensions :",&air1.r.x_dim,10);
960     adr2 = (void *) tw_iglobal("",&air1.r.y_dim,10);
961     adr3 = (void *) tw_iglobal("",&air1.r.z_dim,10);
962     tw_button("set R",set_rfile_proc,&air1);
963     tw_newrow();
964     pfnum = (void *) tw_iglobal("data",&fnum,3);
965     tw_button("dec", dec_fnum_proc, NULL);
966     tw_button("inc", inc_fnum_proc, NULL);
967     tw_newrow();
968     avr1 = (void *) tw_fglobal("Scale:", &xscale, 12);
969     avr2 = (void *) tw_fglobal("", &yscale, 12);
970     avr3 = (void *) tw_fglobal("", &zscale, 12);
971     tw_newrow();
972     ape1 = (void *) tw_fglobal("Er :", &ex[0], 16);
973     ape2 = (void *) tw_fglobal("", &ex[1], 16);
974     ape3 = (void *) tw_fglobal("", &ex[2], 16);
975     tw_newrow();
976     ape4 = (void *) tw_fglobal("   :", &ey[0], 16);
977     ape5 = (void *) tw_fglobal("", &ey[1], 16);
978     ape6 = (void *) tw_fglobal("", &ey[2], 16);
979     tw_newrow();
980     ape7 = (void *) tw_fglobal("   :", &ez[0], 16);
981     ape8 = (void *) tw_fglobal("", &ez[1], 16);
982     ape9 = (void *) tw_fglobal("", &ez[2], 16);
983     tw_newrow();
984     ate1 = (void *) tw_fglobal("Et :", &ec[0], 16);
985     ate2 = (void *) tw_fglobal("", &ec[1], 16);
986     ate3 = (void *) tw_fglobal("", &ec[2], 16);
987     tw_newrow();
988     p_thresh = (void *) tw_fglobal("threshold:",&coreg_threshold,8);
989     p_border = (void *) tw_fglobal("blur scale:",&bscale,8);
990     p_scale = (void *) tw_fglobal("border:",&border,8);
991 
992     tw_end_dialog();
993 }
994 
995 
996 static void advanced_dialog_proc()
997 {
998     static void *advanced_dialog = NULL;
999 
1000     if (advanced_dialog)
1001     {
1002         tw_show_dialog(advanced_dialog);
1003         return;
1004     }
1005     advanced_dialog = (void *)tw_dialog("Advanced Coreg Tool");
1006 
1007    tw_label("Advanced Coreg Functions");
1008    tw_newrow();
1009    tw_choice("Modality : ", mod_choice_proc, 0, "single", "multi", NULL);
1010    tw_choice("MI", mi_choice_proc, 0, "Trad.", "Chi-sqr", NULL);
1011    tw_newrow();
1012    tw_choice("Angles: ", angle_rep_choice_proc,
1013                0, "Quat (Partial)", "Quat (Full)", "Euler", NULL);
1014    tw_newrow();
1015 
1016    tw_choice("Max. Ig.", maxig_switch_proc, 0, "Off", "On", NULL);
1017    tw_choice("Limit bins", bin_limiters_choice_proc, 0, "Off", "On", NULL);
1018    tw_newrow();
1019    tw_choice("Bin size", bin_size_choice_proc, 0, "1 sigma", "0.5 sigma", NULL);
1020    tw_newrow();
1021 
1022    tw_choice("Hist. smooth", smooth_switch_proc, 0, "Off", "Gauss", "Tang", NULL);
1023    tw_fglobal("Power", &smooth_power, 8);
1024    tw_newrow();
1025    tw_choice("Main peak remover", peak_switch_proc, 1, "Off", "On", NULL);
1026    tw_newrow();
1027    tw_fglobal("Source noise", &coreg_noise, 8);
1028    tw_button("Add noise", coreg_noise_caller_proc, NULL);
1029    tw_newrow();
1030    tw_fglobal("Target noise", &coreg_target_noise, 8);
1031    tw_button("Add noise", coreg_target_noise_caller_proc, NULL);
1032    tw_newrow();
1033    tw_button("MI auto",mui_coreg_auto_proc,NULL);
1034    tw_button("MI covar",coreg_covar_proc,NULL);
1035    tw_newrow();
1036    tw_button("Monte-Carlo coreg", bulk_coreg_proc, NULL);
1037    tw_button("Monte-Carlo covar", bulk_covar_proc, NULL);
1038    tw_button("Proccess BN covar", process_bulk_covar_proc, NULL);
1039    tw_newrow();
1040    tw_choice("Step defs:", step_defs_proc, 0, "Brainweb", "Thacker", "Auto", NULL);
1041    tw_newrow();
1042    tw_button("Stack->seq", batch_seq_proc, NULL);
1043    tw_newrow();
1044    tw_sglobal("Output file: ", output_path, 32);
1045 
1046    tw_end_dialog();
1047 }
1048 
1049 
1050 /**************** Coreg Tool Creation and Destruction **********************/
1051 
1052 
1053 void         coreg_tool(int x, int y)
1054 {
1055    static void *tool=NULL;
1056 
1057    if(tool)
1058    {
1059      tw_show_tool(tool);
1060      return;
1061    }
1062 
1063    air_init();
1064    coreg_seq_init_funcs_init();
1065    tool = tw_tool("Coreg Tool", x, y);
1066    tw_help_button("coreg_tool");
1067    tw_button("Advanced", advanced_dialog_proc, NULL);
1068    tw_choice("Tv : ", tv_choice_proc, 0,
1069              "z coreg", "y coreg", "x coreg", NULL);
1070    tw_newrow();
1071 
1072    psx = (void *) tw_fglobal("Centre x :", &xcentre, 9);
1073    psy = (void *) tw_fglobal(" y :", &ycentre, 9);
1074    psz = (void *) tw_fglobal(" z :", &zcentre, 9);
1075    tw_newrow();
1076    tw_button("<", last_slice_proc, NULL);
1077    tw_button(">", next_slice_proc, NULL);
1078    tw_fglobal("factor",&zoomfac,8);
1079    tw_button("zoom", coreg_zoom_proc, NULL);
1080    tw_newrow();
1081 
1082    reslice_choice = tw_choice("Reslice", coreg_reslice_choice_proc, 0,
1083              "nearest", "linear","sinc5","sinc7",
1084              NULL);
1085    tw_button("push", coreg_push_proc, NULL);
1086    tw_newrow();
1087    latch_choice = tw_choice("Latch ", coreg_check_proc, 0,
1088            "image","anaglyph", "edges", "chequer",  NULL);
1089    tw_newrow();
1090    tw_check("Model", coreg_par_choice_notify_proc, mask,
1091             "Trans", "Rot", "Scale", NULL);
1092    tw_newrow();
1093 
1094    tw_button("seq reslice", coreg_batch_reslice, NULL);
1095    tw_button("seq norm", batch_norm_proc, NULL);
1096    tw_button("set scales", set_scales, NULL);
1097    tw_button("Coreg",coreg_auto_proc,NULL);
1098    tw_newrow();
1099 
1100    p_air_file = (void *)tw_sglobal("Air File:", datafile, 36);
1101    tw_newrow();
1102    tw_button("input", read_air_proc, NULL);
1103    tw_button("output", write_air_proc, NULL);
1104    tw_button("swap",coreg_air_invert_proc,NULL);
1105    tw_button("dump", air_dump_proc, &air1);
1106    tw_button("params", air_param_dialog_proc, NULL);
1107    tw_end_tool();
1108 }
1109 

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