src/snark/subreg.cpp

   1 /*
   2  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   3  *                                                                           *
   4  *                              S N A R K   1 4                              *
   5  *                                                                           *
   6  *                     A PICTURE RECONSTRUCTION PROGRAM                      *
   7  *                                                                           *
   8  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   9
  10  subreg.cpp,v 1.2 2008/08/25 16:11:09 jklukowska Exp
  11
  12  Subroutine to compute rough boundary (a bounding rectangle) of
  13  the specified region to reduce the search area that requires
  14  within region test.
  15  Specified region is either a rectangle or an ellipse.
  16  The bounding rectangle is the SMALLEST rectangle, with the same
  17  orientation as the picture, that bounds the nominal rectangle defined
  18  to have the same origin, axes and orientation as the specified region
  19  If this bounding rectangle (AN ESTIMATION) is partially outside the
  20  picture area, a WARNING message is output and the bounding rectangle
  21  is clipped to the picture boundary.
  22  This bounding rectangle is then passed to the DIST routine where a
  23  rigorous within boundary test reduced to this area is performed.
  24
  25  Parameters are passed through COMMON /region/
  26 */
  27
  28
  29 // parameter (MXNREG=40)
  30
  31 #include <cstdio>
  32 #include <cmath>
  33
  34 #include "blkdta.h"
  35 #include "geom.h"
  36 #include "region.h"
  37 #include "uiod.h"
  38 #include "consts.h"
  39
  40 #include "eval.h"
  41
  42 void Eval_class::subreg(INTEGER ObjJ)
  43 {
  44   REAL cx;
  45   REAL cy;
  46   REAL a;
  47   REAL b;
  48
  49   REAL theta;
  50   REAL cs;
  51   REAL ss;
  52
  53   INTEGER ixw;
  54   INTEGER iyw;
  55   INTEGER ic;
  56   INTEGER icx;
  57   INTEGER icy;
  58
  59   REAL tx1;
  60   REAL tx2;
  61   REAL ty1;
  62   REAL ty2;
  63
  64   REAL xw;
  65   REAL yw;
  66
  67   REAL x0;
  68   REAL y0;
  69
  70
  71   REAL xs;
  72   REAL ys;
  73
  74   static BOOLEAN flgout;
  75   //  static BOOLEAN objdon[MXNREG];  moved to eval.h - bug112 - swr - 6/10/05
  76
  77   BOOLEAN idone;
  78
  79   // Set up parameters
  80
  81   cx = Obj[ObjJ].cx;
  82   cy = Obj[ObjJ].cy;
  83   a = Obj[ObjJ].u;
  84   b = Obj[ObjJ].v;
  85   Region.aplus = a + (REAL) 0.0001;
  86   Region.bplus = b + (REAL) 0.0001;
  87   theta = Obj[ObjJ].ang * (REAL) Consts.pi/180.0;
  88   cs = (REAL) cos(theta);
  89   ss = (REAL) sin(theta);
  90
  91   Region.d1 = GeoPar.pixsiz * cs;
  92   Region.d2 = GeoPar.pixsiz * ss;
  93   Region.asq = SQR(a);
  94   Region.bsq = SQR(b);
  95
  96   // Estimate boundary of region
  97
  98   tx1 = (REAL) fabs(a * cs + b * ss);
  99   tx2 = (REAL) fabs(a * cs - b * ss);
 100   ty1 = (REAL) fabs(a * ss - b * cs);
 101   ty2 = (REAL) fabs(a * ss + b * cs);
 102   xw = MAX0(tx1, tx2);
 103   yw = MAX0(ty1, ty2);
 104   ixw = lround(xw / GeoPar.pixsiz);  // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 105   iyw = lround(yw / GeoPar.pixsiz);  // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 106   ic = GeoPar.nelem / 2;
 107   icx = ic + lround(cx / GeoPar.pixsiz);  // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 108   icy = ic - lround(cy / GeoPar.pixsiz);  // changed "(INTEGER)" to "lround" to match Fortran code. Lajos, Jan 20, 2005
 109
 110   Region.ix1 = icx - ixw;
 111   Region.ix2 = icx + ixw;
 112   Region.iy1 = icy - iyw;
 113   Region.iy2 = icy + iyw;
 114
 115   // Ensure estimated boundary is within picture region
 116
 117   idone = objdon[ObjJ];
 118   flgout = FALSE;
 119   if(Region.ix1 < 0) {
 120
 121     if((!flgout) && (idone == 0)) {
 122       fprintf(output, "\n         WARNING: Region %2i might not be entirely within picture", ObjJ);
 123       fprintf(output, "\n         Rough estimate of object boundary:");
 124       fprintf(output, "\n         ix1     ix2     iy1     iy2");
 125       fprintf(output, "\n         %3i     %3i     %3i     %3i", Region.ix1, Region.ix2, Region.iy1, Region.iy2);
 126     }
 127
 128     flgout = TRUE;
 129     Region.ix1 = 0;
 130   }
 131
 132   if(Region.ix2 > GeoPar.nelem - 1) {
 133     if((!flgout) && (idone == 0)) {
 134       fprintf(output, "\n         WARNING: Region %2i might not be entirely within picture", ObjJ);
 135       fprintf(output, "\n         Rough estimate of object boundary:");
 136       fprintf(output, "\n         ix1     ix2     iy1     iy2");
 137       fprintf(output, "\n         %3i     %3i     %3i     %3i", Region.ix1, Region.ix2, Region.iy1, Region.iy2);
 138     }
 139     flgout = TRUE;
 140     Region.ix2 = GeoPar.nelem - 1;
 141   }
 142
 143   if(Region.iy1 < 0) {
 144     if((!flgout) && (idone == 0)) {
 145       fprintf(output, "\n         WARNING: Region %2i might not be entirely within picture", ObjJ);
 146       fprintf(output, "\n         Rough estimate of object boundary:");
 147       fprintf(output, "\n         ix1     ix2     iy1     iy2");
 148       fprintf(output, "\n         %3i     %3i     %3i     %3i", Region.ix1, Region.ix2, Region.iy1, Region.iy2);
 149     }
 150     flgout = TRUE;
 151     Region.iy1 = 0;
 152   }
 153
 154   if(Region.iy2 > GeoPar.nelem - 1) {
 155     if((!flgout) && (idone == 0)) {
 156       fprintf(output, "\n         WARNING: Region %2i might not be entirely within picture", ObjJ);
 157       fprintf(output, "\n         Rough estimate of object boundary:");
 158       fprintf(output, "\n         ix1     ix2     iy1     iy2");
 159       fprintf(output, "\n         %3i     %3i     %3i     %3i", Region.ix1, Region.ix2, Region.iy1, Region.iy2);
 160     }
 161     flgout = TRUE;
 162     Region.iy2 = GeoPar.nelem - 1;
 163   }
 164
 165   objdon[ObjJ] = true;
 166
 167   x0 = (Region.ix1 - ic) * GeoPar.pixsiz;
 168   y0 = (ic - Region.iy1) * GeoPar.pixsiz;
 169   xs = x0 - cx;
 170   ys = y0 - cy;
 171   Region.xd0 = xs * cs + ys * ss;
 172   Region.yd0 = -xs * ss + ys * cs;
 173 }