src/snark/bwray.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/bwray.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  this subroutine is the WRAY for blobs
  11
  12  */
  13
  14 #include <cstdio>
  15 #include <cmath>
  16
  17 #include "consts.h"
  18 #include "blkdta.h"
  19 #include "blob.h"
  20 #include "posit.h"
  21 #include "anglst.h"
  22
  23 #include "uiod.h"
  24
  25 #include "raylen.h"
  26
  27 void Blob_class::bwray(INTEGER np, INTEGER nr, INTEGER* list, REAL* weight,
  28                 INTEGER* numb, REAL* snorm)
  29 {
  30         INTEGER a, b;
  31         INTEGER ac, bc, bc1, ac1;
  32         INTEGER af, al;
  33         INTEGER bf, bl;
  34         INTEGER laf, lal;
  35         INTEGER lbf, lbl;
  36         INTEGER h, m;
  37         INTEGER lmf, lml;
  38         INTEGER h1, m1;
  39         INTEGER mcov;
  40         REAL n;
  41         INTEGER nint;
  42         INTEGER sign;
  43         INTEGER M2, H2;
  44
  45         REAL ax, ay, axid, ayid;
  46         REAL mx, my;
  47         REAL theta, sinth, costh;
  48         REAL prfx, prlx, prfy, prly;
  49
  50         REAL x, y, yup, ylo;
  51         REAL ddeltain, ddeltaout;
  52         REAL dc, dc1;
  53         REAL ydelta; //xdelta not used
  54
  55         REAL slope, islope;
  56         REAL coef, icoef;
  57         REAL ord, cte, ctes;
  58         REAL fcosang;
  59         REAL samp_support;
  60         REAL multf;
  61
  62         *numb = 0;
  63
  64         *snorm = 0;
  65
  66         posit(np, nr, &ax, &ay, &mx, &my);
  67
  68         if (0 >= raylen(SOT_rect, ax, ay, mx, my, 0, 0, GeoPar.nelem * GeoPar.pixsiz / 2.0, GeoPar.nelem * GeoPar.pixsiz / 2.0, 1, 0))
  69         {
  70                 return; // bug 111, if ray(np, nr) miss all pixels in the image, return numb=0, by wei, 6/2005
  71         }
  72
  73         multf = Blob.support / Blob.delta;
  74
  75         samp_support = (REAL) samp / Blob.support;
  76
  77         laf = (INTEGER) (-0.5 * ((int) (Blob.M / 2) + (int) (Blob.H / 2)));
  78         lal =
  79                         (INTEGER) (0.5
  80                                         * ((int) ((Blob.M - 1) / 2) + (int) ((Blob.H - 1) / 2)));
  81
  82         lbf = -1 * (INTEGER) (Blob.H / 2);
  83         lbl = (INTEGER) ((Blob.H - 1) / 2);
  84
  85         lmf = -1 * (INTEGER) (Blob.M / 2);
  86         lml = (INTEGER) ((Blob.M - 1) / 2);
  87
  88         prfx = -0.25 * (Blob.M - 1.0);
  89         prlx = (-1) * prfx;
  90
  91         prfy = -0.25 * Consts.sqrt3 * (Blob.H - 1.0);
  92         prly = (-1) * prfy;
  93
  94         H2 = (INTEGER) (Blob.H / 2);
  95         M2 = (INTEGER) (Blob.M / 2);
  96         axid = ax / Blob.delta;
  97         ayid = ay / Blob.delta;
  98         sinth = my;
  99         costh = mx;
 100
 101         if (mx == 0)
 102                 theta = Consts.pi * 0.5;
 103         else
 104                 theta = (REAL) atan2(my, mx);
 105
 106         if (theta < 0)
 107                 theta += Consts.pi * 2; //use atan2 instead of acos for more acuracy,the range of atan2 is [-pi,pi], by wei 1/2005
 108
 109         //added 7/19, think of theta as always being betw 0 and pi
 110         if (theta > Consts.pi)
 111         {
 112                 theta = theta - Consts.pi;
 113                 sinth = -sinth;
 114                 costh = -costh; //never used anyway
 115                 my = -my;
 116                 mx = -mx;
 117         }
 118
 119         if (fabs(my) > Consts.zero)
 120         {  // non-horizontal lines
 121
 122                 islope = mx / my;
 123                 ord = (axid - islope * ayid);
 124                 coef = 0.5 * (Consts.sqrt3 * islope - 1);
 125
 126                 if (theta < Consts.pi / 6 || theta > 2.0 * Consts.pi / 3.0)
 127                 {  // type 1
 128
 129                         cte = multf / my; // for the "side rays"
 130                         slope = my / mx;
 131                         icoef = 1.0 / coef;
 132                         fcosang = (REAL) fabs(cos(theta + Consts.pi / 6.0));
 133
 134                         mcov = (INTEGER) ceil(multf / fcosang);
 135                         ddeltain = Blob.delta * fcosang;
 136                         ddeltaout = Blob.delta * sinth;
 137
 138                         if (theta < Consts.pi / 6)
 139                         {
 140                                 sign = 1;
 141                         }
 142                         else
 143                         {
 144                                 sign = -1;
 145                         }
 146
 147                         ctes = cte * sign;
 148
 149                         x = prfx;
 150                         yup = (x - ord + ctes) * slope;
 151                         ylo = (x - ord - ctes) * slope;
 152
 153                         if (sign == 1)
 154                         { // + slope
 155
 156                                 if (yup < prfy)
 157                                 {
 158                                         yup = prfy;
 159                                         x = ord - ctes + islope * yup;
 160                                 }
 161
 162                                 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
 163                         }
 164                         else
 165                         { //sign == -1
 166
 167                                 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
 168                                 if (ylo > prly)
 169                                 {
 170                                         ylo = prly;
 171                                         x = ord + ctes + islope * ylo;
 172                                         ac = (INTEGER) floor(x - ylo * Consts.isqrt3);
 173                                 }
 174                         }
 175
 176                         af = MAX0(laf, ac);
 177
 178                         x = prlx;
 179                         yup = (x - ord + ctes) * slope;
 180                         ylo = (x - ord - ctes) * slope;
 181
 182                         if (sign == 1)
 183                         { // + slope
 184
 185                                 if (ylo > prly)
 186                                 {
 187                                         ylo = prly;
 188                                         x = ord + ctes + islope * yup;
 189                                 }
 190                                 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
 191                         }
 192                         else
 193                         { //sign == -1
 194
 195                                 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
 196                                 if (yup < prfy)
 197                                 {
 198                                         yup = prfy;
 199                                         x = ord - ctes + islope * yup;
 200                                         ac = (INTEGER) ceil(x - yup * Consts.isqrt3);
 201                                 }
 202                         }
 203                         al = MIN0(lal, ac);
 204                         a = af;
 205                         bc = (INTEGER) (floor(icoef * (a - ord)));
 206
 207                         bf = MAX0(bc - mcov, lbf);
 208                         bl = MIN0(bc + mcov, lbl);
 209
 210                         x = a + 0.5 * bc;
 211                         y = 0.5 * Consts.sqrt3 * bc;
 212                         dc = Blob.delta * (REAL) fabs(-my * (axid - x) + mx * (ayid - y));
 213
 214                         bc1 = bc + 1;
 215                         dc1 = ddeltain - dc;
 216
 217                         // initial values
 218                         b = bc;
 219                         n = dc;
 220                         m = 2 * a + b;
 221                         h = b;
 222                         m1 = m + M2;
 223                         h1 = h + H2;
 224
 225                         //starting the loops--------------------------------
 226
 227                         while (a <= al)
 228                         {
 229                                 if (bf <= bl)
 230                                 {
 231                                         while (b >= bf)
 232                                         {
 233
 234                                                 if (m >= lmf && m <= lml && b <= lbl)
 235                                                 {
 236
 237                                                         if (n < Blob.support)
 238                                                         {
 239
 240                                                                 list[*numb] = h1 * Blob.M + m1;
 241                                                                 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
 242
 243                                                                 weight[*numb] = Blob.ltable[nint];
 244
 245                                                                 *snorm += weight[*numb] * weight[*numb];
 246                                                                 *numb = *numb + 1;
 247                                                         }
 248                                                 }
 249                                                 b--;
 250                                                 m--;
 251                                                 m1--;
 252                                                 h1--;
 253                                                 n += ddeltain;
 254                                         }
 255
 256                                         b = bc1;
 257                                         n = dc1;
 258                                         m = 2 * a + b;
 259                                         h = b;
 260                                         m1 = m + M2;
 261                                         h1 = h + H2;
 262
 263                                         while (b <= bl)
 264                                         {
 265
 266                                                 if (m >= lmf && m <= lml && b >= lbf)
 267                                                 {
 268                                                         if (n < Blob.support)
 269                                                         {
 270
 271                                                                 list[*numb] = h1 * Blob.M + m1;
 272                                                                 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
 273
 274                                                                 weight[*numb] = Blob.ltable[nint];
 275
 276                                                                 *snorm += weight[*numb] * weight[*numb];
 277                                                                 *numb = *numb + 1;
 278                                                         }
 279                                                 }
 280                                                 b++;
 281                                                 m++;
 282                                                 // h++;
 283                                                 m1++;
 284                                                 h1++;
 285                                                 n += ddeltain;
 286                                         }
 287                                 } //if bf<= bl
 288
 289                                 dc = dc + (REAL) sign * ddeltaout;
 290
 291                                 dc1 = ddeltain - dc;
 292
 293                                 if (dc > ddeltain || dc < 0.0)
 294                                 {
 295                                         bc = bc + sign;
 296                                         dc = (REAL) sign * (fabs(dc) - ddeltain);
 297                                         bc1 = bc + 1;
 298                                         dc1 = ddeltain - dc;
 299                                         bf = MAX0(bc - mcov, lbf);
 300                                         bl = MIN0(bc + mcov, lbl);
 301                                 }
 302
 303                                 a++;
 304                                 b = bc;
 305                                 n = dc;
 306                                 m = 2 * a + b;
 307                                 h = b;
 308                                 m1 = m + M2;
 309                                 h1 = h + H2;
 310                         } //while a
 311                 }   //type 1
 312
 313                 else
 314                 {   // type 2
 315                         fcosang = (REAL) fabs(cos(theta + Consts.pi / 6.0));
 316                         mcov = (INTEGER) ceil(multf / sinth);
 317                         // constants independent of the first grid point
 318                         bf = lbf;
 319                         bl = lbl;
 320                         ddeltain = Blob.delta * sinth;
 321                         ddeltaout = Blob.delta * fcosang;
 322
 323                         if (theta < Consts.pi / 3)
 324                         {
 325                                 sign = 1;
 326                         }
 327                         else
 328                         {
 329                                 sign = -1;
 330                         }
 331
 332                         //----------------------------------------------
 333                         // parameters of the first grid point
 334
 335                         b = bf;
 336
 337                         ac = (INTEGER) (floor(coef * b + ord));
 338
 339                         af = MAX0(ac - mcov, laf);
 340                         al = MIN0(ac + mcov, lal);
 341
 342                         x = ac + 0.5 * b;
 343                         y = 0.5 * Consts.sqrt3 * b;
 344                         dc = Blob.delta * (REAL) fabs(-my * (x - axid) + mx * (y - ayid));
 345
 346                         ac1 = ac + 1;
 347                         dc1 = ddeltain - dc;
 348
 349                         //init values
 350                         a = ac;
 351                         n = dc;
 352                         m = 2 * a + b;
 353                         h = b;
 354                         m1 = m + M2;
 355                         h1 = h + H2;
 356
 357                         // starting the loop .....................
 358
 359                         while (b <= bl)
 360                         {
 361                                 if (af <= al)
 362                                 {
 363                                         while (a >= af)
 364                                         {
 365
 366                                                 if (m >= lmf && m <= lml)
 367                                                 {
 368
 369                                                         if (n < Blob.support)
 370                                                         {
 371                                                                 list[*numb] = h1 * Blob.M + m1;
 372                                                                 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
 373
 374                                                                 weight[*numb] = Blob.ltable[nint];
 375
 376                                                                 *snorm += weight[*numb] * weight[*numb];
 377                                                                 *numb = *numb + 1;
 378                                                         }
 379                                                 }
 380                                                 a--;
 381                                                 m -= 2;
 382                                                 m1 -= 2;
 383                                                 n += ddeltain;
 384                                         }
 385
 386                                         a = ac1;
 387                                         n = dc1;
 388                                         m = 2 * a + b;
 389                                         h = b;
 390                                         m1 = m + M2;
 391                                         h1 = h + H2;
 392
 393                                         while (a <= al)
 394                                         {
 395                                                 if (m >= lmf && m <= lml)
 396                                                 {
 397                                                         if (n < Blob.support)
 398                                                         {
 399
 400                                                                 list[*numb] = h1 * Blob.M + m1;
 401
 402                                                                 nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
 403
 404                                                                 weight[*numb] = Blob.ltable[nint];
 405
 406                                                                 *snorm += weight[*numb] * weight[*numb];
 407                                                                 *numb = *numb + 1;
 408                                                         }
 409                                                 }
 410                                                 a++;
 411                                                 m += 2;
 412                                                 m1 += 2;
 413                                                 n += ddeltain;
 414                                         }
 415                                 } //if af <= al
 416
 417                                 dc = dc + (REAL) sign * ddeltaout;
 418
 419                                 dc1 = ddeltain - dc;
 420
 421                                 if (dc > ddeltain || dc < 0.0)
 422                                 {
 423                                         ac = ac + sign;
 424                                         dc = (REAL) sign * (fabs(dc) - ddeltain);
 425
 426                                         ac1 = ac + 1;
 427                                         dc1 = ddeltain - dc;
 428
 429                                         af = MAX0(ac - mcov, laf);
 430                                         al = MIN0(ac + mcov, lal);
 431                                 }
 432
 433                                 b++;
 434                                 a = ac;
 435                                 n = dc;
 436                                 m = 2 * a + b;
 437                                 h = b;
 438                                 m1 = m + M2;
 439                                 h1 = h + H2;
 440                         } //while b
 441                 }  // type 2
 442         }         // if non- horizontal line.
 443
 444         else
 445         {  //if horizontal line (my == 0)
 446
 447                 fcosang = cos(Consts.pi / 6.0);
 448
 449                 mcov = (INTEGER) ceil(multf / fcosang);
 450
 451                 // constants independent of the first grid point
 452                 ddeltain = Blob.delta * fcosang;
 453
 454                 //parameters of the first grid point
 455                 //find optimum af
 456
 457                 x = prfx;
 458                 yup = ayid + multf;
 459                 ylo = yup - 2 * multf;
 460
 461                 ac = (INTEGER) floor(x - yup * Consts.isqrt3);
 462
 463                 af = MAX0(laf, ac);
 464
 465                 //find optimum al
 466
 467                 x = prlx;
 468                 ac = (INTEGER) ceil(x - ylo * Consts.isqrt3);
 469                 al = MIN0(lal, ac);
 470
 471                 //done with af and al
 472
 473                 a = af;
 474                 bc = (INTEGER) (floor(2 * ayid * Consts.isqrt3));
 475                 bf = MAX0(bc - mcov, lbf);
 476                 bl = MIN0(bc + mcov, lbl);
 477
 478                 y = 0.5 * Consts.sqrt3 * bc;
 479                 ydelta = y * Blob.delta;
 480                 dc = (REAL) fabs(ydelta - ay);
 481                 bc1 = bc + 1;
 482                 dc1 = ddeltain - dc;
 483
 484                 // init values
 485                 b = bc;
 486                 n = dc;
 487                 //n = nc;
 488                 m = 2 * a + b;
 489                 h = b;
 490                 m1 = m + M2;
 491                 h1 = h + H2;
 492
 493                 while (a <= al)
 494                 {
 495                         if (bf <= bl)
 496                         {
 497                                 while (b >= bf)
 498                                 {
 499
 500                                         if (m >= lmf && m <= lml && b <= lbl)
 501                                         {
 502                                                 if (n < Blob.support)
 503                                                 {
 504
 505                                                         list[*numb] = h1 * Blob.M + m1;
 506                                                         nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
 507
 508                                                         weight[*numb] = Blob.ltable[nint];
 509                                                         *snorm += weight[*numb] * weight[*numb];
 510                                                         *numb = *numb + 1;
 511                                                 }
 512                                         }
 513                                         b--;
 514                                         m--;
 515
 516                                         m1--;
 517                                         h1--;
 518                                         n += ddeltain;
 519                                 }
 520
 521                                 b = bc1;
 522                                 n = dc1;
 523                                 m = 2 * a + b;
 524                                 h = b;
 525                                 m1 = m + M2;
 526                                 h1 = h + H2;
 527
 528                                 while (b <= bl)
 529                                 {
 530                                         if (m >= lmf && m <= lml && b >= lbf)
 531                                         {
 532                                                 if (n < Blob.support)
 533                                                 {
 534
 535                                                         list[*numb] = h1 * Blob.M + m1;
 536                                                         nint = (INTEGER) round(samp_support * n); // changed "rint" to "round". Lajos, Jan 25, 2005
 537
 538                                                         weight[*numb] = Blob.ltable[nint];
 539
 540                                                         *snorm += weight[*numb] * weight[*numb];
 541                                                         *numb = *numb + 1;
 542                                                 }
 543                                         }
 544                                         b++;
 545                                         m++;
 546
 547                                         m1++;
 548                                         h1++;
 549                                         n += ddeltain;
 550                                 }
 551                         } // if bf <= bl
 552
 553                         a++;
 554                         b = bc;
 555                         n = dc;
 556                         //n = nc;
 557                         m = 2 * a + b;
 558                         h = b;
 559                         m1 = m + M2;
 560                         h1 = h + H2;
 561
 562                 } //while a
 563         }  // fabs (my)==0
 564
 565         if (trace > 6)
 566         {
 567                 fprintf(output,
 568                                 "\n          bwray    np = %5i  nr = %5i  numb = %5i  snorm = %10.4f",
 569                                 np, nr, *numb, *snorm);
 570         }
 571
 572         if ((*numb > 0) && (trace > 9))
 573         {
 574                 for (int i = 0; i < *numb; i++)
 575                 {
 576
 577                         fprintf(output, "\n                  %5i     %10.4f", list[i],
 578                                         weight[i]);  //bug 139
 579                 }
 580         }
 581         return;
 582 }