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

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedCoreg_view.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 /**
 47 * @file
 48 * @brief Functions for preparing the images on the coreg tool TVs
 49 *
 50 * paul.bromiley@man.ac.uk  30/6/2004
 51 *
 52 * coreg_get_vec and coreg_set_vec:
 53 * These functions now take three arguments: the vector of
 54 * translation, rotation and scaling parameters, a mask that dictates
 55 * which parameters are present (i.e. the model being used in the coregistration)
 56 * and a flag controlling how angles are represented. Flag=0 gives the Tina 4 style
 57 * angular representation, using degenerate quaternions, but storing the bulk of
 58 * the rotation so that only small angles are used in the optimisation.  Flag=1 gives
 59 * a full implementation of quaternions, and works with the full rotation. Flag=2 gives
 60 * a Euler angle representation: this could in theory suffer from gimbal lock, so
 61 * should not in general be used for coregistration, but is useful for covariance
 62 * estimation since it is a non-degenerate representation of the rotation.
 63 *
 64 **/
 65 
 66 #include "tlmedCoreg_view.h"
 67 
 68 #if HAVE_CONFIG_H
 69 #   include <config.h>
 70 #endif
 71 
 72 #include <stdio.h>
 73 #include <math.h>
 74 #include <tina/sys/sysDef.h>
 75 #include <tina/sys/sysPro.h>
 76 #include <tina/math/mathDef.h>
 77 #include <tina/math/mathPro.h>
 78 #include <tina/image/imgDef.h>
 79 #include <tina/image/imgPro.h>
 80 #include <tina/geometry/geomDef.h>
 81 #include <tina/geometry/geomPro.h>
 82 #include <tina/vision/visDef.h>
 83 #include <tina/vision/visPro.h>
 84 #include <tinatool/draw/drawDef.h>
 85 #include <tinatool/draw/drawPro.h>
 86 #include <tinatool/gphx/gphxDef.h>
 87 #include <tinatool/gphx/gphxPro.h>
 88 #include <tinatool/tlbase/tlbasePro.h>
 89 #include <tinatool/tlmedical/tlmedCoreg_tool.h>
 90 
 91 
 92 #define Bit(mask,num)     (mask&(1<<num))
 93 
 94 static Imrect *xcoreg_im=NULL,*xcoreg_disp=NULL;
 95 static Imrect *ycoreg_im=NULL,*ycoreg_disp=NULL;
 96 static Imrect *zcoreg_im=NULL,*zcoreg_disp=NULL;
 97 static Tv* zcoreg_tv = NULL;
 98 static Tv* ycoreg_tv = NULL;
 99 static Tv* xcoreg_tv = NULL;
100 static double xcentre = 128.0;
101 static double ycentre = 128.0;
102 static double zcentre = 40.0;
103 static Imrect *imref=NULL;
104 double xscale = 1.0;
105 double yscale = 1.0;
106 double zscale = 2.0;
107 double nxscale = 1.0;
108 double nyscale = 1.0;
109 double nzscale = 2.0;
110 double ixscale = 1.0;
111 double iyscale = 1.0;
112 double izscale = 2.0;
113 static double min_scale = 1.0;
114 static Bool zoom_on = false;
115 static double zoom_fac = 1.0;
116 static int imptrlx, imptrux,imptrly, imptruy, imptrlz, imptruz;
117 /*Vec3 coreg_ex, coreg_ey, coreg_ez, coreg_ec;
118 */
119 double ex[3],ey[3],ez[3],ec[3];
120 static void ***imptrs=NULL;
121 static int disp_im_type = 0;
122 static Mat3 *store_rot=NULL, *store_inv = NULL;
123 
124 /* Coreg interface flags: see tlmedCoreg_tool for descriptions */
125 
126 extern double coreg_threshold;
127 extern int mi_switch;
128 extern int mod_switch;
129 extern int a_switch;
130 extern int bin_switch;
131 
132 
133 void store_rot_init()
134 {
135 /* Potential memory leaks fixed PAB 13 / 2 / 2003: check to free matrices before reallocating */
136 
137    if(store_rot!=NULL) mat3_free(store_rot);
138    if(store_inv!=NULL) mat3_free(store_inv);
139 
140    store_rot = mat3_alloc();
141    store_inv = mat3_alloc();
142    store_rot->el[0][0] = ex[0];
143    store_rot->el[0][1] = ex[1];
144    store_rot->el[0][2] = ex[2];
145    store_rot->el[1][0] = ey[0];
146    store_rot->el[1][1] = ey[1];
147    store_rot->el[1][2] = ey[2];
148    store_rot->el[2][0] = ez[0];
149    store_rot->el[2][1] = ez[1];
150    store_rot->el[2][2] = ez[2];
151    *store_inv = mat3_inverse(*store_rot);
152 }
153 
154 void store_rot_reset()
155 {
156    mat3_free(store_rot);
157    mat3_free(store_inv);
158 
159    store_rot = NULL;
160    store_inv = NULL;
161 }
162 
163 Imrect *get_zcoreg_im()
164 {
165     return(zcoreg_im);
166 }
167 
168 Imrect *get_ycoreg_im()
169 {
170     return(ycoreg_im);
171 }
172 
173 Imrect *get_xcoreg_im()
174 {
175     return(xcoreg_im);
176 }
177 
178 void set_zcoreg_im(Imrect *im)
179 {
180     zcoreg_im = im;
181 }
182 
183 void set_ycoreg_im(Imrect *im)
184 {
185     ycoreg_im = im;
186 }
187 
188 void set_xcoreg_im(Imrect *im)
189 {
190     xcoreg_im = im;
191 }
192 
193 
194 int coreg_get_vec_store(double *coreg_vec, int mask)
195 {
196     double vec[4];
197     Mat3 mat;
198     int i=0;
199 
200     if (Bit(mask,0))
201     {
202       coreg_vec[i++] = ec[0];
203       coreg_vec[i++] = ec[1];
204       coreg_vec[i++] = ec[2];
205     }
206     if (Bit(mask,1))
207     {
208       mat.el[0][0] = ex[0];
209       mat.el[0][1] = ex[1];
210       mat.el[0][2] = ex[2];
211       mat.el[1][0] = ey[0];
212       mat.el[1][1] = ey[1];
213       mat.el[1][2] = ey[2];
214       mat.el[2][0] = ez[0];
215       mat.el[2][1] = ez[1];
216       mat.el[2][2] = ez[2];
217       mat = mat3_prod(*store_inv,mat);
218 
219       conv_rot_to_quat(&mat,vec);
220       coreg_vec[i++] = vec[1];
221       coreg_vec[i++] = vec[2];
222       coreg_vec[i++] = vec[3];
223     }
224     if (Bit(mask,2))
225     {
226       coreg_vec[i++] = xscale;
227       coreg_vec[i++] = yscale;
228       coreg_vec[i++] = zscale;
229     }
230     return(i);
231 }
232 
233 
234 int coreg_set_vec_store(double *coreg_vec, int mask)
235 {
236     double x[4];
237     Mat3 mat;
238     int i=0;
239 
240     if (Bit(mask,0))
241     {
242       ec[0] = coreg_vec[i++];
243       ec[1] = coreg_vec[i++];
244       ec[2] = coreg_vec[i++];
245     }
246     if (Bit(mask,1))
247     {
248       x[1] = coreg_vec[i++];
249       x[2] = coreg_vec[i++];
250       x[3] = coreg_vec[i++];
251       if ((x[0] = 1.0 - x[1] * x[1] - x[2] * x[2] - x[3] * x[3]) > 0.0)
252           x[0] = sqrt(x[0]);
253       else
254           return (0);
255 
256       conv_quat_to_rot(x,&mat);
257       mat = mat3_prod(*store_rot,mat);
258 
259       ex[0] = mat.el[0][0];
260       ex[1] = mat.el[0][1];
261       ex[2] = mat.el[0][2];
262       ey[0] = mat.el[1][0];
263       ey[1] = mat.el[1][1];
264       ey[2] = mat.el[1][2];
265       ez[0] = mat.el[2][0];
266       ez[1] = mat.el[2][1];
267       ez[2] = mat.el[2][2];
268     }
269     if (Bit(mask,2))
270     {
271       xscale = coreg_vec[i++];
272       yscale = coreg_vec[i++];
273       zscale = coreg_vec[i++];
274     }
275     return(1);
276 }
277 
278 
279 int coreg_get_vec_quat(double *coreg_vec, int mask)
280 {
281     double vec[4];
282     int i=0;
283 
284     if (Bit(mask,0))
285     {
286       coreg_vec[i++] = ec[0];
287       coreg_vec[i++] = ec[1];
288       coreg_vec[i++] = ec[2];
289     }
290     if (Bit(mask,1))
291     {
292       conv_rot_to_quat_pab(ex, ey, ez,vec);
293       coreg_vec[i++] = vec[0];
294       coreg_vec[i++] = vec[1];
295       coreg_vec[i++] = vec[2];
296       coreg_vec[i++] = vec[3];
297     }
298     if (Bit(mask,2))
299     {
300       coreg_vec[i++] = xscale;
301       coreg_vec[i++] = yscale;
302       coreg_vec[i++] = zscale;
303     }
304     return(i);
305 }
306 
307 
308 int coreg_set_vec_quat(double *coreg_vec, int mask)
309 {
310     double x[4];
311     int i=0;
312 
313     if (Bit(mask,0))
314     {
315       ec[0] = coreg_vec[i++];
316       ec[1] = coreg_vec[i++];
317       ec[2] = coreg_vec[i++];
318     }
319     if (Bit(mask,1))
320     {
321       x[0] = coreg_vec[i++];
322       x[1] = coreg_vec[i++];
323       x[2] = coreg_vec[i++];
324       x[3] = coreg_vec[i++];
325       conv_quat_to_rot_pab(x, ex, ey, ez);
326     }
327     if (Bit(mask,2))
328     {
329       xscale = coreg_vec[i++];
330       yscale = coreg_vec[i++];
331       zscale = coreg_vec[i++];
332     }
333     return(1);
334 }
335 
336 
337 int coreg_get_vec_euler(double *coreg_vec, int mask)
338 {
339    /* Next two functions: euler angle representation for coreg angle parameters PAB 13 / 1 / 2004 */
340 
341     double *vec;
342     int i=0;
343 
344     vec = dvector_alloc(0, 3);
345 
346     if (Bit(mask,0))
347     {
348       coreg_vec[i++] = ec[0];
349       coreg_vec[i++] = ec[1];
350       coreg_vec[i++] = ec[2];
351     }
352     if (Bit(mask,1))
353     {
354       conv_rot_to_euler_pab(ex, ey, ez,vec);
355       coreg_vec[i++] = vec[0];
356       coreg_vec[i++] = vec[1];
357       coreg_vec[i++] = vec[2];
358     }
359     if (Bit(mask,2))
360     {
361       coreg_vec[i++] = xscale;
362       coreg_vec[i++] = yscale;
363       coreg_vec[i++] = zscale;
364     }
365 
366     dvector_free(vec, 0);
367     return(i);
368 }
369 
370 int coreg_set_vec_euler(double *coreg_vec, int mask)
371 {
372     double *x;
373     int i=0;
374 
375     x = dvector_alloc(0, 3);
376 
377     if (Bit(mask,0))
378     {
379       ec[0] = coreg_vec[i++];
380       ec[1] = coreg_vec[i++];
381       ec[2] = coreg_vec[i++];
382     }
383     if (Bit(mask,1))
384     {
385       x[0] = coreg_vec[i++];
386       x[1] = coreg_vec[i++];
387       x[2] = coreg_vec[i++];
388       conv_euler_to_rot_pab(x, ex, ey, ez);
389     }
390     if (Bit(mask,2))
391     {
392       xscale = coreg_vec[i++];
393       yscale = coreg_vec[i++];
394       zscale = coreg_vec[i++];
395     }
396 
397     dvector_free(x, 0);
398     return(1);
399 }
400 
401 
402 int coreg_get_vec(double *coreg_vec, int mask)
403 {
404     int return_val=0;
405 
406     if(a_switch==0)
407     {
408         return_val = coreg_get_vec_store(coreg_vec, mask);
409     }
410     if(a_switch==1)
411     {
412         return_val = coreg_get_vec_quat(coreg_vec, mask);
413     }
414     if(a_switch==2)
415     {
416         return_val = coreg_get_vec_euler(coreg_vec, mask);
417     }
418     return return_val;
419 }
420 
421 
422 int coreg_set_vec(double *coreg_vec, int mask)
423 {
424     int return_val=0;
425 
426     if(a_switch==0)
427     {
428         return_val = coreg_set_vec_store(coreg_vec, mask);
429     }
430     if(a_switch==1)
431     {
432         return_val = coreg_set_vec_quat(coreg_vec, mask);
433     }
434     if(a_switch==2)
435     {
436         return_val = coreg_set_vec_euler(coreg_vec, mask);
437     }
438     return return_val;
439 }
440 
441 
442 Vec3 coreg_ex()
443 {
444    return(vec3(ex[0],ex[1],ex[2]));
445 }
446 
447 Vec3 coreg_ey()
448 {
449    return(vec3(ey[0],ey[1],ey[2]));
450 }
451 
452 Vec3 coreg_ez()
453 {
454    return(vec3(ez[0],ez[1],ez[2]));
455 }
456 
457 Vec3 coreg_ec()
458 {
459    return(vec3(ec[0],ec[1],ec[2]));
460 }
461 
462 void coreg_comp(double *a, Vec3 vec)
463 {
464     a[0] = (double)vec.el[0];
465     a[1] = (double)vec.el[1];
466     a[2] = (double)vec.el[2];
467 }
468 
469 void reset_tv_coords()
470 {
471     if (zcoreg_tv!=NULL)
472     {
473         zcoreg_tv->centre3 = coreg_ec();
474         zcoreg_tv->centre2.el[0] = xcentre;
475         zcoreg_tv->centre2.el[1] = ycentre;
476         zcoreg_tv->ex3 = coreg_ex();
477         zcoreg_tv->ey3 = coreg_ey();
478         zcoreg_tv->ez3 = coreg_ez();
479         zcoreg_tv->scalex = xscale/min_scale;
480         zcoreg_tv->scaley = yscale/min_scale;
481     }
482 
483     if (ycoreg_tv!=NULL)
484     {
485         ycoreg_tv->centre3 = coreg_ec();
486         ycoreg_tv->centre2.el[0] = xcentre;
487         ycoreg_tv->centre2.el[1] = zcentre;
488         ycoreg_tv->ex3 = coreg_ex();
489         ycoreg_tv->ey3 = coreg_ez();
490         ycoreg_tv->ez3 = vec3_minus(coreg_ey());
491         ycoreg_tv->scalex = xscale/min_scale;
492         ycoreg_tv->scaley = zscale/min_scale;
493     }
494     if (xcoreg_tv!=NULL)
495     {
496         xcoreg_tv->centre3 = coreg_ec();
497         xcoreg_tv->centre2.el[0] = zcentre;
498         xcoreg_tv->centre2.el[1] = ycentre;
499         xcoreg_tv->ex3 = coreg_ez();
500         xcoreg_tv->ey3 = coreg_ey();
501         xcoreg_tv->ez3 = vec3_minus(coreg_ex());
502         xcoreg_tv->scalex = zscale/min_scale;
503         xcoreg_tv->scaley = yscale/min_scale;
504     }
505 }
506 
507 void xcoreg_tv_set(Tv *tv)
508 {
509     xcoreg_tv = tv;
510 }
511 
512 Tv *xcoreg_tv_get()
513 {
514    return(xcoreg_tv);
515 }
516 
517 void ycoreg_tv_set(Tv *tv)
518 {
519     ycoreg_tv = tv;
520 }
521 
522 Tv *ycoreg_tv_get()
523 {
524    return(ycoreg_tv);
525 }
526 
527 void zcoreg_tv_set(Tv *tv)
528 {
529     zcoreg_tv = tv;
530 }
531 
532 Tv *zcoreg_tv_get()
533 {
534    return(zcoreg_tv);
535 }
536 
537 static void cmap_set_anaglyph(Tv *tv)
538 {
539     if (tv && tv->tv_screen)
540       {
541           Bool (*cmap_create_fn) () = NULL;
542 
543           tv->cmap_create_fn = tv_anag_cmap_create;
544           tv_screen_cmap_find_and_install(tv, tv->tv_screen);
545           tv_screen_color_set(tv->tv_screen, tv->cmap_data_visible->std_lut[tv->color]);
546       }
547 }
548 
549 static void cmap_set_standard(Tv *tv)
550 {
551     if (tv && tv->tv_screen )
552       {
553           Bool (*cmap_create_fn) () = NULL;
554 
555           tv->cmap_create_fn = tv_standard_cmap_create;
556           tv_screen_cmap_find_and_install(tv, tv->tv_screen);
557           tv_screen_color_set(tv->tv_screen, tv->cmap_data_visible->std_lut[tv->color]);
558       }
559 }
560 
561 
562 void coreg_coords(Vec3 *nex, Vec3 *ney, Vec3 *nez)
563 {
564     *nex = coreg_ex();
565     *ney = coreg_ey();
566     *nez = coreg_ez();
567 }
568 
569 void coreg_tcoords(Vec3 *nex, Vec3 *ney, Vec3 *nez)
570 {
571     nex->el[0] = (float)ex[0];
572     nex->el[1] = (float)ey[0];
573     nex->el[2] = (float)ez[0];
574     ney->el[0] = (float)ex[1];
575     ney->el[1] = (float)ey[1];
576     ney->el[2] = (float)ez[1];
577     nez->el[0] = (float)ex[2];
578     nez->el[1] = (float)ey[2];
579     nez->el[2] = (float)ez[2];
580 }
581 
582 
583 void coreg_centre(double *x, double *y, double *z, Vec3*c)
584 {
585    *x = xcentre;
586    *y = ycentre;
587    *z = zcentre;
588    *c = coreg_ec();
589 }
590 
591 void coreg_scales(double *x, double *y, double *z, double *nx, double *ny, double *nz)
592 {
593     *x = xscale;
594     *y = yscale;
595     *z = zscale;
596     *nx = nxscale;
597     *ny = nyscale;
598     *nz = nzscale;
599 }
600 
601 
602 void set_coreg_zoom(double zfac, int mask)
603 {
604     if (mask >0 && (mask&4)) zoom_on = true; 
605     if (mask >0 && !(mask&4)) zoom_on = false; 
606     
607     zoom_fac = zfac;
608     min_scale = zoom_fac *MIN(MIN(xscale,yscale),zscale);
609     if(zcoreg_tv!=NULL) zcoreg_tv->scalex = xscale/min_scale;
610     if(zcoreg_tv!=NULL) zcoreg_tv->scaley = yscale/min_scale;
611     if(ycoreg_tv!=NULL) ycoreg_tv->scalex = xscale/min_scale;
612     if(ycoreg_tv!=NULL) ycoreg_tv->scaley = zscale/min_scale;
613     if(xcoreg_tv!=NULL) xcoreg_tv->scalex = zscale/min_scale;
614     if(xcoreg_tv!=NULL) xcoreg_tv->scaley = yscale/min_scale;
615 }
616 
617 void coreg_coords_init()
618 {
619     coreg_comp(ex, vec3_ex());
620     coreg_comp(ey, vec3_ey());
621     coreg_comp(ez, vec3_ez());
622     xscale = ixscale;
623     yscale = iyscale;
624     zscale = izscale;
625     if (disp_im_type == 0)
626         coreg_comp(ec, vec3(xcentre,ycentre,zcentre));
627     else
628         coreg_comp(ec, vec3(xcentre*nxscale/xscale,ycentre*nyscale/yscale,zcentre*nzscale/zscale));
629 }
630 
631 Vec3 coreg_bproj(double x, double y, double z)
632 /* compute equivalent location of transformed data point in original
633    data set */
634 {
635     Vec3 posi, post;
636     Vec3 ex,ey,ez;
637     coreg_tcoords(&ex,&ey,&ez);
638 
639     posi.el[0] = (x - xcentre)*nxscale;
640     posi.el[1] = (y - ycentre)*nyscale;
641     posi.el[2] = (z - zcentre)*nzscale;
642     post.el[0] = ec[0] + vec3_dot(posi,ex)/xscale;
643     post.el[1] = ec[1] + vec3_dot(posi,ey)/yscale;
644     post.el[2] = ec[2] + vec3_dot(posi,ez)/zscale;
645     return(post);
646 }
647 
648 
649 Vec3 coreg_proj(double x, double y, double z)
650 /* compute equivalent location of original data set in new coordinates */
651 {
652     Vec3 posi, post;
653 
654     post.el[0] = (x - ec[0])*xscale;
655     post.el[1] = (y - ec[1])*yscale;
656     post.el[2] = (z - ec[2])*zscale;
657 
658     posi.el[0] = vec3_dot(post,coreg_ex())/nxscale + xcentre;
659     posi.el[1] = vec3_dot(post,coreg_ey())/nyscale + ycentre;
660     posi.el[2] = vec3_dot(post,coreg_ez())/nzscale + zcentre;
661     
662     return(posi);
663 }
664 
665 static Imrect *imf_edges(Imrect *im)
666 {
667     static double sigma = 1.0, precision = 0.01;
668     static double low_edge_thres = 20.0, up_edge_thres = 35.0;
669     static int len_thres = 5;
670     Imrect *edgeim;
671 
672     edgeim = canny(im, sigma, precision, low_edge_thres, up_edge_thres,
673                    len_thres); 
674 
675     return(edgeim);
676 }
677 
678 
679 static void tv_crosshair(Tv * tv, Ipos pos)
680 {
681     int x, y;
682     int w = tv->width;
683     int h = tv->height;
684 
685     x = ipos_x(pos);
686     y = ipos_y(pos);
687     tv_linexy(tv, 0, y, w, y);
688     tv_linexy(tv, x, 0, x, h);
689 }
690 
691 
692 static Imrect *imf_interlace(Imrect *ima, Imrect *imb, Imregion *roi)
693 {
694     Imrect *im2,*im3,*im4,*im5,*im6,*im7;
695     im2 = imf_checquer(roi->ux-roi->lx,roi->uy-roi->ly,
696                        16 , 16);
697     imf_times_inplace(1.0/256.0,im2);
698     im3 = imf_times(-1.0,im2);
699     im4 = imf_add(1.0,im3);
700     im_free(im3);
701     im5 = imf_prod(im2, ima);
702     im_free(im2);
703     im6 = imf_prod(im4,imb);
704     im_free(im4);
705     im7 = imf_sum(im5,im6);
706     im_free(im5);
707     im_free(im6);
708     return(im7);
709 }
710 
711 
712 static Imrect *apply_threshold(Imrect *im1, int image_code)
713 {
714 /* Function to apply the threshold from the air parameters dialog box to the second image in the
715 anaglyph display mode */
716 
717         Imrect *im2=NULL;
718         Imregion *roi;
719         int lx, ux, ly, uy, i, j;
720         int tot_pix=0, pix_left=0;
721         float *row1, *row2, percentage, min, max;
722 
723         roi = im1->region;
724         lx = roi->lx;
725         ly = roi->ly;
726         ux = roi->ux;
727         uy = roi->uy;
728         row1 = fvector_alloc(lx, ux);
729         row2 = fvector_alloc(lx, ux);
730 
731         imf_minmax_nzero(im1, &min, &max);
732         im2 = im_alloc(im1->height, im1->width, roi, float_v);
733 
734         for(i=ly; i<uy; i++)
735         {
736                 im_get_rowf(row1, im1, i, lx, ux);
737                 for(j=lx; j<ux; j++)
738                 {
739                         tot_pix++;
740                         if(row1[j]>coreg_threshold)
741                         {
742                                 row2[j]=row1[j];
743                                 pix_left++;
744                         }       
745                         else
746                         {
747                                 row2[j] = min;
748                         }
749                 }
750                 im_put_rowf(row2, im2, i, lx, ux);
751         }
752         percentage = 100* pix_left/tot_pix;
753         if(image_code==0)
754         {
755                 format("Threshold retains %f%% of z image pixels.\n", percentage);
756         }
757         else if(image_code==1)
758         {
759                 format("Threshold retains %f%% of y image pixels.\n", percentage);
760         }
761         else
762         {
763                 format("Threshold retains %f%% of x image pixels.\n", percentage);
764         }
765         fvector_free(row1, lx);
766         fvector_free(row2, lx);
767         return im2;
768 }
769 
770 
771 static void zpartdraw(Tv * tv)
772 {
773     Imrect *im=NULL,*im1,*im2;
774     double storex,storey;
775     double min_nscale = zoom_fac*MIN(MIN(nxscale,nyscale),nzscale);
776     Imregion *roi;
777 
778     tv_erase(tv);
779     if (imptrs == NULL) return;
780     storex = tv->scalex;
781     storey = tv->scaley;
782     tv->scalex = nxscale/min_nscale;
783     tv->scaley = nyscale/min_nscale;
784     if (zcoreg_im==NULL)
785     {
786        im = seq_slicez(zcentre,NULL,coreg_bproj);
787        cmap_set_standard(tv);
788        im1 = imf_scale_nzero(im, 0.0, 255.0);
789        tv_imrect2(tv, im1);
790        im_free(im1);
791     }
792     else if(zcoreg_disp!=NULL && zcoreg_disp->vtype!= ptr_v)
793     {
794        roi = zcoreg_disp->region;
795        im = seq_slicez(zcentre,roi, coreg_bproj);
796        im1 = imf_scale_nzero(im, 0.0, 255.0); 
797        if (disp_im_type == 1)
798        {
799           cmap_set_anaglyph(tv);  
800           tv_imrect_anaglyph(tv, im1, zcoreg_disp, 0.0, 255.0);
801        }
802        else
803        {
804           im2 = imf_interlace(im1,zcoreg_disp, roi);
805           cmap_set_standard(tv);
806           tv_imrect2(tv, im2);
807           im_free(im2);
808        }
809        im_free(im1); 
810     }
811     else
812     {
813        im = seq_slicez(zcentre,NULL, coreg_bproj);
814        cmap_set_standard(tv);
815        im1 = imf_scale_nzero(im, 0.0, 255.0);
816        tv_imrect2(tv, im1);
817        im_free(im1);
818        tv_edges_conn(zcoreg_tv, zcoreg_disp);       
819     }
820     tv->scalex = storex;
821     tv->scaley = storey;
822     if (im!=NULL) im_free(im);
823     tv_color_set(tv, blue);
824     tv_crosshair(tv,ipos(tv->width/2,tv->height/2));
825     tv_reset_draw(tv);
826 }
827 
828 
829 static void ypartdraw(Tv * tv)
830 {
831     Imrect *im=NULL,*im1,*im2;
832     double storex,storey;
833     double min_nscale = zoom_fac*MIN(MIN(nxscale,nyscale),nzscale);
834     Imregion *roi;
835 
836     tv_erase(tv);
837     if (imptrs==NULL) return;
838     storex = tv->scalex;
839     storey = tv->scaley;
840     tv->scalex = nxscale/min_nscale;
841     tv->scaley = nzscale/min_nscale;
842     if (ycoreg_im==NULL)
843     {
844        im = seq_slicey(ycentre,NULL, coreg_bproj);
845        cmap_set_standard(tv);
846        im1 = imf_scale_nzero(im, 0.0, 255.0);
847        tv_imrect2(tv, im1);
848        im_free(im1);
849     }
850     else if(ycoreg_disp!=NULL && ycoreg_disp->vtype != ptr_v)
851     {
852        roi = ycoreg_disp->region;
853        im = seq_slicey(ycentre,roi, coreg_bproj);
854        im1 = imf_scale_nzero(im, 0.0, 255.0);
855        if (disp_im_type == 1)
856        {
857           cmap_set_anaglyph(tv);
858           tv_imrect_anaglyph(tv, im1, ycoreg_disp, 0.0, 255.0);
859        }
860        else
861        {
862           im2 = imf_interlace(im1,ycoreg_disp, roi);
863           cmap_set_standard(tv);
864           tv_imrect2(tv, im2);
865           im_free(im2);
866        }
867        im_free(im1);
868     }
869     else
870     {
871        im = seq_slicey(ycentre,NULL, coreg_bproj);
872        cmap_set_standard(tv);
873        im1 = imf_scale_nzero(im, 0.0, 255.0);
874        tv_imrect2(tv, im1);
875        im_free(im1);
876        tv_edges_conn(ycoreg_tv, ycoreg_disp);       
877     }
878     tv->scalex = storex;
879     tv->scaley = storey;
880     if (im!=NULL) im_free(im);
881     tv_color_set(tv, blue);
882     tv_crosshair(tv,ipos(tv->width/2,tv->height/2));
883     tv_reset_draw(tv);
884 }
885 
886 static void xpartdraw(Tv * tv)
887 {
888     Imrect *im=NULL,*im1,*im2;
889     double storex,storey;
890     double min_nscale = zoom_fac*MIN(MIN(nxscale,nyscale),nzscale);
891     Imregion *roi;
892 
893     tv_erase(tv);
894     if (imptrs == NULL) return;
895     storex = tv->scalex;
896     storey = tv->scaley;
897     tv->scalex = nzscale/min_nscale;
898     tv->scaley = nyscale/min_nscale;
899     if (xcoreg_im==NULL)
900     {
901        im = seq_slicex(xcentre,NULL, coreg_bproj);
902        im1 = imf_scale_nzero(im, 0.0, 255.0);
903        cmap_set_standard(tv);
904        tv_imrect2(tv, im1);
905        im_free(im1);
906     }
907     else if(xcoreg_disp!=NULL && xcoreg_disp->vtype != ptr_v )
908     {
909        roi = xcoreg_disp->region;
910        im = seq_slicex(xcentre,roi, coreg_bproj);
911        im1 = imf_scale_nzero(im, 0.0, 255.0);
912        if (disp_im_type == 1)
913        {
914           cmap_set_anaglyph(tv);
915           tv_imrect_anaglyph(tv, im1, xcoreg_disp, 0.0, 255.0);
916        }
917        else
918        {
919           im2 = imf_interlace(im1,xcoreg_disp, roi);
920           cmap_set_standard(tv);
921           tv_imrect2(tv, im2);
922           im_free(im2);
923        }
924        im_free(im1);
925        
926     }
927     else 
928     {
929        im = seq_slicex(xcentre,NULL, coreg_bproj);
930        cmap_set_standard(tv);
931        im1 = imf_scale_nzero(im, 0.0, 255.0);
932        tv_imrect2(tv, im1);
933        im_free(im1);
934        tv_edges_conn(xcoreg_tv, xcoreg_disp);
935     }
936     tv->scalex = storex;
937     tv->scaley = storey;
938     if (im!=NULL) im_free(im);
939     tv_color_set(tv, blue);
940     tv_crosshair(tv,ipos(tv->width/2,tv->height/2));
941     tv_reset_draw(tv);
942 }
943 
944 static void xfulldraw(Tv *tv)
945 {
946     double shiftcx;
947     double shiftcy;
948 /*
949     double min_scale = zoom_fac*MIN(MIN(xscale,yscale),zscale);
950 */
951     double min_nscale = zoom_fac*MIN(MIN(nxscale,nyscale),nzscale);
952     static int first= 0;
953 
954     shiftcx = ((float)(tv->width / 2) - tv->cx)*min_nscale/nzscale;
955     shiftcy = ((float)(tv->height / 2) - tv->cy)*min_nscale/nyscale;
956 
957     tv_erase(tv);
958     coreg_comp(ex, vec3_minus(tv->ez3));
959     coreg_comp(ey, tv->ey3);
960     coreg_comp(ez, tv->ex3);
961     tv->centre3 = coreg_bproj(xcentre,
962                               shiftcy+ycentre,
963                               shiftcx+zcentre);
964 
965     if (xcoreg_im==NULL)
966     {
967        zcentre = tv->centre3.el[2];
968        ycentre = tv->centre3.el[1];
969        tv->centre2.el[0] = zcentre;
970        tv->centre2.el[1] = ycentre;
971     }
972     tv->cx = (float)(tv->width / 2);
973     tv->cy = (float)(tv->height / 2);
974     coreg_comp(ec, tv->centre3);
975     if (first++>1) /* stop X from arbitrarily rescaling first time */
976     {
977        zscale = tv->scalex*min_scale;
978        yscale = tv->scaley*min_scale;
979     }
980     reset_tv_coords();
981     xpartdraw(tv);
982     if (ycoreg_tv!=NULL) ypartdraw(ycoreg_tv);   
983     if (zcoreg_tv!=NULL) zpartdraw(zcoreg_tv);   
984     reset_air();
985     reset_centre(xcentre,ycentre,zcentre);
986 }
987 
988 static void yfulldraw(Tv *tv)
989 {
990     double shiftcx;
991     double shiftcy;
992 /*
993     double min_scale = zoom_fac*MIN(MIN(xscale,yscale),zscale);
994 */
995     double min_nscale = zoom_fac*MIN(MIN(nxscale,nyscale),nzscale);
996     static int first = 0 ;
997 
998     shiftcx = ((float)(tv->width / 2) - tv->cx)*min_nscale/nxscale;
999     shiftcy = ((float)(tv->height / 2) - tv->cy)*min_nscale/nzscale;
1000 
1001     tv_erase(tv);
1002     coreg_comp(ex, tv->ex3);
1003     coreg_comp(ey, vec3_minus(tv->ez3));
1004     coreg_comp(ez, tv->ey3);
1005     tv->centre3 = coreg_bproj(shiftcx+xcentre
1006                             , ycentre
1007                             , shiftcy+zcentre);
1008     if (ycoreg_im==NULL)
1009     {
1010        xcentre = tv->centre3.el[0];
1011        zcentre = tv->centre3.el[2];
1012        tv->centre2.el[0] = xcentre;
1013        tv->centre2.el[1] = zcentre;
1014     }
1015 
1016     tv->cx = (float)(tv->width / 2);
1017     tv->cy = (float)(tv->height / 2);
1018     coreg_comp(ec , tv->centre3);
1019     if (first++>1) /* stop X from arbitrarily rescaling first time */
1020     {
1021        xscale = tv->scalex*min_scale;
1022        zscale = tv->scaley*min_scale;
1023     }
1024     reset_tv_coords();
1025     ypartdraw(tv);
1026     if (zcoreg_tv!=NULL) zpartdraw(zcoreg_tv);
1027     if (xcoreg_tv!=NULL) xpartdraw(xcoreg_tv);
1028     reset_air();
1029     reset_centre(xcentre,ycentre,zcentre);
1030 }
1031 
1032 static void zfulldraw(Tv *tv)
1033 {
1034     double shiftcx;
1035     double shiftcy;
1036 /*
1037     double min_scale = zoom_fac*MIN(MIN(xscale,yscale),zscale);
1038 */
1039     double min_nscale = zoom_fac*MIN(MIN(nxscale,nyscale),nzscale);
1040     static int first = 0;
1041 
1042     shiftcx = ((float)(tv->width / 2) - tv->cx)*min_nscale/nxscale;
1043     shiftcy = ((float)(tv->height / 2) - tv->cy)*min_nscale/nyscale;
1044 
1045     tv_erase(tv);
1046     coreg_comp(ex, tv->ex3);
1047     coreg_comp(ey, tv->ey3);
1048     coreg_comp(ez, tv->ez3);
1049     tv->centre3 = coreg_bproj(shiftcx+xcentre,
1050                               shiftcy+ycentre,
1051                               zcentre);
1052 
1053     if (zcoreg_im==NULL)
1054     {
1055        xcentre = tv->centre3.el[0];
1056        ycentre = tv->centre3.el[1];
1057        tv->centre2.el[0] = xcentre;
1058        tv->centre2.el[1] = ycentre;
1059     }
1060     tv->cx = (float)(tv->width / 2);
1061     tv->cy = (float)(tv->height / 2);
1062     coreg_comp(ec , tv->centre3);
1063     if (first++>1)/* stop X from arbitrarily rescaling first time */
1064     {
1065        xscale = tv->scalex*min_scale;
1066        yscale = tv->scaley*min_scale;
1067     }
1068     reset_tv_coords();
1069     zpartdraw(tv);
1070     if (ycoreg_tv!=NULL) ypartdraw(ycoreg_tv);
1071     if (xcoreg_tv!=NULL) xpartdraw(xcoreg_tv);
1072     reset_air();
1073     reset_centre(xcentre,ycentre,zcentre);
1074 }
1075 
1076 
1077 static Ipos proj3_orthx(Tv * tv, Vec3 v)
1078 {
1079     Ipos    pos = {Ipos_id};
1080     double  x, y;
1081     Vec3    vx;
1082 
1083     v = vec3_diff(v, tv->centre3);
1084     vx.el[0] = v.el[0]*ixscale;
1085     vx.el[1] = v.el[1]*iyscale;
1086     vx.el[2] = v.el[2]*izscale;
1087     x = tv->cx + tv->scalex*vec3_dot(vx, tv->ex3)/izscale;
1088     y = tv->cy + tv->scaley*vec3_dot(vx, tv->ey3)/iyscale;
1089     pos.x = ROUND(x);
1090     pos.y = ROUND(y);
1091     return (pos);
1092 }
1093 
1094 static Ipos proj3_orthy(Tv * tv, Vec3 v)
1095 {
1096     Ipos    pos = {Ipos_id};
1097     double  x, y;
1098     Vec3    vx;
1099 
1100     v = vec3_diff(v, tv->centre3);
1101     vx.el[0] = v.el[0]*ixscale;
1102     vx.el[1] = v.el[1]*iyscale;
1103     vx.el[2] = v.el[2]*izscale;
1104     x = tv->cx + tv->scalex*vec3_dot(vx, tv->ex3)/ixscale;
1105     y = tv->cy + tv->scaley*vec3_dot(vx, tv->ey3)/izscale;
1106     pos.x = ROUND(x);
1107     pos.y = ROUND(y);
1108     return (pos);
1109 }
1110 
1111 static Ipos proj3_orthz(Tv * tv, Vec3 v)
1112 {
1113     Ipos    pos = {Ipos_id};
1114     double  x, y;
1115     Vec3    vx;
1116 
1117     v = vec3_diff(v, tv->centre3);
1118     vx.el[0] = v.el[0] * ixscale;
1119     vx.el[1] = v.el[1] * iyscale;
1120     vx.el[2] = v.el[2] * izscale;
1121     x = tv->cx + tv->scalex*vec3_dot(vx, tv->ex3)/ixscale;
1122     y = tv->cy + tv->scaley*vec3_dot(vx, tv->ey3)/iyscale;
1123     pos.x = ROUND(x);
1124     pos.y = ROUND(y);
1125     return (pos);
1126 }
1127 
1128 static void linexy(Tv * tv, int x1, int y1, int x2, int y2)
1129 {
1130     tv_line3(tv, vec3((float) x1 , (float) y1 , tv->centre3.el[2]),
1131                  vec3((float) x2 , (float) y2 , tv->centre3.el[2]));
1132 }
1133 
1134 static void linexz(Tv * tv, int x1, int y1, int x2, int y2)
1135 {
1136     tv_line3(tv, vec3((float) x1 , tv->centre3.el[1], (float) y1),
1137                  vec3((float) x2 , tv->centre3.el[1], (float) y2));
1138 }
1139 
1140 static void linezy(Tv * tv, int x1, int y1, int x2, int y2)
1141 {
1142     tv_line3(tv, vec3(tv->centre3.el[0], (float) y1, (float)x1),
1143                  vec3(tv->centre3.el[0], (float) y2, (float)x2));
1144 }
1145 
1146 static void    zskeldraw(Tv * tv)
1147 {
1148     int     lx, ly, ux, uy;
1149     int     dx, dy;
1150     int     x, y;
1151 
1152     if (tv == NULL || tv->tv_screen == NULL)
1153         return;
1154 
1155     if (!zoom_on) set_coreg_zoom(zoom_fac,-1);
1156 
1157     imptrs = seq_limits(&imptrlz,&imptruz,&imptrly,
1158                           &imptruy,&imptrlx,&imptrux);
1159     lx = imptrlx ;
1160     ux = imptrux ;
1161     ly = imptrly ;
1162     uy = imptruy ;
1163 
1164     dx = (ux - lx) / 16;
1165     dy = (uy - ly) / 16;
1166     if (dx<1) dx = 1; /* BUG fix NAT */
1167     if (dy<1) dy = 1; /* 26/5/95     */
1168 
1169     for (x = lx; x < ux; x += dx)
1170         linexy(tv, x, ly, x, uy);
1171     linexy(tv, ux, ly, ux, uy);
1172     for (y = ly; y < uy; y += dy)
1173         linexy(tv, lx, y, ux, y);
1174     linexy(tv, lx, y, ux, y);
1175 }
1176 
1177 static void    yskeldraw(Tv * tv)
1178 {
1179     int     lx, ly, ux, uy;
1180     int     dx, dy;
1181     int     x, y;
1182 
1183     if (tv == NULL || tv->tv_screen == NULL)
1184         return;
1185 
1186     if (!zoom_on) set_coreg_zoom(zoom_fac,-1);
1187 
1188     imptrs = seq_limits(&imptrlz,&imptruz,&imptrly,
1189                           &imptruy,&imptrlx,&imptrux);
1190 
1191     lx = imptrlx ;
1192     ux = imptrux ;
1193     ly = imptrlz ;
1194     uy = imptruz ;
1195 
1196     dx = (ux - lx) / 16;
1197     dy = (uy - ly) / 16;
1198     if (dx<1) dx = 1; /* BUG fix NAT */
1199     if (dy<1) dy = 1; /* 26/5/95     */
1200 
1201     for (x = lx; x < ux; x += dx)
1202         linexz(tv, x, ly, x, uy);
1203     linexz(tv, ux, ly, ux, uy);
1204     for (y = ly; y < uy; y += dy)
1205         linexz(tv, lx, y, ux, y);
1206     linexz(tv, lx, y, ux, y);
1207 }
1208 
1209 static void    xskeldraw(Tv * tv)
1210 {
1211     int     lx, ly, ux, uy;
1212     int     dx, dy;
1213     int     x, y;
1214 
1215     if (tv == NULL || tv->tv_screen == NULL)
1216         return;
1217 
1218     if (!zoom_on) set_coreg_zoom(zoom_fac,-1);
1219 
1220     imptrs = seq_limits(&imptrlz,&imptruz,&imptrly,
1221                           &imptruy,&imptrlx,&imptrux);
1222     lx = imptrlz ;
1223     ux = imptruz ;
1224     ly = imptrly ;
1225     uy = imptruy ;
1226 
1227     dx = (ux - lx) / 16;
1228     dy = (uy - ly) / 16;
1229     if (dx<1) dx = 1; /* BUG fix NAT */
1230     if (dy<1) dy = 1; /* 26/5/95     */
1231 
1232     for (x = lx; x < ux; x += dx)
1233         linezy(tv, x, ly, x, uy);
1234     linezy(tv, ux, ly, ux, uy);
1235     for (y = ly; y < uy; y += dy)
1236         linezy(tv, lx, y, ux, y);
1237     linezy(tv, lx, y, ux, y);
1238 }
1239 
1240 static void zinit(Tv * tv)
1241 {
1242     Sequence *seq = seq_get_current();
1243     List *frame;
1244     frame = get_seq_start_el(seq);
1245     if (frame!=NULL) imref = frame->to;
1246 
1247     tv_camera3(tv, vec3(xcentre, ycentre, zcentre), 1.0, 3.0, vec3_ez(), vec3_ey());
1248     if (imref != NULL)
1249     {
1250        tv_camera2(tv, vec2(xcentre, ycentre), MAX(imref->width,imref->height)/2.0, vec2_ey());
1251     }
1252     else
1253     {
1254        tv_camera2(tv, vec2(xcentre, ycentre), 256, vec2_ey());
1255     }
1256     xscale = ixscale;
1257     yscale = iyscale;
1258     zscale = izscale;
1259     min_scale = zoom_fac *MIN(MIN(xscale,yscale),zscale);
1260     tv->proj3 = proj3_orthz;
1261     tv->scalex = xscale/min_scale;
1262     tv->scaley = yscale/min_scale;
1263 }
1264 
1265 static void yinit(Tv * tv)
1266 {
1267     Sequence *seq = seq_get_current();
1268     List *frame=NULL;
1269     frame = get_seq_start_el(seq);
1270     if (frame!=NULL) imref = frame->to;
1271 
1272     tv_camera3(tv, vec3(xcentre,ycentre,zcentre), 1.0, 3.0, vec3_minus(vec3_ey()), vec3_ez());
1273     if (imref != NULL)
1274     {
1275        tv_camera2(tv, vec2(xcentre, zcentre), MAX(imref->height,imref->width)/2.0, vec2_ey());
1276     }
1277     else
1278     {
1279        tv_camera2(tv, vec2(xcentre, zcentre), 256, vec2_ey());
1280     }
1281     xscale = ixscale;
1282     yscale = iyscale;
1283     zscale = izscale;
1284     min_scale = zoom_fac*MIN(MIN(xscale,yscale),zscale);
1285     tv->proj3 = proj3_orthy;
1286     tv->scalex = xscale/min_scale;
1287     tv->scaley = zscale/min_scale;
1288 }
1289 
1290 static void xinit(Tv * tv)
1291 {
1292     Sequence *seq = seq_get_current();
1293     List *frame=NULL;
1294     frame = get_seq_start_el(seq);
1295     if (frame!=NULL) imref = frame->to;
1296 
1297     tv_camera3(tv, vec3(xcentre,ycentre,zcentre), 1.0, 3.0, vec3_minus(vec3_ex()), vec3_ey());
1298     if (imref != NULL)
1299     {
1300        tv_camera2(tv, vec2(zcentre, ycentre), MAX(imref->height,imref->width)/2.0, vec2_ey());
1301     }
1302     else
1303     {
1304        tv_camera2(tv, vec2(zcentre, ycentre), 256, vec2_ey());
1305     }
1306     xscale = ixscale;
1307     yscale = iyscale;
1308     zscale = izscale;
1309     min_scale = zoom_fac*MIN(MIN(xscale,yscale),zscale);
1310     tv->proj3 = proj3_orthx;
1311     tv->scalex = zscale/min_scale;
1312     tv->scaley = yscale/min_scale;
1313 }
1314 
1315 
1316 void coreg_init(Tv *tv)
1317 {
1318     if(zcoreg_tv!=NULL) zinit(zcoreg_tv);
1319     if(ycoreg_tv!=NULL) yinit(ycoreg_tv);
1320     if(xcoreg_tv!=NULL) xinit(xcoreg_tv);
1321 }
1322 
1323 void coreg_redraw()
1324 {
1325     if(zcoreg_tv!=NULL) zfulldraw(zcoreg_tv);
1326     else if(ycoreg_tv!=NULL) yfulldraw(ycoreg_tv);
1327     else if(xcoreg_tv!=NULL) xfulldraw(xcoreg_tv);
1328 }
1329 
1330 Tv     *xcoreg_tv_make(void)
1331 {
1332     static Tv *tv;
1333 
1334     tv = tv_create("Xcoreg");
1335     xcoreg_tv = tv;
1336     xinit(tv);
1337     (void) tv_set_backdraw(tv, xfulldraw);
1338     (void) tv_set_skeldraw(tv, xskeldraw);
1339     tv_set_init(tv, coreg_init);
1340     (void) tv_set_zoomlevel(tv, IMZOOM);
1341     return (tv);
1342 }
1343 
1344 Tv     *ycoreg_tv_make(void)
1345 {
1346     static Tv *tv;
1347 
1348     tv = tv_create("Ycoreg");
1349     ycoreg_tv = tv;
1350     yinit(tv);
1351     (void) tv_set_backdraw(tv, yfulldraw);
1352     (void) tv_set_skeldraw(tv, yskeldraw);
1353     tv_set_init(tv, coreg_init);
1354     (void) tv_set_zoomlevel(tv, IMZOOM);
1355     
1356     return (tv);
1357 }
1358 
1359 Tv     *zcoreg_tv_make(void)
1360 {
1361     static Tv *tv;
1362 
1363     tv = tv_create("Zcoreg");
1364     zcoreg_tv = tv;
1365     zinit(tv);
1366     (void) tv_set_backdraw(tv, zfulldraw);
1367     (void) tv_set_skeldraw(tv, zskeldraw);
1368     tv_set_init(tv, coreg_init);
1369     (void) tv_set_zoomlevel(tv, IMZOOM);
1370     return (tv);
1371 }
1372 
1373 void set_coreg_centre(double slicex, double slicey, double slicez)
1374 {
1375     xcentre = slicex;
1376     ycentre = slicey;
1377     zcentre = slicez;
1378     ec[0] = xcentre;
1379     ec[1] = ycentre;
1380     ec[2] = zcentre;
1381 }
1382 
1383 void set_coreg_trans(struct air16  air1)
1384 {
1385     Vec3 posi,post;
1386     Vec3 nex,ney,nez;
1387     double min_nscale;
1388 
1389     nxscale = air1.s.x_size;
1390     nyscale = air1.s.y_size;
1391     nzscale = air1.s.z_size;
1392     xscale = air1.r.x_size;
1393     yscale = air1.r.y_size;
1394     zscale = air1.r.z_size;
1395     ixscale = xscale;
1396     iyscale = yscale;
1397     izscale = zscale;
1398     min_nscale = MIN(MIN(nxscale,nyscale),nzscale);
1399     min_scale = MIN(MIN(xscale,yscale),zscale);
1400  
1401     ex[0] = air1.e[0][0]*air1.s.x_size/min_nscale;  
1402     ex[1] = air1.e[0][1]*air1.s.y_size/min_nscale;  
1403     ex[2] = air1.e[0][2]*air1.s.z_size/min_nscale;  
1404     ey[0] = air1.e[1][0]*air1.s.x_size/min_nscale;  
1405     ey[1] = air1.e[1][1]*air1.s.y_size/min_nscale;  
1406     ey[2] = air1.e[1][2]*air1.s.z_size/min_nscale;  
1407     ez[0] = air1.e[2][0]*air1.s.x_size/min_nscale;  
1408     ez[1] = air1.e[2][1]*air1.s.y_size/min_nscale;  
1409     ez[2] = air1.e[2][2]*air1.s.z_size/min_nscale;  
1410     post.el[0] = air1.e[3][0]*air1.s.x_size/min_nscale;
1411     post.el[1] = air1.e[3][1]*air1.s.y_size/min_nscale;
1412     post.el[2] = air1.e[3][2]*air1.s.z_size/min_nscale;
1413  
1414 /* now transform centre location */
1415     coreg_tcoords(&nex,&ney,&nez);
1416  
1417     posi.el[0] = xcentre*nxscale;
1418     posi.el[1] = ycentre*nyscale;
1419     posi.el[2] = zcentre*nzscale;
1420     ec[0]  = (post.el[0] + vec3_dot(posi,nex))/xscale;  
1421     ec[1]  = (post.el[1] + vec3_dot(posi,ney))/yscale;  
1422     ec[2]  = (post.el[2] + vec3_dot(posi,nez))/zscale;  
1423     reset_tv_coords();
1424     coreg_redraw();
1425 }
1426     
1427 
1428 void get_coreg_trans(struct air16 *air1)
1429 {
1430     Vec3 posi,post;
1431     Vec3 nex,ney,nez;
1432     double min_scale = MIN(MIN(nxscale,nyscale),nzscale);
1433  
1434     air1->e[0][0] = ex[0]*min_scale/nxscale;
1435     air1->e[0][1] = ex[1]*min_scale/nyscale;
1436     air1->e[0][2] = ex[2]*min_scale/nzscale;
1437     air1->e[1][0] = ey[0]*min_scale/nxscale;
1438     air1->e[1][1] = ey[1]*min_scale/nyscale;
1439     air1->e[1][2] = ey[2]*min_scale/nzscale;
1440     air1->e[2][0] = ez[0]*min_scale/nxscale;
1441     air1->e[2][1] = ez[1]*min_scale/nyscale;
1442     air1->e[2][2] = ez[2]*min_scale/nzscale;
1443     air1->s.x_size = nxscale;
1444     air1->s.y_size = nyscale;
1445     air1->s.z_size = nzscale;
1446     air1->r.x_size = xscale;
1447     air1->r.y_size = yscale;
1448     air1->r.z_size = zscale;
1449  
1450 /* now transform centre location */
1451     coreg_tcoords(&nex,&ney,&nez);
1452  
1453     posi.el[0] = -xcentre*nxscale;
1454     posi.el[1] = -ycentre*nyscale;
1455     posi.el[2] = -zcentre*nzscale;
1456     post.el[0] = ec[0]*xscale + vec3_dot(posi,nex);
1457     post.el[1] = ec[1]*yscale + vec3_dot(posi,ney);
1458     post.el[2] = ec[2]*zscale + vec3_dot(posi,nez);
1459     air1->e[3][0] = post.el[0]*min_scale/nxscale;
1460     air1->e[3][1] = post.el[1]*min_scale/nyscale;
1461     air1->e[3][2] = post.el[2]*min_scale/nzscale;
1462 }
1463  
1464 
1465 float display_scaler(Imrect *pre_thresh_im, Imrect *post_thresh_im)
1466 {
1467         /* Determines the scaling parameters for imf_scale_nzero such that 
1468         the same scaling that would have been applied to the pre threshold
1469         image is applied to the post threshold image for dispaly in the tv */
1470 
1471         float old_min, old_max, new_min, new_max, old_scale, new_scale_low;
1472 
1473         imf_minmax_nzero(pre_thresh_im, &old_min, &old_max);
1474         imf_minmax_nzero(post_thresh_im, &new_min, &new_max);
1475         
1476         old_scale = (256.0)/(old_max-old_min);
1477         new_scale_low = (new_min-old_min)*old_scale;
1478         return(new_scale_low);
1479 }
1480 
1481 
1482 void latch_slices(int im_type)
1483 {
1484     Imrect *im1, *im2, *im3;
1485     extern Imregion *tv_get_im_roi(); 
1486     Imregion *roi;
1487     
1488     /* PAB 21/04/2004: memory leak fixed: tv_get_im_roi and seq_slicex,y,z (calls im_alloc) both alloc a roi */
1489 
1490     disp_im_type = im_type;
1491     if (imptrs!=NULL)
1492     {
1493       if(zcoreg_im==NULL)
1494       {
1495          nxscale = xscale;
1496          nyscale = yscale;
1497          nzscale = zscale;
1498          roi = tv_get_im_roi(zcoreg_tv_get()); 
1499          zcoreg_im = seq_slicez(zcentre,roi,coreg_bproj);
1500          rfree((void *)roi);
1501          roi = tv_get_im_roi(ycoreg_tv_get());
1502          ycoreg_im = seq_slicey(ycentre,roi,coreg_bproj);
1503          rfree((void *)roi);
1504          roi = tv_get_im_roi(xcoreg_tv_get());
1505          xcoreg_im = seq_slicex(xcentre,roi,coreg_bproj);
1506          rfree((void *)roi);
1507       }
1508       if (zcoreg_disp!=NULL) im_free(zcoreg_disp);
1509       if (ycoreg_disp!=NULL) im_free(ycoreg_disp);
1510       if (xcoreg_disp!=NULL) im_free(xcoreg_disp);
1511       zcoreg_disp=NULL;
1512       ycoreg_disp=NULL;
1513       xcoreg_disp=NULL;
1514       if (im_type==2)
1515       {
1516          zcoreg_disp = imf_edges(zcoreg_im);
1517          ycoreg_disp = imf_edges(ycoreg_im);
1518          xcoreg_disp = imf_edges(xcoreg_im);
1519       }
1520       else if(im_type==1||im_type==3)
1521       {
1522          im1 = apply_threshold(zcoreg_im, 0);
1523          /* scale_low = display_scaler(zcoreg_im, im1);  */
1524          zcoreg_disp = imf_scale_nzero(im1, 0.0, 255.0);
1525          im_free(im1);
1526 
1527          im2 = apply_threshold(ycoreg_im, 1);
1528          /* scale_low = display_scaler(ycoreg_im, im2);  */
1529          ycoreg_disp = imf_scale_nzero(im2, 0.0, 255.0);
1530          im_free(im2);
1531 
1532          im3 = apply_threshold(xcoreg_im, 2);
1533          /* scale_low = display_scaler(xcoreg_im, im3);  */
1534          xcoreg_disp = imf_scale_nzero(im3, 0.0, 255.0);
1535          im_free(im3);
1536       }
1537     }
1538     coreg_init(NULL);
1539     coreg_redraw();
1540 }
1541 
1542 
1543 
1544 void ***coreg_slice_init(Sequence *seq, Vec3 *newiscale)
1545 {
1546     List *frame;
1547     Imregion roi;
1548 
1549     imptrs = seq_slice_init(seq);
1550 
1551     if (newiscale !=NULL)
1552     {
1553        ixscale = xscale = newiscale->el[0];
1554        iyscale = yscale = newiscale->el[1];
1555        izscale = zscale = newiscale->el[2];
1556     }
1557     frame = get_seq_start_el(seq);
1558     if((frame==NULL)||(imref = frame->to)==NULL) return(NULL);
1559     roi = *(imref->region);
1560     seq_limits(&imptrlz,&imptruz,&imptrly,
1561                  &imptruy,&imptrlx,&imptrux);
1562 
1563     if (zcoreg_tv!=NULL) tv_set_roi(zcoreg_tv,&roi);
1564     roi.ly = imptrlz;
1565     roi.uy = imptruz;
1566     if (ycoreg_tv!=NULL) tv_set_roi(ycoreg_tv,&roi);
1567     roi.ly = imptrly;
1568     roi.uy = imptruy;
1569     roi.lx = imptrlz;
1570     roi.ux = imptruz;
1571     if (xcoreg_tv!=NULL) tv_set_roi(xcoreg_tv,&roi);
1572     coreg_init(NULL);
1573 
1574     return(imptrs);
1575 }
1576 
1577 
1578 void latch_clear()
1579 {
1580     Sequence *seq = (Sequence *)seq_get_current();
1581     if (zcoreg_im!=NULL) im_free(zcoreg_im);
1582     zcoreg_im = NULL;
1583     if (ycoreg_im!=NULL) im_free(ycoreg_im);
1584     ycoreg_im = NULL;
1585     if (xcoreg_im!=NULL) im_free(xcoreg_im);
1586     xcoreg_im = NULL;
1587     if (zcoreg_disp!=NULL) im_free(zcoreg_disp);
1588     zcoreg_disp = NULL;
1589     if (ycoreg_disp!=NULL) im_free(ycoreg_disp);
1590     ycoreg_disp = NULL;
1591     if (xcoreg_disp!=NULL) im_free(xcoreg_disp);
1592     xcoreg_disp = NULL;
1593     coreg_slice_init(seq, seq_dim_to_vec3());
1594     coreg_coords_init();
1595     nxscale = xscale;
1596     nyscale = yscale;
1597     nzscale = zscale;
1598     set_coreg_zoom(zoom_fac, -1);
1599     coreg_init(NULL);
1600     coreg_redraw();
1601 }
1602 

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