home *** CD-ROM | disk | FTP | other *** search
/ Encyclopedia of Graphics File Formats Companion / GFF_CD.ISO / software / unix / saoimage / sao1_07.tar / crdsynth.c < prev    next >
C/C++ Source or Header  |  1990-04-20  |  9KB  |  258 lines
  1. #ifndef lint
  2. static char SccsId[] = "%W%  %G%";
  3. #endif
  4.  
  5. /* Module:    crdsynth.c (Coordinate Synthesis)
  6.  * Purpose:    Support the computing of transform matrices and shortcuts
  7.  * Subroutine:    invert_transform()        returns: void
  8.  * Subroutine:    compute_iadd()            returns: void
  9.  * Subroutine:    combine_transform()        returns: void
  10.  * Subroutine:    set_trans_speed()        returns: void
  11.  * Xlib calls:    none
  12.  * Copyright:    1988 Smithsonian Astrophysical Observatory
  13.  *        You may do anything you like with this file except remove
  14.  *        this copyright.  The Smithsonian Astrophysical Observatory
  15.  *        makes no representations about the suitability of this
  16.  *        software for any purpose.  It is provided "as is" without
  17.  *        express or implied warranty.
  18.  * Modified:    {0} Michael VanHilst    initial version            16 March 1988
  19.  *        {n} <who> -- <does what> -- <when>
  20.  */
  21.  
  22. #include <stdio.h>        /* stderr, NULL, etc */
  23. #include <math.h>        /* get trig functions, sqrt, and fabs */
  24. #include "hfiles/coord.h"    /* coord structs */
  25.  
  26. /*
  27.  * Subroutine:    invert_transform
  28.  * Purpose:    Create matrix to perform opposite transformation
  29.  *        Given AtoB (old) and Bsys ioff, compute BtoA (new)
  30.  */
  31. void invert_transform ( new, old, ioff )
  32.      Transform *new, *old;    /* pointers for new and original transforms */
  33.      float ioff;        /* integer coord to real translation factor */
  34. {
  35.   void invert_matrix(), compute_iadd_invert(), set_trans_speed();  
  36.  
  37. #ifdef DEBUG
  38.   void exit_errmsg();
  39.   /* check for undifined coordinate transforms */
  40.   if( ((old->inx_outx ==0) && (old->iny_outx ==0)) ||
  41.       ((old->inx_outy ==0) && (old->iny_outy ==0)) ) {
  42.     exit_errmsg("indeterminate coordinate transform!");
  43.   }
  44. #endif
  45.   /* check for simplified cases of orthogonal transform */
  46.   if( (old->iny_outx ==0) && (old->inx_outy ==0) ) {
  47.     new->iny_outx = 0.0;
  48.     new->inx_outy = 0.0;
  49.     new->inx_outx = 1.0 / old->inx_outx;
  50.     new->iny_outy = 1.0 / old->iny_outy;
  51.     new->add_outx = -old->add_outx * new->inx_outx;
  52.     new->add_outy = -old->add_outy * new->iny_outy;
  53.     new->iadd_outx = (-old->add_outx + ioff) * new->inx_outx;
  54.     new->iadd_outy = (-old->add_outy + ioff) * new->iny_outy;
  55.   } else if( (old->inx_outx ==0) && (old->iny_outy ==0) ) {
  56.     new->inx_outx = 0.0;
  57.     new->iny_outy = 0.0;
  58.     new->iny_outx = 1.0 / old->inx_outy;
  59.     new->inx_outy = 1.0 / old->iny_outx;
  60.     new->add_outx = -old->add_outy * new->iny_outx;
  61.     new->add_outy = -old->add_outx * new->inx_outy;
  62.     new->iadd_outx = (-old->add_outy + ioff) * new->iny_outx;
  63.     new->iadd_outy = (-old->add_outx + ioff) * new->inx_outy;
  64.   } else {
  65.     /* compute parameters by matrix inversion */
  66.     invert_matrix(old, new);
  67.     compute_iadd_invert(old, new, ioff);
  68.   }
  69.   /* set rotation flag */
  70.   new->no_rot = old->no_rot;
  71.   /* set parameters for speedy integer calculations */
  72.   set_trans_speed(new);
  73. }
  74.  
  75. #ifdef NOTNEEDED /* %% not currently needed */
  76. /*
  77.  * Subroutine:    compute_iadd
  78.  * Purpose:    Compute the offsets used for integer transforms
  79.  */
  80. void compute_iadd ( old, new, ioff )
  81.      Transform *old, *new;
  82.      float ioff;
  83. {
  84.   void compute_iadd_invert();
  85.  
  86.   /* check for simplified cases of orthogonal transform */
  87.   if( (new->iny_outx ==0) && (new->inx_outy ==0) ) {
  88.     new->iadd_outx = new->add_outx + (ioff * new->inx_outx);
  89.     new->iadd_outy = new->add_outy + (ioff * new->iny_outy);
  90.   } else if( (new->inx_outx ==0) && (new->iny_outy ==0) ) {
  91.     new->iadd_outx = new->add_outx + (ioff * new->iny_outx);
  92.     new->iadd_outy = new->add_outy + (ioff * new->inx_outy);
  93.   } else {
  94.   /* if no simple solution, carry offset through from inverse transform */
  95.     compute_iadd_invert(old, new, ioff);
  96.   }
  97. }
  98. #endif
  99.  
  100. /*
  101.  * Subroutine:    combine_transform
  102.  * combine two sets of transform parameters into a single set
  103.  * first and second apply to sequence starting with the input values
  104.  * algorithm is 3x3 matrix multiplication (3rd row is 0 0 1)
  105.  * combination[i][j] is sum of products of column j of 1st and row i of 2nd
  106.  */
  107. void combine_transform ( new, first, second )
  108.      Transform *new, *first, *second;
  109. {
  110.   void set_trans_speed();
  111.  
  112.   new->inx_outx = ((second->inx_outx * first->inx_outx) +
  113.            (second->iny_outx * first->inx_outy));
  114.   new->iny_outx = ((second->inx_outx * first->iny_outx) +
  115.            (second->iny_outx * first->iny_outy));
  116.   new->add_outx = ((second->inx_outx * first->add_outx) +
  117.            (second->iny_outx * first->add_outy) +
  118.            second->add_outx);
  119.   new->inx_outy = ((second->inx_outy * first->inx_outx) +
  120.            (second->iny_outy * first->inx_outy));
  121.   new->iny_outy = ((second->inx_outy * first->iny_outx) +
  122.            (second->iny_outy * first->iny_outy));
  123.   new->add_outy = ((second->inx_outy * first->add_outx) +
  124.            (second->iny_outy * first->add_outy) +
  125.            second->add_outy);
  126.   /* compute offset factor for transform starting with an integer value */
  127.   /* use integer input offset of first and float offset of second */
  128.   new->iadd_outx = ((second->inx_outx * first->iadd_outx) +
  129.             (second->iny_outx * first->iadd_outy) +
  130.             second->add_outx);
  131.   new->iadd_outy = ((second->inx_outy * first->iadd_outx) +
  132.             (second->iny_outy * first->iadd_outy) +
  133.             second->add_outy);
  134.  
  135.   /* set rotation flag, most cases don't have cross x-y relationships */
  136.   new->no_rot = first->no_rot && second->no_rot;
  137.   /* set up parameters needed for speedy calculation of transform */
  138.   /* uses integer multiply or divide on orthogonal rotations */
  139.   set_trans_speed(new);
  140.   return;
  141. }
  142.  
  143. #define EPSILON 0.0001
  144. #define NEGONE -0.9999
  145. #define CLOSE(a,b) (fabs((double)((a)-((float)(b))))<EPSILON)
  146. /*
  147.  * Subroutine:    set_trans_speed
  148.  * Purpose:    Set parameters for fast integer computation
  149.  */
  150. void set_trans_speed ( trans )
  151.      Transform *trans;
  152. {
  153.   int xzm, yzm;
  154.   static int integer_test();
  155.  
  156.   trans->ixzoom = 0;
  157.   trans->iyzoom = 0;
  158.   /* is it an unflipped transform? */
  159.   if( CLOSE(trans->iny_outx, 0) && CLOSE(trans->inx_outy, 0) ) {
  160.     /* straight forward x->x, y->y transform */
  161.     trans->no_rot = 1;
  162.     trans->flip = 0;
  163.     xzm = integer_test(trans->inx_outx, &trans->ixzoom);
  164.     yzm = integer_test(trans->iny_outy, &trans->iyzoom);
  165.   } else {
  166.     trans->no_rot = 0;
  167.     if( CLOSE(trans->inx_outx, 0) && CLOSE(trans->iny_outy, 0) ) {
  168.       trans->flip = 1;
  169.       xzm = integer_test(trans->iny_outx, &trans->ixzoom);
  170.       yzm = integer_test(trans->inx_outy, &trans->iyzoom);
  171.     } else {
  172.       trans->int_math = 0;
  173.       return;
  174.     }
  175.   }
  176.   /* if both zooms can be handled by integer operations ... */
  177.   /* ... when given integer inputs, use faster operations */
  178.   if( trans->ixzoom && trans->iyzoom ) {
  179.     if( (trans->ixzoom == 1) && (trans->iyzoom == 1) ) {
  180.       /* if both are zoom 1 */
  181.       trans->int_math = 1;
  182.       trans->zoom = 0;
  183.     } else if( (xzm >= 0) && (yzm >= 0) ) {
  184.       /* if both zooms can be integer multiplies, use integer multiply */
  185.       trans->int_math = 1;
  186.       trans->multiply = 1;
  187.       trans->zoom = 1;
  188.     } else if( (xzm <= 0) && (yzm <= 0) ) {
  189.       /* if both zooms can be integer divides, use float divide */
  190.       trans->int_math = 1;
  191.       trans->multiply = 0;
  192.       trans->zoom = 1;
  193.     } else {
  194.       /* if its not so simple (mixed operations) don't bother using shortcut */
  195.       trans->int_math = 0;
  196.     }
  197.   } else {
  198.     trans->int_math = 0;
  199.   }
  200. }
  201.  
  202. /*
  203.  * Subroutine:    integer_test
  204.  * Returns:    0 if not possible, 1 for zoom up, -1 for zoom down
  205.  * Purpose:    Test zoom to determine whether it could be computed
  206.  *        with integer math
  207.  */
  208. static int integer_test ( fzoom, izoom )
  209.      double fzoom;    /* i: zoom factor (>1, ==1, or <1) */
  210.      int *izoom;    /* o: integer zoom factor (>1 or ==1) */
  211. {
  212.   int ival;
  213.  
  214.   /* test for unity scaling */
  215.   if( CLOSE(fzoom, 1) ) {
  216.     *izoom = 1;
  217.     return( 0 );
  218.   }
  219.   /* test for even integer multiplies */
  220.   if( fzoom > 1.0 ) {
  221.     /* round and test for close match */
  222.     ival = fzoom + 0.5;
  223.     if( CLOSE(fzoom, ival) ) {
  224.       *izoom = ival;
  225.       return( 1 );
  226.     }
  227.   } else if( fzoom < NEGONE ) {
  228.     ival = fzoom - 0.5;
  229.     if( CLOSE(fzoom, ival) ) {
  230.       *izoom = ival;
  231.       return( 1 );
  232.     }
  233.   } else if( fzoom != 0.0 ) {
  234.     /* test for even integer divides */
  235.     fzoom = 1.0 / fzoom;
  236.     if( fzoom > 0 ) {
  237.       /* round and test for close match */
  238.       ival = fzoom + 0.5;
  239.       if( CLOSE(fzoom, ival) ) {
  240.     *izoom = ival;
  241.     return( -1 );
  242.       }
  243.     } else {
  244.       ival = fzoom - 0.5;
  245.       if( CLOSE(fzoom, ival) ) {
  246.     *izoom = ival;
  247.     return( -1 );
  248.       }
  249.     }
  250.   }
  251. #ifdef DEBUG
  252.   else
  253.     (void)fprintf(stderr, "Warning: Zero zoom factor (CoordSynth.c)\n");
  254. #endif
  255.   *izoom = 0;
  256.   return( 0 );
  257. }
  258.