segment.cc

来自「用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵」· CC 代码 · 共 797 行 · 第 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 <assert.h>
#include <stdio.h>
#include <string.h>
#include <fstream>
#include <iostream>

#include "types.hh"
#include "array.hh"
#include "util.hh"
#include "exception.hh"
#include "string.hh"
#include "configure.hh"
#include "kmeans.hh"
#include "segmentation.hh"

#include "smatrix.hh"
#include "dmatrix.hh"

#include "affinity.hh"
#include "pb.hh"
#include "ic.hh"
#include "region.hh"


//read/write cached Pb edge map?
static bool readLattice;
static const char* latticeFileName;

//write out intermediate results
static const char* evecFileName;
static const char* egradFileName;

//cues
static Group::cueModeType cueMode;
static const char* const cueModeStrings[] = {"ic","patch","both",NULL};
static bool useColor;
static int dthresh;
static float sigma;

//computing discrete ncut
static const char** segFiles;
static int numSegFiles;
static int dustThresh;

//affinity matrix
static const char* smatrixFileName;

//eigenvectors
static bool filterEvec;
static int numVectors;

//clustering
static int numSuperpixels;
static const char* segFileName;
static const char* graphFileName;

//static int numGroups;
//static const char* segFileName;

//the image
static const char* imageFileName;

static void
registerConfig ()
{
  Util::registerConfig();
  Group::Pb::registerConfig();
  Group::Region::registerConfig();

  Configure::registerString(
      "image", NULL,
      "Image file to proccess");

  Configure::registerBool(
      "readlattice", false,
      "Read in pb data?");

  Configure::registerString(
      "latticefile", NULL,
      "Pb lattice data file name");

  Configure::registerInt(
      "dthresh",5,"Connectivity distance threshold in pixels");

  Configure::registerFloat(
      "sigma",0.1,"Sigma non-linearity");

  Configure::registerBool(
      "usecolor", true,
      "Use color info.");

  Configure::registerEnum(
     "cues",cueModeStrings,Group::ic,
     "What cues to use in similarity computation");

  Configure::registerString(
      "smatrixfile", NULL,
      "Write out affinity matrix");

  Configure::registerInt(
      "numvectors", 16,
      "Number of eigenvectors.");

  Configure::registerBool(
      "filterevec", true,
      "Median filter eigenvectors.");

  Configure::registerString(
      "evecfile", NULL,
      "Write out eigenvector data");

  Configure::registerString(
      "egradfile", NULL,
      "Write out eigenvector gradient data");

  Configure::registerInt(
      "numsuperpixels", 50,
      "Number of superpixels.");

  Configure::registerString(
      "segfile", NULL,
      "Write out superpixel segmentation");

  Configure::registerString(
      "graphfile", NULL,
      "Write out collapsed graph data");

  Configure::registerInt(
      "dustthresh", 10,
      "How small before we merge a disconnected component.");
}

static void
init (const int index, const int argc, const char** argv)
{
  Util::init();

  segFiles = &argv[index];
  numSegFiles = argc - index;

  imageFileName = Configure::getString("image");

  latticeFileName = Configure::getString("latticefile"); 
  readLattice = Configure::getBool("readlattice");

  dthresh = Configure::getInt("dthresh");
  sigma = Configure::getFloat("sigma");
  useColor = Configure::getBool("usecolor");
  cueMode = (Group::cueModeType)Configure::getEnum("cues");
  smatrixFileName = Configure::getString("smatrixfile");

  numVectors = Configure::getInt("numvectors");
  filterEvec = Configure::getBool("filterevec");
  egradFileName = Configure::getString("egradfile");
  evecFileName = Configure::getString("evecfile");

  graphFileName = Configure::getString("graphfile");
  numSuperpixels = Configure::getInt("numsuperpixels");
  dustThresh = Configure::getInt("dustthresh");
  segFileName = Configure::getString("segfile");

  if (imageFileName == NULL)
  {
    throw Util::Exception (Util::String ("image filename is required"));
  }
}

//////////////////////////////////////////////////////////////////////////////
//
// compute gradients on telescoping sets of eigenvectors and run
// non-maxima supression on the results.
//
void computeEigengradients(const Util::ImageStack& evec, Util::ImageStack& pbs)
{
  const int numVectors = evec.size(0);
  const int width = evec.size(1);
  const int height = evec.size(2);

  Util::Image g1(width,height); g1.init(0);
  Util::Image g2(width,height); g2.init(0);
  Util::Image g12(width,height); g12.init(0);
  Util::Image mag(width,height); 
  Util::Image ori(width,height); 
  Util::Image dx(width,height); 
  Util::Image dy(width,height);
  Util::Image gradient(width,height);

  pbs.resize(numVectors,width,height);

  //create filters
  double smooth = 0.5;
  Util::Image hFilt;
  Util::Image vFilt;
  Group::FilterBank::createGaussianKernel(smooth, smooth, 10.0, 1.5708, 1, false,hFilt);
  Group::FilterBank::createGaussianKernel(smooth, smooth, 10.0, 0, 1, false,vFilt);

  //filter each set of vectors up to i
  for (int i = 0; i < numVectors; i++)
  {
    getFiltered(*evec.slice(i),hFilt,dx);
    getFiltered(*evec.slice(i),vFilt,dy);
    for (int x = 0; x < width; x++)
    {
      for (int y = 0; y < height; y++)
      {
        double f1 = g1(x,y) + dx(x,y)*dx(x,y);
        double f2 = g2(x,y) + dy(x,y)*dy(x,y);
        double f12 = g12(x,y) + dx(x,y)*dy(x,y);
        double orival = 0.5*atan(2*f12 / (f1 - f2));
        if ((f1-f2) < 0)
        {
          if (2*f12>0)
          {
            orival+=M_PI/2;
          }
          else
          {
            orival-=M_PI/2;
          }
        }
        mag(x,y) =  0.5*(f1 + f2 + sqrt((f1-f2)*(f1-f2) + 4*f12*f12));
        ori(x,y) =  fmod(-orival+(M_PI/2),M_PI);
        g1(x,y) = f1;
        g2(x,y) = f2;
        g12(x,y) = f12;
      }
    }
    Group::nonmaxSuppress(mag,ori,gradient);
    for (int x = 0; x < width; x++)
    {
      for (int y = 0; y < height; y++)
      {
        pbs(i,x,y) = gradient(x,y);
      }
    }
  }
}


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

void 
densify(const Util::Array2D<float>& W, const Util::Segmentation& seg, const int desiredNumSP, 
        Util::Segmentation& coarseSeg, Util::Array2D<float>& Wcoarse)
{
    static int level = 0;
    const int width = seg.getWidth();
    const int height = seg.getHeight();
    const int numP = W.size(0);
    level++;

    Util::Message::debug("computing eigenvectors of G");
    DMatrix* DW = new DMatrix(W);
    double* rowSum = DW->getRowSum();
    DW->normalizedLaplacian(rowSum);
    double* devec;
    double* deval;
    DW->eig(numVectors+1,&devec,&deval);
    DW->undoNormalizedLaplacian(rowSum);

    Util::Message::startBlock("finding superpixels");
    Util::Message::debug("transposing eigenvectors");
    Util::Array2D<float> embed(numP,numVectors);
    for (int vnum = 0; vnum < numVectors; vnum++)
    {
      for (int pnum = 0; pnum < numP; pnum++)
      {
        embed(pnum,vnum) = devec[(vnum+1)*numP + pnum] / devec[0*numP + pnum];
      }
    }

    Util::Message::startBlock(50,Util::String("running kmeans k=%d",desiredNumSP));
    Util::Image means(desiredNumSP, numVectors);
    Util::Array1D<int> bestMembership(numP);
    Util::kmeans(embed, desiredNumSP, 1e-4, 30, true, means, bestMembership);
    float bestNcut = DW->computeNCut(rowSum,bestMembership);
    for (int iter = 0; iter < 50; iter++)
    {
      Util::Message::stepBlock();
      Util::Array1D<int> membership(width*height);
      Util::kmeans(embed, desiredNumSP, 1e-4, 30, true, means, membership);
      float ncut = DW->computeNCut(rowSum,membership);
      if (ncut < bestNcut)
      {
        ncut = bestNcut;
        bestMembership = membership;
      }
    }
    Util::Message::endBlock();

    for (int x = 0; x < width; x++)
    {
      for (int y = 0; y < height; y++)
      {
        coarseSeg(x,y) = bestMembership(seg(x,y));
      }
    }

    Util::Message::debug("post-processing segmentation");
    coarseSeg.fragment();
    coarseSeg.mergeDust(dustThresh);
    Util::Array1D<int> segmentSizes;
    coarseSeg.segSizes(segmentSizes);
    int numSP = segmentSizes.size(0);
    Util::Message::debug(Util::String("found %d superpixels",numSP));
    Util::Message::endBlock();

    Util::Message::debug("coarsening graph");
    Util::Array1D<int> p2sp(numP);
    for (int x = 0; x < width; x++)
    {
      for (int y = 0; y < height; y++)
      {
        p2sp(seg(x,y)) = coarseSeg(x,y);
      }
    }
    Util::Array2D<float> Wsum(numSP,numSP); 
    Wsum.init(0);
    Util::Array2D<float> D(numSP,numSP); 
    D.init(0);
    for (int i = 0; i < numP; i++)
    {
      for (int j = 0; j < numP; j++)
      {
        //collapse the edge from i to j into the superpixel edge
        Wsum(p2sp(i),p2sp(j)) += W(i,j);
        D(p2sp(i),p2sp(j))++; //TODO???? should we be adding a different degree here?
      } 
    }

    Util::Message::debug("densifying graph");
    Wcoarse.resize(numSP,numSP);
    for (int i = 0; i < numSP; i++)
    {
      for (int j = 0; j < numSP; j++)
      {
        double dense = 0;
        if (D(i,j) != 0)
        {
          dense = segmentSizes(i)*(Wsum(i,j) / D(i,j))*segmentSizes(j);
        }
        Wcoarse(i,j) = dense;
      } 
    }

    if (evecFileName != NULL)
    {
      Util::ImageStack evecImage(numVectors,width,height);
      for (int vnum = 0; vnum < numVectors; vnum++)
      {
        for (int x = 0; x < width; x++)
        {
          for (int y = 0; y < height; y++)
          {
            int pnum = seg(x,y);
            evecImage(vnum,x,y) = devec[vnum*numP + pnum];
          }
        }
      }
      char filename[1000];
      sprintf(filename,"%s.%d",evecFileName,level);
      Util::Message::debug("writing eigenvector file");
      std::ofstream out(filename);
      if (!out.good()) {
        throw Util::Exception (Util::String ("Error opening '%s' for writing.", evecFileName));
      }
      out << evecImage;
      out.close();
    }

    if (graphFileName != NULL)
    {
      char filename[1000];
      sprintf(filename,"%s.%d",graphFileName,level);
      Util::Message::debug("writing coarsened graph file");
      std::ofstream out(filename);
      out << seg << W << Wsum << D << Wcoarse << coarseSeg;
      out.close();
    }
}