fish/surf/mapcalc.c

   1 #include <aros/oldprograms.h>
   2 #include <stdio.h>
   3 #include "mytypes.h"
   4 #include "revolve.h"       /* need to get scrnpair from here */
   5 #include "mapstuff.h"
   6 #include "menuexp.h"
   7
   8 #ifdef TEST
   9 #undef DebugOn
  10 #define DebugOn 1
  11 #endif /* TEST */
  12
  13 extern struct Library * MathBase;
  14 extern struct Library * MathTransBase;
  15
  16 #define DEBUG
  17 /*
  18  * this section of code derived from:
  19  * "The essential algorithms of ray tracing" by Eric Haines
  20  * presented in Sigraph proceedings on Ray Tracing
  21  * my major change has been to simplify it for two dimensions
  22  */
  23
  24 typedef struct {
  25     float x, y;
  26 } Vector;
  27
  28 static float DotVector(a,b)
  29     Vector *a, *b;
  30 {
  31     return( a->x * b->x + a->y * b->y);
  32 }
  33
  34 static void DivVector(in, scale, out)
  35     Vector *in, *out;
  36     float scale;
  37 {
  38     if ( fabs(scale) < SingleTinyVal ) {
  39         out->x = SingleLargeVal;
  40         out->y = SingleLargeVal;
  41     }
  42     else {
  43         out->x = in->x / scale;
  44         out->y = in->y / scale;
  45     }
  46 }
  47
  48
  49
  50 static Vector Na, Nb, Nc;
  51 static float Du0, Du1, Du2,
  52              Dv0, Dv1, Dv2;
  53 static Vector Qux, Quy,
  54               Qvx, Qvy;
  55 static float Dux, Duy,
  56              Dvx, Dvy;
  57 static bool IsQuadu, IsQuadv;
  58
  59 void CalcMapConsts(vp)
  60     register ScrnPair *vp;
  61 #define p00 vp[0]
  62 #define p01 vp[1]
  63 #define p11 vp[2]
  64 #define p10 vp[3]
  65 {
  66     Vector Pa, Pb, Pc, Pd;
  67
  68     Pa.x = p00.x - p10.x + p11.x - p01.x;
  69     Pa.y = p00.y - p10.y + p11.y - p01.y;
  70
  71     Pb.x = p10.x - p00.x;
  72     Pb.y = p10.y - p00.y;
  73
  74     Pc.x = p01.x - p00.x;
  75     Pc.y = p01.y - p00.y;
  76
  77     Pd.x = p00.x;
  78     Pd.y = p00.y;
  79
  80     Na.x = Pa.y; Na.y = -Pa.x;
  81     Nc.x = Pc.y; Nc.y = -Pc.x;
  82     Nb.x = Pb.y; Nb.y = -Pb.x;
  83
  84     Du0 = DotVector(&Nc, &Pd);
  85     Du1 = DotVector(&Na, &Pd) + DotVector(&Nc, &Pb);
  86     Du2 = DotVector( &Na, &Pb);
  87
  88     if( fabs( Du2 ) > SingleTinyVal ) {
  89         float TwoDu2;
  90
  91         IsQuadu = true;
  92         TwoDu2 = 2.0 * Du2;
  93         DivVector( &Na, TwoDu2, &Qux );
  94         DivVector( &Nc, -Du2, &Quy );
  95         Duy = Du0/Du2;
  96         Dux = -Du1/TwoDu2;
  97     }
  98     else {
  99         IsQuadu = false;
 100     }
 101
 102     Dv0 = DotVector( &Nb, &Pd);
 103     Dv1 = DotVector(&Na, &Pd) + DotVector(&Nb, &Pc);
 104     Dv2 = DotVector( &Na, &Pc);
 105     if( fabs( Dv2 ) > SingleTinyVal ) {
 106         float TwoDv2;
 107
 108         IsQuadv = true;
 109         TwoDv2 = 2.0 * Dv2;
 110         DivVector( &Na, TwoDv2, &Qvx);
 111         DivVector( &Nb, -Dv2, &Qvy);
 112 /*      DivVector( &Nc, -Dv2, &Qvy); */
 113         Dvx = - Dv1/TwoDv2;
 114         Dvy = Dv0/Dv2;
 115     }
 116     else {
 117         IsQuadv = false;
 118     }
 119 #ifdef DEBUG
 120     if( DebugOn ) {
 121         printf("du2 %f, du1 %f, du0 %f\n", Du2, Du1, Du0 );
 122         printf("dv2 %f, dv1 %f, dv0 %f\n", Dv2, Dv1, Dv0 );
 123         printf("Na = (%f, %f), Nb = (%f,%f), Nc = (%f,%f)\n",
 124         Na.x, Na.y, Nb.x, Nb.y, Nc.x, Nc.y );
 125         printf("IsQuad =(%c, %c)\n", IsQuadu?'t':'f', IsQuadv? 't': 'f' );
 126     }
 127 #endif /* DEBUG */
 128
 129 }
 130
 131
 132 /*
 133  * given points px,py in the quadrilateral, map them to points inside a
 134  * unit square
 135  */
 136 void  MapXYRatio(px, py, outx, outy, SweepCode)
 137     float px, py;
 138     float *outx, *outy;
 139     short SweepCode;
 140 {
 141     float resu, resv;
 142     Vector Ri;
 143
 144
 145     Ri.x = px; Ri.y = py;
 146
 147     if( !IsQuadu ) {
 148         float denom;
 149         denom = (Du1 - DotVector(&Na, &Ri));
 150         if( fabs(denom) < SingleTinyVal )
 151             resu = 2.0;
 152         else
 153             resu = (DotVector(&Nc, &Ri) - Du0)/denom;
 154     } else {
 155         float Ka, Kb;
 156         float discrim;
 157
 158         Ka = Dux + DotVector( &Qux, &Ri);
 159         Kb = Duy + DotVector( &Quy, &Ri);
 160         discrim = sqrt(fabs(Ka * Ka - Kb));
 161         resu = Ka + ((discrim > Ka)? discrim: (-discrim));
 162 #ifdef DEBUG
 163         if( DebugOn ) {
 164             printf("dux=%f, duy = %f, ka = %f, kb = %f\n",
 165                 Dux, Duy, Ka, Kb );
 166         }
 167 #endif /* DEBUG */
 168
 169     }
 170
 171     if( !IsQuadv ) {
 172         float denom;
 173         denom = (Dv1 - DotVector(&Na, &Ri));
 174         if( fabs(denom) < SingleTinyVal )
 175             resv = 2.0;
 176         else
 177             resv = (DotVector(&Nb, &Ri) - Dv0)/denom;
 178     } else {
 179         float Ka, Kb;
 180         float discrim;
 181
 182         Ka = Dvx + DotVector( &Qvx, &Ri);
 183         Kb = Dvy + DotVector( &Qvy, &Ri);
 184         discrim = sqrt(fabs( Ka * Ka - Kb));
 185         resv = Ka + ((discrim > Ka)? discrim: (-discrim));
 186 #ifdef DEBUG
 187         if( DebugOn ) {
 188             printf("dvx=%f, dvy = %f, ka = %f, kb = %f\n",
 189                 Dvx, Dvy, Ka, Kb );
 190         }
 191 #endif /* DEBUG */
 192     }
 193
 194 #ifdef DEBUG
 195     if( DebugOn ) {
 196         printf("(%f,%f) -> (%f, %f)\n", px, py, resu, resv );
 197     }
 198 #endif /* DEBUG */
 199
 200     if( resu > 1.0 || resu < 0.0 ) {
 201         resu = ( SweepCode & MP_XMAX)? 1.0: 0.0;
 202     }
 203     if( resv > 1.0 || resv < 0.0 ) {
 204         resv = ( SweepCode & MP_YMAX)? 1.0: 0.0;
 205     }
 206
 207     *outx = resu; *outy = resv;
 208 }
 209
 210
 211
 212 /*
 213  * here abides testcode
 214  */
 215 #ifdef TEST
 216 #include <stdio.h>
 217
 218 ReadScrnPair(set, a)
 219     char *set;
 220     ScrnPair *a;
 221 {
 222     int tx, ty;
 223     printf("enter screen pair %s\n",set);
 224     scanf("%d %d", &tx, &ty);
 225     a->x = tx; a->y = ty;
 226 }
 227
 228 ReadLimits(a)
 229 ScrnPair a[];
 230 {
 231     ReadScrnPair("a", &a[0]);
 232     ReadScrnPair("b", &a[1]);
 233     ReadScrnPair("c", &a[2]);
 234     ReadScrnPair("d", &a[3]);
 235 }
 236
 237 main() {
 238     float inx, iny;
 239     float outy, outx;
 240     ScrnPair pts[4];
 241
 242     ReadLimits(pts);
 243     CalcMapConsts(pts);
 244     while(!feof(stdin)) {
 245         printf("enter quadrilateral points\n");
 246         scanf("%f %f", &inx, &iny );
 247         MapXYRatio( inx, iny, &outx, &outy, 0);
 248         printf("p(%f,%f)->p(%f,%f)\n", inx, iny, outx, outy);
 249     }
 250 }
 251 #endif /* TEST */