ic.cc

来自「用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵」· CC 代码 · 共 601 行 · 第 1/2 页

          }
          else if ((xi < oldx) && (yi > oldy)) //down to left
          {
            intersected = Util::max(boundaries.H(oldx,yi),intersected); 
            intersected = Util::max(boundaries.H(xi,yi),intersected); 
            intersected = Util::max(boundaries.V(oldx,oldy),intersected); 
            intersected = Util::max(boundaries.V(oldx,yi),intersected); 
          }
          else if ((xi < oldx) && (yi < oldy)) //up to left
          {
            intersected = Util::max(boundaries.H(oldx,oldy),intersected); 
            intersected = Util::max(boundaries.H(xi,oldy),intersected); 
            intersected = Util::max(boundaries.V(oldx,oldy),intersected); 
            intersected = Util::max(boundaries.V(oldx,yi),intersected); 
          }
        }
        maxpb = Util::max(maxpb,intersected);
        oldx = xi;
        oldy = yi;

        // if the approximation is not good, then skip this point
        if (!good) { continue; }

        // if the accumulated pb is too high, then stop
        if (maxpb > thresh)
        {
          break;
        }

        // record this connection
        PointIC p;
        p.x = xi;
        p.y = yi;
        p.sim = 1.0f - maxpb;
        points(count) = p;
        count++;
      }

      // add our list of points to scanLines; we have to reverse
      // the order in octants 2,3,4,5
      switch (octant)
      {
        case 0:
        case 1:
        case 6:
        case 7:
          for (int i = 0; i < count; i++)
          {
            const int yind = points(i).y - y0 + wr;
            scanLines(yind,scanCount(yind)++) = points (i);
          }
          break;
        case 2:
        case 3:
        case 4:
        case 5:
          for (int i = count - 1; i >= 0; i--)
          {
            const int yind = points(i).y - y0 + wr;
            scanLines(yind,scanCount(yind)++) = points (i);
          }
          break;
        default:
            assert (0);
      }
  }


  //
  // given a pb image, a pixel (x0,y0), and a pb threshold, compute the intervening-contour 
  // weight to all pixels inside a given box of radius "wr" subject to the threshold "thresh".
  // pb is max-accumulated along bresenham lines.  the result is stored in scanline order as
  // a list of points and their pb value.
  //
  void interveningContour(const DualLattice& boundaries, const float thresh,
                          const int x0, const int y0, const int wr,
                          Util::Array1D<PointIC> &adj, int &count)
  {
      const int width = boundaries.width;
      const int height = boundaries.height;

      // make sure (x0,y0) is valid
      assert (x0 >= 0 && x0 < width);
      assert (y0 >= 0 && y0 < height);

      // make sure adj array is big enough
      adj.resize((2*wr+1)*(2*wr+1));

      // allocate space for lists of pixels; this operation is O(1)
      // since the space need not be initialized.
      Util::Array2D <PointIC>scanLines(2*wr+1,2*wr+1);

      // we need to keep track of the length of the scan lines
      Util::Array1D <int>scanCount(2*wr+1);
      scanCount.init(0);

      // scratch space for ic_walk() function
      Util::Array1D<PointIC> scratch(4*wr+2);

      // the rectangle of interest, a square with edge of length
      // 2*wr+1 clipped to the image dimensions
      const int rxa = Util::max(0,x0-wr);
      const int rya = Util::max(0,y0-wr);
      const int rxb = Util::min(x0+wr,width-1);
      const int ryb = Util::min(y0+wr,height-1);

      // walk around the boundary, collecting points in the scanline array
      // first walk around the rectangle boundary clockwise for theta = [pi,0]
      //std::cerr << "[" << x0 << "," << y0 << "]  ";
      //std::cerr << "(" << rxa << "," << rya << ")-(" << rxb << "," << ryb << ")" << std::endl;
      if (x0 > rxa) // left 
      {
        if ((y0 > 0) && (y0 < ryb))
        {
          ic_walk(boundaries, thresh, x0,y0, rxa,y0-1, rxa,y0, rxa,y0+1, wr, scratch,scanCount,scanLines);
        }
        for (int y = y0-1; y > rya; y--)
        {
          ic_walk(boundaries, thresh, x0,y0, rxa,y-1, rxa,y, rxa,y+1, wr, scratch,scanCount,scanLines);
        }
      }

      if (x0 > rxa+1 || y0 > rya+1 || ((x0 > rxa) && (y0 > rya)) ) // top-left
      {
        ic_walk(boundaries, thresh, x0,y0, rxa,rya+1, rxa,rya, rxa+1,rya, wr, scratch, scanCount, scanLines);
      }
      if ( ((x0 == rxa) && (y0 == rya+1)) || ((x0 == rxa+1) && (y0 == rya)) )
      {
        PointIC pnt;
        pnt.x = rxa;
        pnt.y = rya;
        pnt.sim = 1.0f;
        const int yind = pnt.y - y0 + wr;
        scanLines(yind,scanCount(yind)++) = pnt;
      }

      if (y0 > rya) // top
      {
        for (int x = rxa+1; x < rxb; x++)
        {
          ic_walk(boundaries, thresh, x0,y0, x-1,rya, x,rya, x+1,rya, wr, scratch, scanCount, scanLines);
        }
      }

      if ((y0 > rya+1) || (x0 < rxb-1) || ((y0 > rya) && (x0 < rxb)) ) // top-right
      {
        ic_walk(boundaries, thresh, x0,y0, rxb-1,rya, rxb,rya, rxb,rya+1, wr, scratch, scanCount, scanLines);
      }
      if ( ((x0 == rxb-1) && (y0 == rya)) || ((x0 == rxb) && (y0 == rya+1)) )
      {
        PointIC pnt;
        pnt.x = rxb;
        pnt.y = rya;
        pnt.sim = 1.0f;
        const int yind = pnt.y - y0 + wr;
        scanLines(yind,scanCount(yind)++) = pnt;
      }


      if (x0 < rxb) // right
      {
        for (int y = rya+1; y < y0; y++)
        {
          ic_walk(boundaries, thresh, x0,y0, rxb,y-1, rxb,y, rxb,y+1, wr, scratch, scanCount, scanLines);
        }
      }

      // now counterclockwise for theta = (pi,0)
      if (x0 > rxa) // left
      {
        for (int y = y0+1; y < ryb; y++)
        {
          ic_walk(boundaries, thresh, x0,y0, rxa,y-1, rxa,y, rxa,y+1, wr, scratch, scanCount, scanLines);
        }
      }

      if ((x0 > rxa+1) || (y0 < ryb-1) || ((x0 > rxa) && (y0 < ryb))) // bottom-left
      {
        ic_walk(boundaries, thresh, x0,y0, rxa,ryb-1, rxa,ryb, rxa+1,ryb, wr, scratch, scanCount, scanLines);
      }
      if ( ((x0 == rxa) && (y0 == ryb-1)) || ((x0 == rxa+1) && (y0 == ryb)) )
      {
        PointIC pnt;
        pnt.x = rxa;
        pnt.y = ryb;
        pnt.sim = 1.0f;
        const int yind = pnt.y - y0 + wr;
        scanLines(yind,scanCount(yind)++) = pnt;
      }
      if (y0 < ryb) // bottom
      {
        for (int x = rxa+1; x < rxb; x++)
        {
          ic_walk(boundaries, thresh, x0,y0, x-1,ryb, x,ryb, x+1,ryb, wr, scratch, scanCount, scanLines);
        }
      }

      if ((y0 < ryb-1) || (x0 < rxb-1) || ((y0 < ryb) && (x0 < rxb))) // bottom-right
      {
        ic_walk(boundaries, thresh, x0,y0, rxb-1,ryb, rxb,ryb, rxb,ryb-1, wr, scratch, scanCount, scanLines);
      }
      if ( ((x0 == rxb-1) && (y0 == ryb)) || ((x0 == rxb) && (y0 == ryb-1)) )
      {
        PointIC pnt;
        pnt.x = rxb;
        pnt.y = ryb;
        pnt.sim = 1.0f;
        const int yind = pnt.y - y0 + wr;
        scanLines(yind,scanCount(yind)++) = pnt;
      }

      if (x0 < rxb) // right
      {
        for (int y = ryb-1; y > y0; y--)
        {
          ic_walk(boundaries, thresh, x0,y0, rxb,y-1, rxb,y, rxb,y+1, wr, scratch, scanCount, scanLines);
        }
        if ((y0 > 0) && (y0 < ryb))
        {
          ic_walk(boundaries, thresh, x0,y0, rxb,y0-1, rxb,y0, rxb,y0+1, wr, scratch, scanCount, scanLines);
        }
      }


      for (int y = 0; y < 2*wr+1; y++)
      {
        int len = scanCount(y);
        assert (len >= 0 && len <= 2*wr+1);

        if (y + y0 - wr < 0)
        {
          assert(len == 0);    
        }

        // check that pixels are in the right row
        for (int i = 0; i < len; i++)
        {
          assert(scanLines(y,i).y == y+y0-wr);
        }

        // check that pixels in each row are in increasing order
        for (int i = 0; i < len-1; i++)
        {
          assert(scanLines(y,i).x < scanLines(y,i+1).x);
        }
      }

      // construct the adjacency list
      count = 0;
      for (int y = 0; y < 2*wr+1; y++)
      {
        int len = scanCount(y);
        for (int i = 0; i < len; i++)
        {
          adj(count++) = scanLines(y,i);
        }
      }
  }


  // compute the max over lattice energies on a straightline 
  // path connecting p1 and p2
/*
  void interveningContour (const DualLattice& boundaries, const int x1, const int y1, 
                            const int x2, const int y2, float& icsim)
  {
      float maxpb = 0;

      const int width = boundaries.pb.size(0);
      const int height = boundaries.pb.size(1);
      const int dx = x2 - x1;
      const int dy = y2 - y1;
      const int steps = Util::max (abs (dx), abs (dy));
      if (steps == 0) { return; }

      const float xincr = (float) dx / (float) steps;
      const float yincr = (float) dy / (float) steps;

      float x = x1;
      float y = y1;
      float olddist = boundaries.dist(x1,y1);
      for (int k = 0; k < steps; k++)
      {
        x += xincr;
        y += yincr;
        const int xi = (int) rint (x);
        const int yi = (int) rint (y);
        float newdist = boundaries.dist(xi,yi);
        if ( ((olddist >= 0) && (newdist < 0)) || 
              ((olddist <= 0) && (newdist > 0)) )
        {
          maxpb = Util::max(maxpb, boundaries.pb(xi,yi));
        }
        olddist = newdist;
      }
      icsim = 1-maxpb;
  }
*/
} //namespace Group