src/snark/raylen.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/raylen.cpp $
   5  $LastChangedRevision: 91 $
   6  $Date: 2014-07-02 17:34:27 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  RAYLEN RETURNS THE LENGTH OF THE RAY IN AN ELEMENTAL OBJECT
  11
  12  SEE SNARK05 MANUAL FOR A FULL DESCRIPTION OF THE FUNCTIONING OF
  13  THE ROUTINE
  14
  15  TYPE = 1 FOR ELLIPSE/CIRCLE
  16  TYPE = 2 FOR RECTANGLE/SQUARE
  17  TYPE = 3 FOR ISOS. TRIANGLE
  18  TYPE = 4 FOR SEGMENT OF CIRCLE
  19  TYPE = 5 FOR SECTOR OF CIRCLE
  20  AX,AY AND MX,MY DEFINE THE RAY:
  21  (AX,AY) COORDINATES OF ANY POINT ON THE RAY;
  22  (MX,MY) ARE THE X- AND Y- DIRECTION COSINES OF THE RAY
  23  IF THE RAY MAKES AN ANGLE THETA WITH X-AXIS THEN MX = COS(THETA)
  24  AND MY = SIN(THETA)
  25
  26  (CX,CY) ARE THE COORDINATES OF THE CENTER OF THE OBJECT
  27
  28  (U,V) ARE THE HALF- AXIS LENGTHS OF THE OBJECT
  29
  30  (PX,PY) ARE THE COSINE AND SINE OF THE ANGLE OF ROTATION OF THE
  31  OBJECT ( THE ANGLE U MAKES WITH THE X-AXIS)
  32  */
  33
  34 #include <cstdlib>
  35 #include <cstdio>
  36 #include <cmath>
  37
  38 #include "blkdta.h"
  39 #include "raylen.h"
  40 #include "consts.h"
  41
  42 REAL raylen(SNARK_Object_type type, REAL ax, REAL ay, REAL mx, REAL my, REAL cx,
  43                 REAL cy, REAL u, REAL v, REAL px, REAL py)
  44 {
  45         REAL raylen_temp;
  46         REAL xo, yo;
  47
  48         REAL xnew, ynew, unew, vnew, xdc, ydc;
  49
  50         REAL diam;
  51
  52         REAL temp, slope;
  53         REAL ycept;
  54         REAL xcept;
  55
  56         REAL x1, x2;
  57
  58         REAL a, b, c, d, bb;
  59
  60         BOOLEAN rotate;
  61
  62         REAL vtou, vsu;
  63
  64         raylen_temp = 0.0;
  65
  66         // SHIFT THE ORIGIN TO THE CENTER OF THE OBJECT
  67
  68         xo = ax - cx;
  69         yo = ay - cy;
  70
  71         // ROTATE ABOUT ORIGIN BY AN ANGLE = (U,X-AXIS) SO THAT U AND V ARE
  72         // ALONG THE NEW X AND Y AXES RESPECTIVELY
  73
  74         xnew = xo * px + yo * py;
  75         ynew = -xo * py + yo * px;
  76         unew = u;
  77         vnew = v;
  78         xdc = mx * px + my * py;
  79         ydc = -mx * py + my * px;
  80
  81
  82         // FIND THE CIRCLE COMPLETELY ENCLOSING OBJECT
  83
  84         diam = (REAL) 2.1 * MAX0(unew, vnew);
  85
  86         // CHECK IF RAY HAS SLOPE .LE. 1.0  FOR COMPUTATIONAL PURPOSES
  87         // REMEMBER THIS FOR TRIA, SEGM, AND SECT SINCE THEY DO NOT
  88         // HAVE TWO-FOLD SYMMETRY
  89
  90         rotate = FALSE;
  91         if (!((REAL) fabs(xdc) >= (REAL) fabs(ydc)))
  92         {
  93
  94                 // ROTATE ABOUT ORIGIN BY 90 DEG TO INSURE THAT SLOPE IS LE 1.0 IN
  95                 // MAGNITUDE
  96
  97                 rotate = TRUE;
  98                 temp = xnew;
  99                 xnew = ynew;
 100                 ynew = -temp;
 101
 102                 temp = xdc;
 103                 xdc = ydc;
 104                 ydc = -temp;
 105
 106                 temp = unew;
 107                 unew = vnew;
 108                 vnew = temp;
 109         }
 110
 111         // SEE IF RAY IS ANYWHERE NEAR OBJECT CIRCLE
 112
 113         slope = ydc / xdc;
 114         if ((REAL) fabs(slope) <= (REAL) 0.000001)
 115         {
 116                 slope = 0;
 117         }
 118
 119         ycept = ynew - slope * xnew;
 120
 121         if ((REAL) fabs(ycept) < (REAL) 0.000001)
 122         {
 123                 ycept = 0.0;
 124         }
 125
 126         xcept = 0.0;
 127
 128         if ((REAL) fabs(slope) > Consts.zero)
 129         {
 130                 xcept = xnew - ynew / slope;
 131         }
 132
 133         if (((REAL) fabs(xcept) > diam) && ((REAL) fabs(ycept) > diam))
 134         {
 135                 return raylen_temp;
 136         }
 137
 138         // if RECTANGLE or CIRCLE/ELLIPSE
 139
 140         if (type == SOT_rect)
 141         {
 142                 // CHECK IF RAY THROUGH RECTANGLE/SQUARE IS NEARLY PARALLEL TO X AXIS
 143                 if (!((REAL) fabs(ydc) > (REAL) 0.0001))
 144                 {
 145                         if ((REAL) fabs(ynew) <= vnew)
 146                                 raylen_temp = (REAL) 2.0 * unew;
 147
 148                         return raylen_temp;
 149                 }
 150                 // GENERAL LINE INSIDE A RECTANGLE
 151
 152 #ifdef FFCOMPARE
 153
 154                 // there is still difference here between fortran
 155                 raylen_temp = (REAL) fabs(
 156                                 (
 157                                                 MAX0(-unew, MIN0(unew, ((REAL) (vnew - ynew) / slope) + xnew))
 158                                                 - MAX0(-unew, MIN0(unew, ((REAL) -(vnew + ynew) / slope) + xnew))
 159                                 )
 160                                 / xdc
 161                 );
 162 #else
 163
 164                 raylen_temp =
 165                                 (REAL) fabs(
 166                                                 (
 167                                                 MAX0(-unew, MIN0(unew, (vnew - ynew) / slope + xnew))
 168                                                                 - MAX0(-unew,
 169                                                                                 MIN0(unew, -(vnew + ynew) / slope + xnew)))
 170                                                                 / xdc);
 171 #endif
 172
 173                 return raylen_temp;
 174         }
 175
 176         // GENERAL LINE INSIDE ELLIPSE
 177         if (type == SOT_elip)
 178         {
 179                 vsu = vnew * vnew + (unew * slope) * (unew * slope);
 180                 d = vsu - ycept * ycept;
 181
 182                 // IF D POSITIVE, RAY INTERSECTS THE CIRCLE/ELLIPSE
 183                 if (d > Consts.zero)
 184                 {
 185                         raylen_temp = (REAL) 2.0 * unew * vnew * (REAL) sqrt(d)
 186                                         / (REAL) fabs(vsu * xdc);
 187 #ifdef FFCOMPARE
 188                         REAL zzz = (REAL) sqrt(d);
 189                         REAL xxx = (REAL) 2.0 * unew * vnew * zzz;
 190                         REAL yyy = (REAL) fabs(vsu * xdc);
 191                         REAL aaa = (REAL) xxx / yyy;
 192                         raylen_temp = (REAL) xxx / yyy;
 193 #endif
 194                 }
 195                 return raylen_temp;
 196         }
 197
 198         // TRIANGLE, SEGMENT, SECTOR
 199
 200         if (!(type == SOT_segm))
 201         {
 202
 203                 // TRIANGLE AND TRIANGULAR PORTION OF SECTOR
 204
 205                 // SPECIAL CASE FOR TRIANGLE AND SECTOR
 206                 // CHECK IF RAY IS NEARLY PARALLEL TO ONE OF THE SIDES
 207
 208                 vtou = vnew / unew;
 209                 if (!((REAL) fabs(vtou - (REAL) fabs(slope)) >= (REAL) 0.0001))
 210                 {
 211                         if (!(rotate))
 212                         {
 213                                 if (!((REAL) fabs(xcept) > unew))
 214                                 {
 215                                         if (slope >= 0.0)
 216                                                 raylen_temp = (REAL) 0.5
 217                                                                 * (REAL) fabs((unew - xcept) / xdc);
 218                                         if (slope < 0.0)
 219                                                 raylen_temp = (REAL) 0.5
 220                                                                 * (REAL) fabs((unew + xcept) / xdc);
 221                                 }
 222                         }
 223                         else
 224                         {
 225                                 if (!((REAL) fabs(ycept) > vnew))
 226                                 {
 227                                         if (slope >= 0.0)
 228                                                 raylen_temp = (REAL) 0.5
 229                                                                 * (REAL) fabs((vnew - ycept) / ydc);
 230                                         if (slope < 0.0)
 231                                                 raylen_temp = (REAL) 0.5
 232                                                                 * (REAL) fabs((vnew + ycept) / ydc);
 233                                 }
 234                         }
 235                 }
 236                 else
 237                 {
 238                         // TRIANGLE  GENERAL CASE
 239
 240                         if (!(rotate))
 241                         {
 242                                 x1 = (vnew - ycept) / (slope + vtou);
 243                                 if ((x1 <= -Consts.zero) || (x1 > unew))  //bug 273, jklukowska
 244                                         x1 = xcept;
 245                                 x2 = (vnew - ycept) / (slope - vtou);
 246                                 if ((x2 < -unew) || (x2 >= Consts.zero))  //bug 273, jklukowska
 247                                         x2 = xcept;
 248
 249                                 // RAY THROUGH VERTEX OF ISOSCELES TRIANGLE
 250                                 if (((REAL) fabs(x1) < (REAL) 0.00001)
 251                                                 && ((REAL) fabs(x2) < (REAL) 0.00001)
 252                                                 && ((REAL) fabs(xcept) <= unew))
 253                                         x1 = xcept;
 254                                 raylen_temp = (REAL) fabs((x1 - x2) / xdc);
 255                         }
 256                         else
 257                         {
 258
 259                                 x1 = (vnew - ycept) / (slope + vtou);
 260                                 if ((x1 < 0.0) || (x1 >= unew + 0.00001))
 261                                         x1 = 0.0;
 262                                 x2 = (-vnew - ycept) / (slope - vtou);
 263                                 if ((x2 < 0.0) || (x2 >= unew + 0.00001))
 264                                         x2 = 0.0;
 265                                 // RAY THROUGH VERTEX OF ISOSCELES TRIANGLE
 266
 267                                 if (((REAL) fabs(x1 - unew) < (REAL) 0.000001)
 268                                                 && ((REAL) fabs(x2 - unew) < (REAL) 0.000001)
 269                                                 && ((REAL) fabs(ycept) <= vnew))
 270                                         x2 = 0.0;
 271                                 raylen_temp = (REAL) fabs((x1 - x2) / xdc);
 272                         }
 273
 274                 }
 275
 276                 if (type == SOT_tria)
 277                 {
 278                         return raylen_temp;
 279                 }
 280         }
 281
 282         // SEGMENT OR THE SEGMENT PORTION OF THE SECTOR OF A CIRCLE
 283
 284         if (!rotate)
 285         {
 286                 if (((REAL) fabs(xcept) > unew) && (ycept > 0.0))
 287                 {
 288                         return raylen_temp;
 289                 }
 290
 291                 a = (REAL) 1.0 + slope * slope;
 292                 bb = ycept - vnew;
 293                 b = slope * bb;
 294                 c = bb * bb - (unew * unew + vnew * vnew);
 295                 d = b * b - a * c;
 296                 if (d <= 0.0)
 297                 {
 298                         return raylen_temp;
 299                 }
 300
 301                 // THE RAY DOES INTERSECT THE CIRCLE
 302                 d = (REAL) sqrt(d);
 303                 if (b >= 0.0)
 304                         x1 = -(b + d) / a;
 305                 if (b < 0.0)
 306                         x1 = (-b + d) / a;
 307                 x2 = c / (a * x1);
 308                 if (x1 * slope + ycept > 0.0)
 309                         goto L111;
 310                 if (x2 * slope + ycept > 0.0)
 311                         goto L112;
 312
 313                 // X COMPONENT OF INTERSECTION LENGTH IS DIFFERENCE OF ROOTS OF QUADR
 314                 if (raylen_temp > 0.0)
 315                 {
 316                         return raylen_temp;
 317                 }
 318                 raylen_temp += (REAL) 2.0 * d / (REAL) fabs(a * xdc);
 319                 return raylen_temp;
 320
 321                 // RAY INTERSECTS LINEAR BOUNDARY OF SEGMENT
 322                 L111: x1 = xcept;
 323                 if (x2 * slope + ycept <= 0.0)
 324                         goto L113;
 325
 326                 L112: x2 = xcept;
 327
 328                 L113: raylen_temp += (REAL) fabs((x1 - x2) / xdc);
 329                 return raylen_temp;
 330         }
 331
 332         if (((REAL) fabs(ycept) > vnew) && (xcept > 0.0))
 333         {
 334                 return raylen_temp;
 335         }
 336
 337         a = (REAL) 1.0 + slope * slope;
 338         b = slope * ycept - unew;
 339         c = ycept * ycept - vnew * vnew;
 340         d = b * b - a * c;
 341
 342         if (d <= 0.0)
 343         {
 344                 return raylen_temp;
 345         }
 346
 347         d = (REAL) sqrt(d);
 348         if (b >= 0.0)
 349                 x1 = -(b + d) / a;
 350         if (b < 0.0)
 351                 x1 = (-b + d) / a;
 352         x2 = c / (a * x1);
 353         if (x1 > 0.0)
 354                 goto L121;
 355         if (x2 > 0.0)
 356                 goto L122;
 357
 358         // X COMPONENT OF INTERSECTION LENGTH IS DIFFERENCE OF ROOTS OF QUADR
 359         if (raylen_temp > 0.0)
 360         {
 361                 return raylen_temp;
 362         }
 363         raylen_temp += (REAL) 2.0 * d / (REAL) fabs(a * xdc);
 364         return raylen_temp;
 365
 366         // RAY INTERSECTS LINEAR BOUNDARY OF SEGMENT
 367         L121: x1 = 0.0;
 368         if (x2 <= 0.0)
 369                 goto L123;
 370
 371         L122: x2 = 0.0;
 372
 373         L123: raylen_temp += (REAL) fabs((x1 - x2) / xdc);
 374         return raylen_temp;
 375 }