segmentation.cc

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

    //apply renumbering
    for (int x = 0; x < _map.size(0); x++)
    {
      for (int y = 0; y < _map.size(1); y++)
      {
        _map(x,y) = renum(_map(x,y));
        assert(_map(x,y) >= 0);
        assert(_map(x,y) < _numSegments);
      }
    }
  }


  //
  // break a segmentation into connected components so that
  // every segment label refers to exactly one connected 
  // component.
  //
  void
  Segmentation::fragment (float radius)
  {
    int nextId = _numSegments;
    for (int x = 0; x < _map.size(0); x++) 
    {
      for (int y = 0; y < _map.size(1); y++) 
      {
        int oldId = _map(x,y);
        if (oldId < _numSegments)
        {
          _replace(x,y,oldId,nextId,radius);
          nextId++;
        }
      }
    }
    _renumber();
  }

  //
  // merge segments containing fewer than dustthresh pixels
  // into neighboring regions
  //
  void
  Segmentation::mergeDust(int dustthresh)
  {
    fragment(1.0f);

    Array1D<int> sizes;
    segSizes(sizes);

    //if segment s is dustlike, merge with a neighbor
    Util::Message::debug(String("merging dust starting with %d segments",_numSegments),1);

    Array1D<int> nbrcnt(_numSegments);
    for (int s = 0; s < _numSegments; s++)
    {
      if (sizes(s) < dustthresh)
      {
        nbrcnt.init(0);

        //find how many nbring pixels this pixel has from each segment
        for (int x = 1; x < _map.size(0)-1; x++)
        {
          for (int y = 1; y < _map.size(1)-1; y++)
          {
            //if we are a neigbhor of s, add 1 to our nbrcnt
            if ( (_map(x  ,y+1) == s) || (_map(x+1,y+1) == s) ||
                 (_map(x+1,y  ) == s) || (_map(x+1,y-1) == s) ||
                 (_map(x  ,y-1) == s) || (_map(x-1,y-1) == s) ||
                 (_map(x-1,y  ) == s) || (_map(x-1,y+1) == s) )
            {
              nbrcnt(_map(x,y))++;   
            }
          }
        }

        //pick the most commonly occuring neighbor
        int bigestnbr = 0;
        for (int ss = 0; ss < _numSegments; ss++)
        {
          if (nbrcnt(ss) > nbrcnt(bigestnbr))
          {
            bigestnbr = ss;
          }
        }
        
        //and merge with that nbr
        for (int x = 0; x < _map.size(0); x++)
        {
          for (int y = 0; y < _map.size(1); y++)
          {
            if (_map(x,y) == s)
            {
              _map(x,y) = bigestnbr;
            }
          }
        }
      }
    }

    //lastly, renumber
    _renumber();
    Util::Message::debug(String("merged down to %d segments",_numSegments),1);
  }


  ////////////////////////////////////////////////////////////////////////////////
  //
  // segment -> boundary conversions
  //


  void Segmentation::computeBoundaryMap(Util::Array2D<bool>& boundaryMap) const
  {
    const int w = getWidth();
    const int h = getHeight();
    boundaryMap.resize(w,h);
    for (int x = 0; x < w; x++) {
        for (int y = 0; y < h; y++) {
            const int p = _map(x,y);
            if (x < w-1 && p != _map(x+1,y) ||
                y < h-1 && p != _map(x,y+1) ||
                x < w-1 && y < h-1 && p != _map(x+1,y+1))
            {
                boundaryMap(x,y) = true;
            }
            else
            {
                boundaryMap(x,y) = false;
            }
        }
    }
  }

  void Segmentation::computeBoundaryMapHalf(Util::Array2D<bool>& boundaryMap) const
  {
    const int w = getWidth();
    const int h = getHeight();
    boundaryMap.resize(w/2,h/2);
    for (int x = 0; x/2 < w/2; x+=2) {
        for (int y = 0; y/2 < h/2; y+=2) {
            const int p = _map(x,y);
            if (x < w-1 && p != _map(x+1,y) ||
                y < h-1 && p != _map(x,y+1) ||
                x < w-1 && y < h-1 && p != _map(x+1,y+1))
            {
                boundaryMap(x/2,y/2) = true;
            }
            else
            {
                boundaryMap(x/2,y/2) = false;
            }
        }
    }
  }

  class Point {
    public:
      Point () { x = y = -1; }
      Point (int x, int y) { this->x = x; this->y = y; }
      int x, y;
  };

  class DistNode 
  {
    public:
      DistNode () { dSqr = FLT_MAX; epoch = -1; }
      DistNode (int x, int y, float dSqr, int epoch) 
          : p(x,y) { this->dSqr = dSqr; this->epoch = epoch; }
      bool valid () { return dSqr != FLT_MAX; }
      Point p;
      float dSqr;
      int epoch;
  };

  // compute the distance from every point in the image to the boundary.
  void Segmentation::computeDistanceMap (Image& distanceMap) const
  {
      const int width = getWidth();
      const int height = getHeight();
      distanceMap.resize(width,height);

      Util::Array2D<bool> bmap;
      computeBoundaryMap(bmap);

      Util::Array1D<Point>* frontier = 
          new Util::Array1D<Point> (width*height);
      Util::Array1D<Point>* newFrontier = 
          new Util::Array1D<Point> (width*height);

      int epoch = 0;
      // initialize the frontier set
      Util::Array2D<DistNode> dmap (width,height);
      int frontierSize = 0;
      for (int x = 0; x < width; x++) 
      {
        for (int y = 0; y < height; y++) 
        {
            if (bmap(x,y) == true) 
            {
                (*frontier)(frontierSize++) = Point(x,y);
                dmap(x,y) = DistNode(x,y,0,epoch);
            }
        }
      }

      // enlarge frontier until all pixels are covered
      while (frontierSize > 0) {
          //std::cerr << "frontier size = " << frontierSize << std::endl;
          epoch++;
          int newFrontierSize = 0;
          for (int i = 0; i < frontierSize; i++) {
              const int x = (*frontier)(i).x;
              const int y = (*frontier)(i).y;
              assert (dmap(x,y).valid());
              assert (epoch==1 || bmap(x,y) != true);
              assert (dmap(x,y).epoch >= epoch-1);
              const int hx = dmap(x,y).p.x;
              const int hy = dmap(x,y).p.y;
              assert (dmap(hx,hy).valid());
              assert (bmap(hx,hy) == true);
              for (int v = -1; v <= 1; v++) {
                  for (int u = -1; u <= 1; u++) {
                      const int ix = x + u;
                      const int iy = y + v;
                      if (ix < 0 || ix >= width) { continue; }
                      if (iy < 0 || iy >= height) { continue; }
                      if (u == 0 && v == 0) { continue; }
                      const int dx = hx - ix;
                      const int dy = hy - iy;
                      const float dSqr = dx*dx + dy*dy;
                      if (dSqr < dmap(ix,iy).dSqr) {
                          if (dmap(ix,iy).epoch < epoch) {
                              (*newFrontier)(newFrontierSize++) = 
                                  Point(ix,iy);
                          }
                          dmap(ix,iy) = DistNode(hx,hy,dSqr,epoch);
                      }
                  }
              }
          }
          frontierSize = newFrontierSize;
          Util::swap(frontier,newFrontier);
      }

      // record distance map and homes
      for (int x = 0; x < width; x++) {
          for (int y = 0; y < height; y++) {
              const float d = sqrt(dmap(x,y).dSqr);
              distanceMap(x,y) = d;
          }
      }

      // clean up
      delete frontier;
      delete newFrontier;
  }

  // compute orientation of maximum variance of a list of points, [0,PI)
  float
  primaryOrient (const Array2D<float> points, const int numPoints) 
  {
      assert (numPoints > 0);
      assert (points.size(1) == 2);

      // compute means
      double meanx=0, meany=0;
      for (int i = 0; i < numPoints; i++) 
      {
          const double x = points(i,0);
          const double y = points(i,1);
          meanx += x;
          meany += y;
      }    
      meanx /= numPoints;
      meany /= numPoints;
     
      // compute (unnormalized) covariance matrix
      //	a b
      //  c d
      double a=0, b=0, c=0, d=0;
      for (int i = 0; i < numPoints; i++) 
      {
          const double x = points(i,0);
          const double y = points(i,1);
          a += (x - meanx) * (x - meanx);
          b += (x - meanx) * (y - meany);
          c += (y - meany) * (x - meanx);
          d += (y - meany) * (y - meany);
      }

      // compute eigenvalues
      const double k = sqrt (a*a + 4*b*c - 2*a*d + d*d);
      const double e0 = (a + d - k) / 2;
      const double e1 = (a + d + k) / 2;

      // compute eigenvectors
      const double v0[2] = { -b , a - e0 };
      const double v1[2] = { -b , a - e1 };

      // validate eigenpairs
      const double res1 = a*v0[0] + b*v0[1] - e0*v0[0];
      const double res2 = c*v0[0] + d*v0[1] - e0*v0[1];
      const double res3 = a*v1[0] + b*v1[1] - e1*v1[0];
      const double res4 = c*v1[0] + d*v1[1] - e1*v1[1];
      assert (fabs(res1) < 1e-9);
      assert (fabs(res2) < 1e-9);
      assert (fabs(res3) < 1e-9);
      assert (fabs(res4) < 1e-9);

      // get eigenvector corresponding to largest eigenvalue
      double v[2];
      if (fabs(e0) > fabs(e1)) {
          v[0] = v0[0];
          v[1] = v0[1];
      } else {
          v[0] = v1[0];
          v[1] = v1[1];
      }

      // compute direction of eigenvector
      float theta = atan2(v[1],v[0]);
      assert(finite(theta));
      if (theta < 0) { theta += M_PI; }
      if (theta < 0) { theta = 0; }
      if (theta >= M_PI) { theta = 0; }
      assert (finite(theta));
      assert (theta >= 0 && theta < M_PI);
      return theta;
  }

  // given a boundary map (bmap), estimate the boundary orientation
  // at each boundary point using a window of the given radius
  void Segmentation::computeOrientationMap (const float radius, Image& orientationMap) const
  {
      const int width = getWidth();
      const int height = getHeight();
      const int winrad = (int) ceil (radius);
      Util::Array2D<float> points(width*height,2);

      Util::Array2D<bool> bmap;
      computeBoundaryMap(bmap);

      orientationMap.resize(width,height);
      for (int x = 0; x < width; x++) {
          for (int y = 0; y < height; y++) {
              if (bmap(x,y) != true) { continue; }
              int numPoints = 0;
              for (int v = -winrad; v <= winrad; v++) {
                  for (int u = -winrad; u <= winrad; u++) {
                      int ix = x + u;
                      int iy = y + v;
                      if (ix < 0 || ix >= width) { continue; }
                      if (iy < 0 || iy >= height) { continue; }
                      const float dSqr = u*u + v*v;
                      if (dSqr > radius*radius) { continue; }
                      if (bmap(ix,iy) != true) { continue; }
                      points(numPoints,0) = ix;
                      points(numPoints,1) = iy;
                      numPoints++;
                  }
              }
              const float theta = primaryOrient(points,numPoints);
              orientationMap(x,y) = theta;
          }
      }
  }

  /////////////////////////////////////////////////////////////////////////

  int 
  Segmentation::getWidth () const
  {
      return _map.size(0);
  }

  int 
  Segmentation::getHeight () const
  {
      return _map.size(1);
  }

  int 
  Segmentation::getNumSegments() const
  {
      return _numSegments;
  }

  //
  // get a list of the segment sizes
  //
  void
  Segmentation::segSizes(Array1D<int>& sizes)
  {
    _renumber();
    sizes.resize(_numSegments);
    sizes.init(0);
    for (int x = 0; x < _map.size(0); x++) 
    {
      for (int y = 0; y < _map.size(1); y++) 
      {
        sizes(_map(x,y))++;
      }
    }
  }

  /////////////////////////////////////////////////////////////////////////

  const char*
  Segmentation::getDate () const
  {
      return _date;
  }

  const char* 
  Segmentation::getUser () const
  {
      return _user;
  }

  const char* 
  Segmentation::getImage () const
  {
      return _image;
  }

  int 
  Segmentation::getImageID () const
  {
      return _imageID;
  }

  int 
  Segmentation::getUserID () const
  {
      return _userID;
  }

  int 
  Segmentation::getPresentation () const
  {
      return _presentation;
  }
  bool 
  Segmentation::isColor () const
  {
      return ! isGray ();
  }

  bool 
  Segmentation::isGray () const
  {
      return _presentation & presGray;
  }

  bool 
  Segmentation::isInverted () const
  {
      return _presentation & presInvert;
  }

  bool 
  Segmentation::isFlipFloped () const
  {
      return _presentation & presFlipFlop;
  }

}