pb.cc

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

// Copyright (C) 2002 Charless C. Fowlkes <[e-mail]>
// Copyright (C) 2002 David R. Martin <[e-mail]>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
// 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html.


#include <float.h>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <math.h>
#include <string.h>

#include "configure.hh"
#include "exception.hh"
#include "message.hh"
#include "string.hh"
#include "array.hh"
#include "util.hh"

#include "image.hh"
#include "filterbank.hh"
#include "nonmax.hh"
#include "texture.hh"
#include "color.hh"
#include "pb.hh"

namespace Group
{

  //////////////////////////////////////////////////////////////////////////////
  // 
  // configuration structure and interface
  //

  bool Pb::registerConfigCalled = false;

  void Pb::registerConfig()
  {
    if(registerConfigCalled)
    {
      return;
    }

    
    // filterbank for tg
    Configure::registerInt("Pb::tg::filt::norient",
                           6, "Number of orientations for filterbank.");
    Configure::registerInt("Pb::tg::filt::nscales", 1,
                           "Number of scales for filterbank.");
    Configure::registerFloat("Pb::tg::filt::scalefactor", sqrt(2),
                             "Scale factor for filterbank.");
    Configure::registerFloat("Pb::tg::filt::startscale", 0.007,
                             "Beginning scale for filterbank, as fraction of image diagonal.");

    // textons
    Configure::registerString("Pb::textons::file", NULL,
                              "If provided, textons are read from the file rather than computed.");
    Configure::registerInt("Pb::textons::startchannels", 12,
                           "Beginning number of texton channels.");
    Configure::registerFloat("Pb::textons::channelfactor", 2,
                             "Increase number of textons by this factor for each scale.");

    // texture gradient
    Configure::registerBool("Pb::tg::scalelinked", false,
                            "Use textons computed from different scale filters for different "
                            "scale texture gradients?");
    Configure::registerInt("Pb::tg::norient", 8,
                           "Number of orientations for texture gradient.");
    Configure::registerInt("Pb::tg::nscales", 1,
                           "Number of scales for texture gradient.");
    Configure::registerFloat("Pb::tg::startscale", 0.0198,
                             "Starting scale for texture gradient.");
    Configure::registerFloat("Pb::tg::scalefactor", sqrt(2),
                             "Ratio between successive texture scales.");
    Configure::registerBool("Pb::tg::gaussDeriv", false,
                            "Use gaussian derivative estimates (instead of Savitsky-Golay filtering)?");
    Configure::registerFloat("Pb::tg::epsilon", 0.01,
                             "Epsilon for needleness function.");

    // color gradient
    Configure::registerInt("Pb::cg::norient", 8,
                           "Number of orientations for CG.");
    Configure::registerInt("Pb::cg::nscales", 1,
                           "Number of scales for CG.");
    Configure::registerFloat("Pb::cg::startscale", 0.0198,
                             "Starting scale for CG.");
    Configure::registerFloat("Pb::cg::scalefactor", sqrt(2),
                             "Ratio between successive CG scales.");
    Configure::registerInt("Pb::cg::bins::a", 25,
                           "Number of bins for a channel.");
    Configure::registerInt("Pb::cg::bins::b", 25,
                           "Number of bins for b channel.");
    Configure::registerFloat("Pb::cg::sigma", 0.10,
                             "Sigma for CG histogram kernel.");
    Configure::registerFloat("Pb::cg::support", 2,
                             "Support of CG histogram kernel, in units of sigma.");
    Configure::registerFloat("Pb::cg::zoom", 5,
                             "Resolution of CG histogram kernel, in units of points/sigma.");

    // brightness gradient
    Configure::registerInt("Pb::bg::norient", 8,
                           "Number of orientations for BG.");
    Configure::registerInt("Pb::bg::nscales", 1,
                           "Number of scales for BG.");
    Configure::registerFloat("Pb::bg::startscale", 0.0106,
                             "Starting scale for BG.");
    Configure::registerFloat("Pb::bg::scalefactor", sqrt(2),
                             "Ratio between successive BG scales.");
    Configure::registerInt("Pb::bg::bins::L", 12,
                           "Number of bins for L channel.");
    Configure::registerFloat("Pb::bg::sigma", 0.20,
                             "Sigma for BG histogram kernel.");
    Configure::registerFloat("Pb::bg::support", 2,
                             "Support of BG histogram kernel, in units of sigma.");
    Configure::registerFloat("Pb::bg::zoom", 5,
                             "Resolution of BG histogram kernel, in units of points/sigma.");


    registerConfigCalled = true;
  }


  void Pb::getDefaultConfig(Config & config)
  {

    // filterbank for tg
    config.tgFiltNumOrient = Configure::getInt("Pb::tg::filt::norient");
    config.tgFiltNumScales = Configure::getInt("Pb::tg::filt::nscales");
    config.tgFiltScaleFactor =
      Configure::getFloat("Pb::tg::filt::scalefactor");
    config.tgFiltStartScale =
      Configure::getFloat("Pb::tg::filt::startscale");
    config.gaussDeriv = Configure::getBool("Pb::tg::gaussDeriv");
    config.tgEpsilon = Configure::getFloat("Pb::tg::epsilon");

    // textons
    config.textonFile = Configure::getString("Pb::textons::file");
    config.textonStartChannels =
      Configure::getInt("Pb::textons::startchannels");
    config.textonChannelFactor =
      Configure::getFloat("Pb::textons::channelfactor");

    // texture gradient
    config.tgScaleLinked = Configure::getBool("Pb::tg::scalelinked");
    config.tgNumOrient = Configure::getInt("Pb::tg::norient");
    config.tgNumScales = Configure::getInt("Pb::tg::nscales");
    config.tgStartScale = Configure::getFloat("Pb::tg::startscale");
    config.tgScaleFactor = Configure::getFloat("Pb::tg::scalefactor");

    // color gradient
    config.cgNumOrient = Configure::getInt("Pb::cg::norient");
    config.cgNumScales = Configure::getInt("Pb::cg::nscales");
    config.cgStartScale = Configure::getFloat("Pb::cg::startscale");
    config.cgScaleFactor = Configure::getFloat("Pb::cg::scalefactor");
    config.cgBinsA = Configure::getInt("Pb::cg::bins::a");
    config.cgBinsB = Configure::getInt("Pb::cg::bins::b");
    config.cgSigma = Configure::getFloat("Pb::cg::sigma");
    config.cgSupport = Configure::getFloat("Pb::cg::support");
    config.cgZoom = Configure::getFloat("Pb::cg::zoom");

    // brightness gradient
    config.bgNumOrient = Configure::getInt("Pb::bg::norient");
    config.bgNumScales = Configure::getInt("Pb::bg::nscales");
    config.bgStartScale = Configure::getFloat("Pb::bg::startscale");
    config.bgScaleFactor = Configure::getFloat("Pb::bg::scalefactor");
    config.bgBinsL = Configure::getInt("Pb::bg::bins::L");
    config.bgSigma = Configure::getFloat("Pb::bg::sigma");
    config.bgSupport = Configure::getFloat("Pb::bg::support");
    config.bgZoom = Configure::getFloat("Pb::bg::zoom");

  }

  void Pb::checkConfig()
  {
    // TODO: if any of the configuration parameters are invalid, then 
    // throw an exception describing the problem.
  }

  const Pb::Config& Pb::getConfig () const
  {
    return config;
  }

  Pb::Pb()
  {
    getDefaultConfig(config);
  }

  Pb::Pb(const Config & config)
  {
    this->config = config;
  }

  Pb::~Pb()
  {
  }

  //////////////////////////////////////////////////////////////////////////////
  //
  // initialize does all the computationally expensive
  // image processing in order to compute the various
  // gradient features.  it can be called multiple times
  // with different images if desired....
  //
  
  void Pb::initialize(const Util::ImageStack& im, const bool useColor)
  {
    checkConfig();

    m_useColor = useColor;
    m_width = im.size(1);
    m_height = im.size(2);
    m_idiag = sqrt(m_width * m_width + m_height * m_height);
    
    if (im.size(0) == 3)
    {
      rgb2lab(im,m_lab);
      labNormalize(m_lab);
    } 
    else
    {
      if (m_useColor == true)
      {
        Util::Message::debug("WARNING: cowardly refusing to compute color pb for grayscale image");  
      }
      m_lab = im;
      m_useColor = false;
    }

    Util::Message::startBlock("Pb initialization");

    Util::Message::startBlock("brightness gradient computation");
    computeBG();
    Util::Message::endBlock();
    if (m_useColor)
    {
      Util::Message::startBlock("color gradient computation");
      computeCG();
      Util::Message::endBlock();
    }

    Util::Message::startBlock("texture gradient computation");
    computeTG();
    Util::Message::endBlock();

    Util::Message::endBlock();

    m_initialized = true;
  }

  //////////////////////////////////////////////////////////////////////////////
  //
  // compute the pb image from the cached gradients
  //

  //
  // logistic fits.  this is of course only valid for the default 
  // parameter settings!!  changing relative scales etc. can make
  // these fits "non-optimal"
  //
  static const float bg_stdev = 0.13;
  static const float cg_stdev = 0.077;
  static const float tg_stdev = 0.063;

  static float cpb_model(const Util::Array1D<float> features)
  {
    float val = -3.08;
    val += features(0) / bg_stdev * 0.31;  //bg
    val += features(1) / tg_stdev * 0.53;  //cg
    val += features(2) / cg_stdev * 0.44;  //tg
    float post = 1.0 / (1.0 + exp(-val));
    if(!finite(post)) { return 0; }
    post = Util::minmax(0.0f,post,1.0f);
    return post;
  }

  static float gpb_model(const Util::Array1D<float> features)
  {
    float val = -2.81;
    val += features(0) / bg_stdev * 0.50;  // bg 
    val += features(1) / tg_stdev * 0.52;  // tg 
    float post = 1.0 / (1.0 + exp(-val));
    if(!finite(post)) { return 0; }
    post = Util::minmax(0.0f,post,1.0f);
    return post;
  }
  
  static int closestOrient(const float theta, const int norient)
  {
    assert(theta >= 0);
    int orient = (int)rint((theta / M_PI) * norient);
    assert(orient >= 0);
    return (orient % norient);
  }

  //
  // compute probability of a boundary at a given location and
  // orientation.  this assumes the default configuration parameters.  
  // in particular, it only uses the zeroth scale.
  //
  void Pb::computeOrientedPb(const int numOrient, Util::ImageStack& pbs)
  {
    assert(m_initialized);
    pbs.resize(numOrient,m_width,m_height);
    Util::Array1D<float> features;
    if(m_useColor)
    {
      features.resize(3);
    }
    else
    {
      features.resize(2);
    }

    features.init(0);
    for(int orient = 0; orient < numOrient; orient++)
    {
      const float theta = (M_PI*orient / numOrient);
      const int bgorient = closestOrient(theta, config.bgNumOrient);
      const int cgorient = closestOrient(theta, config.cgNumOrient);
      const int tgorient = closestOrient(theta, config.tgNumOrient);
      for(int x = 0; x < m_width; x++)
      {
        for(int y = 0; y < m_height; y++)
        {
          float p = 0;
          if (m_useColor)
          {
            features(0) = m_bg(bgorient,x,y);
            features(1) = m_tg(tgorient,x,y);
            features(2) = m_cg(cgorient,x,y);
            p = cpb_model(features);
          }
          else
          {
            features(0) = m_bg(bgorient,x,y);
            features(1) = m_tg(tgorient,x,y);
            p = gpb_model(features);
          }
          assert(p >= 0 && p <= 1);
          pbs(orient,x,y) = p;
        }
      }
    }
  }

  //
  // compute probability of a boundary at a given location and