pb.cc

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

  // orientation.  this assumes the default configuration parameters.  
  // in particular, it only uses the zeroth scale.
  //
  void Pb::computeOrientedNMSPb(const int numOrient, Util::ImageStack& pbs)
  {
    assert(m_initialized);
    computeOrientedPb(numOrient,pbs);

    const float oStep = M_PI / numOrient;
    const float r = 2.5;
    for (int orient = 0; orient < numOrient; orient++)
    {
      float theta = orient*oStep;
      assert(theta >= 0 && theta < M_PI);
      ParabolaSet parabs;
      fitSubpixelParabolas(*pbs.slice(orient), theta, r, parabs);
      Util::Image mag = parabs.smooth;

      Util::Image nms(pbs.size(1),pbs.size(2));
      nonmaxSuppress(mag,theta,nms);
      for (int x = 0; x < pbs.size(1); x++)
      {
        for (int y = 0; y < pbs.size(2); y++)
        {
          pbs(orient,x,y) = nms(x,y); 
        }
      }
    }                 
  }

  //
  // compute pb maxed over orientations
  //
  void Pb::computePb(const int numOrient, DualLattice& boundaries)
  {
    assert(m_initialized);
    Util::ImageStack pbs(numOrient,m_width,m_height);
    computeOrientedPb(numOrient,pbs);

    // non-maximal suppression and sub-pixel localization
    Util::Message::debug("fitting sub-pixel parabolas");
    //const float fitRad = 0.0104*m_idiag;
    const float fitRad = 0.01*m_idiag;
    ParabolaSet parabs;
    fitSubpixelParabolas(pbs,fitRad,parabs);

    // intersect needles to get final pb
    Util::Message::debug("projecting parabolas down to pixel grid");
    intersectPixels(parabs, boundaries);
  }

  void Pb::computePb(const int numOrient, Util::Image& boundaries)
  {
    assert(m_initialized);
    Util::ImageStack pbs(numOrient,m_width,m_height);
    computeOrientedPb(numOrient,pbs);

    // non-maximal suppression and sub-pixel localization
    Util::Message::debug("fitting sub-pixel parabolas");
    //const float fitRad = 0.0104*m_idiag;
    const float fitRad = 0.01*m_idiag;
    ParabolaSet parabs;
    fitSubpixelParabolas(pbs,fitRad,parabs);

    // intersect needles to get final pb
    Util::Message::debug("projecting parabolas down to pixel grid");
    intersectPixels(parabs, boundaries);
  }

  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // access to raw gradient info and texton maps
  //

  const TextonMap& Pb::getTextons(const int scale) const
  {
    assert(m_initialized);
    assert(scale >=0 && scale < config.tgNumScales);
    return m_textons(scale);
  }
  
  const Util::Image& Pb::getTG(const int scale, const int orient) const
  {
    assert(m_initialized);
    assert(scale >= 0 && scale < config.tgNumScales);
    assert(orient >= 0 && orient < config.tgNumOrient);
    return *m_tg.slice(scale*config.tgNumOrient + orient);
  }

  const Util::Image& Pb::getCG(const int scale, const int orient) const
  {
    assert(m_initialized);
    assert(scale >= 0 && scale < config.cgNumScales);
    assert(orient >= 0 && orient < config.cgNumOrient);
    if (m_useColor)
    {
      return *m_cg.slice(scale*config.cgNumOrient + orient);
    }
    else
    {
      throw Util::Exception("no color gradient available");
    }
  }

  const Util::Image& Pb::getBG(const int scale, const int orient) const
  {
    assert(m_initialized);
    assert(scale >= 0 && scale < config.bgNumScales);
    assert(orient >= 0 && orient < config.bgNumOrient);
    return *m_bg.slice(scale*config.bgNumOrient + orient);
  }


  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // hints for accessing TG,CG,BG
  //
  
  inline int Pb::getTGNumScales() const
  {
    return config.tgNumScales;
  }
  inline int Pb::getTGNumOrient() const
  {
    return config.tgNumOrient;
  }
  inline int Pb::getCGNumScales() const
  {
    return config.cgNumScales;
  }
  inline int Pb::getCGNumOrient() const
  {
    return config.cgNumOrient;
  }
  inline int Pb::getBGNumScales() const
  {
    return config.bgNumScales;
  }
  inline int Pb::getBGNumOrient() const
  {
    return config.bgNumOrient;
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // helper functions for computing features and pb
  //

  //
  // compute "needleness" which is a function of the derivatives at a point
  //
  void Pb::computeNeedleness(const float epsilon, const Util::Image& a, const Util::Image& b, 
                                 const Util::Image& c, Util::Image& needleness)
  {
    const int width = a.size(0);
    const int height = a.size(1);
    needleness.resize(width,height);
    for(int x = 0; x < width; x++)
    {
      for(int y = 0; y < height; y++)
      {
        needleness(x,y) = Util::max(0.0f,a(x,y)) * Util::max(0.0f,-c(x,y)) / (fabs(b(x,y))+epsilon);
      }
    }
  }

  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // core computation
  //

  //
  // compute texture gradient
  //
  void Pb::computeTG()
  {
    m_tg.resize(config.tgNumScales*config.tgNumOrient,m_width,m_height);
    m_textons.resize(config.tgNumScales);
    Util::Image tgSmooth(m_width,m_height);
    Util::Image tgSlope(m_width,m_height);
    Util::Image tgCurv(m_width,m_height);

    // compute TG at each scale
    for(int tgScale = 0; tgScale < config.tgNumScales; tgScale++)
    {
      const float tgScaleVal = m_idiag * config.tgStartScale * pow(config.tgScaleFactor, tgScale);
      float filtStartScale = m_idiag * config.tgFiltStartScale;
      if(config.tgScaleLinked)
      {
        filtStartScale *= pow(config.tgScaleFactor, tgScale);
      }

      // create the filterbank
      FilterBank fbank(config.tgFiltNumOrient, config.tgFiltNumScales,
                       filtStartScale, config.tgFiltScaleFactor);

      Util::Message::debug(Util::String("Filterbank %d %d %f %f",
         config.tgFiltNumOrient,config.tgFiltNumScales, filtStartScale,config.tgFiltScaleFactor),2); 

      //extract textons
      m_textons(tgScale).map.resize(m_lab.size(1),m_lab.size(2));
      m_textons(tgScale).numChannels = 0;

      if(config.textonFile == "")
      {
        const int nchannels = (int)rint(config.textonStartChannels * pow(config.textonChannelFactor,tgScale));
        computeTextons(*m_lab.slice(Util::LAB_L),fbank,nchannels,m_textons(tgScale));
      }
      else
      {
        computeTextons(*m_lab.slice(Util::LAB_L),fbank,config.textonFile,m_textons(tgScale));
      }
      Util::Message::debug(Util::String("Used %d textons",m_textons(tgScale).numChannels),2);

      // compute the TG
      Util::ImageStack tgraw;
      Group::computeTG(m_textons(tgScale),tgScaleVal,config.tgNumOrient,tgraw);

      // compute localized TG
      Util::Message::startBlock(config.tgNumOrient,Util::String("TG %d smoothing",tgScale).text());
      for(int orient = 0; orient < config.tgNumOrient; orient++)
      {
        Util::Message::stepBlock();
        int index = tgScale * config.tgNumOrient + orient;
        const float r = tgFitRadius(tgScale);
        const float theta = M_PI * orient / config.tgNumOrient;
        if (config.gaussDeriv)
        {
          derivs(theta, r, *tgraw.slice(index), tgSmooth, tgSlope, tgCurv);
        }
        else
        {
          lowess(theta, r, *tgraw.slice(index), tgSmooth, tgSlope, tgCurv);
        }

        Util::Image needle;
        computeNeedleness(config.tgEpsilon, tgSmooth, tgSlope, tgCurv, needle);
        for (int x = 0; x < m_width; x++)
        {
          for (int y = 0; y < m_height; y++)
          {
            m_tg(index,x,y) = needle(x,y);
          }   
        }
      }
      Util::Message::endBlock();
    }
  }


  // 
  // compute color gradient
  //
  void Pb::computeCG()
  {
    assert(m_useColor);

    // allocate
    m_cg.resize(config.cgNumScales*config.cgNumOrient,m_width,m_height);
    Util::ImageStack cgA;
    Util::ImageStack cgB;

    // compute CG at each scale
    for(int cgScale = 0; cgScale < config.cgNumScales; cgScale++)
    {
      const float cgScaleVal = m_idiag * config.cgStartScale * pow(config.cgScaleFactor, cgScale);
      Group::computeCG(*m_lab.slice(Util::LAB_A), config.cgBinsA, cgScaleVal, config.cgNumOrient,
                config.cgSigma, config.cgSupport, config.cgZoom, cgA);
      Group::computeCG(*m_lab.slice(Util::LAB_B), config.cgBinsB, cgScaleVal, config.cgNumOrient,
                config.cgSigma, config.cgSupport, config.cgZoom, cgB);

      for(int orient = 0; orient < config.cgNumOrient; orient++)
      {
        const int index = cgScale * config.cgNumOrient + orient;
        for(int x = 0; x < m_width; x++)
        {
          for(int y = 0; y < m_height; y++)
          {
            m_cg(index,x,y) = 0.5 * (cgA(orient,x,y) + cgB(orient,x,y));
          }
        }
      }
    }

  }


  //
  // compute brightness gradient
  //
  void Pb::computeBG()
  {
    // allocate
    m_bg.resize(config.bgNumScales * config.bgNumOrient,m_width,m_height);
    Util::ImageStack bg;

    // compute BG at each scale
    for(int bgScale = 0; bgScale < config.bgNumScales; bgScale++)
    {
      const float bgScaleVal = m_idiag * config.bgStartScale * pow(config.bgScaleFactor, bgScale);
      Group::computeCG(*m_lab.slice(Util::LAB_L),config.bgBinsL, bgScaleVal, config.bgNumOrient,
                config.bgSigma, config.bgSupport, config.bgZoom, bg);

      for(int orient = 0; orient < config.bgNumOrient; orient++)
      {
        const int index = bgScale * config.bgNumOrient + orient;
        for(int x = 0; x < m_width; x++)
        {
          for(int y = 0; y < m_height; y++)
          {
            m_bg(index,x,y) = bg(orient,x,y);
          }
        }
      }
    }
  }


  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // helper functions for computing features and pb
  //

  float Pb::tgFitRadius(const int scale) const
  {
    const float alpha = 1.0;
    const float beta = 1.0;
    assert(scale >= 0 && scale <= config.tgNumScales);

    // compute disc radius used to compute TG at this scale 
    const float rDisc = config.tgStartScale
      * pow(config.tgScaleFactor, scale);

    // get the max sigma over filterbank used for TG
    const float sigmaFilt = config.tgFiltStartScale
      * pow(config.tgFiltScaleFactor, config.tgFiltNumScales);

    // return radius for fit, in pixels
    const float r = m_idiag * (alpha * sigmaFilt + beta * rDisc);
    return Util::max(2.1f, r);
  }

  float Pb::cgFitRadius(const int scale) const
  {
    const float r = m_idiag * config.cgStartScale
      * pow(config.cgScaleFactor, scale);
     return Util::max(2.1f, r);
  }

  float Pb::bgFitRadius(const int scale) const
  {
    const float r = m_idiag * config.bgStartScale
      * pow(config.bgScaleFactor, scale);
     return Util::max(2.1f, r);
  }

}  //namespace Group