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