segment.cc

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

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

static int
_main ()
{
    Util::Message::startBlock("segment_main");
    Util::Message::debug(Util::String("reading image %s",imageFileName));
    Util::ImageStack image;
    if (!Util::readJpegFile(imageFileName,image))
    {
      throw Util::Exception (Util::String ("Could not read JPEG file '%s'.", imageFileName));
    }

    const int width = image.size(1);
    const int height = image.size(2);
    const int idiag = (int)sqrt((double)(width*width+height*height));
    const int numPixels = width*height;
    Util::Message::debug(Util::String("image size [%d x %d]",width,height),1);

    Group::DualLattice boundaries; 
    Util::Message::startBlock("acquiring contour features");
    if (readLattice)
    {
      if (latticeFileName != NULL)
      {
        Util::Message::debug(Util::String("reading latticefile %s",latticeFileName));
        std::ifstream in(latticeFileName);
        if (!in.good()) {
          throw Util::Exception (Util::String ("Error opening '%s' for reading.", latticeFileName));
        }
        in >> boundaries.H;
        in >> boundaries.V;
        boundaries.width = boundaries.H.size(0);
        boundaries.height = boundaries.V.size(1);
        in.close();
      }      
    }
    else
    {
      Group::Pb pbFactory;
      pbFactory.initialize(image,useColor);
      pbFactory.computePb(8,boundaries);
      if (latticeFileName != NULL)
      {
        Util::Message::debug("writing lattice file");
        std::ofstream out(latticeFileName);
        if (!out.good()) {
          throw Util::Exception(Util::String ("Error opening '%s' for writing.", latticeFileName));
        }
        out << boundaries.H;
        out << boundaries.V;
        out.close();
      }
    }

    Group::Region patchFactory;
    int smapRadius = 0;
    if (cueMode == Group::ic)
    {
      smapRadius = dthresh;    
    }
    else
    {
      smapRadius = Util::max(patchFactory.getSupportMapRadius(idiag),dthresh);
    }
    Util::Message::debug(Util::String("running intervening contour with map radius = %d", smapRadius));
    Group::SupportMap ic;
    Group::computeSupport(boundaries,smapRadius,1.0f,ic);
    Util::Message::endBlock();

    Util::Message::startBlock("acquiring region features");
    if ((cueMode == Group::patch) || (cueMode == Group::both))
    {
      Util::Message::debug("region processing");
      patchFactory.initialize(image,ic);
    }
    Util::Message::endBlock();

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

    Util::Message::debug("building affinity matrix");
    SMatrix* affinities = NULL;
    Group::computeAffinities(ic,patchFactory,cueMode,sigma,dthresh,useColor,&affinities);
    assert(affinities != NULL);

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

    if (smatrixFileName != NULL)
    {
      Util::Message::debug("writing affinity matrix");
      FILE* fp = fopen(smatrixFileName,"w");
      affinities->dump(fp);
      fclose(fp);
    }

    if ((evecFileName != NULL) || (egradFileName != NULL) || (graphFileName != NULL) || (segFileName != NULL))
    {
      Util::Message::debug("computing normalized laplacian");
      double* rowSum = affinities->getRowSum();
      affinities->normalizedLaplacian(rowSum);
      
      Util::Message::debug("running eignesolver");
      double* eval;
      double* evec;
      affinities->eigs(numVectors+1,&eval,&evec);

      //output raw eigenvectors if desired
      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 = y*width+x;
              evecImage(vnum,x,y) = evec[vnum*numPixels + pnum];
            }
          }
        }

        Util::Message::debug("writing eigenvector file");
        std::ofstream out(evecFileName);
        if (!out.good()) {
          throw Util::Exception (Util::String ("Error opening '%s' for writing.", evecFileName));
        }
        out << evecImage;
        out.close();
      }

      //divide thru by leading eigenvector, then median filter as desired
      Util::ImageStack evecImage(numVectors,width,height);
      Util::Message::debug("scaling and filtering eigenvectors");
      for (int vnum = 1; vnum < (numVectors+1); vnum++)
      {
        Util::Image im(width,height);
        for (int x = 0; x < width; x++)
        {
          for (int y = 0; y < height; y++)
          {
            int pnum = y*width+x;
            im(x,y) = evec[vnum*numPixels + pnum] / evec[0*numPixels + pnum];
          }
        }
        float norm = 0.0f;
        if (filterEvec)
        { 
          Util::Image mim(width,height);
          getPctFiltered(im,3,0.5,mim);
          for (int x = 0;  x < width; x++)
          {
            for (int y = 0; y < height; y++)
            {
              evecImage(vnum-1,x,y) = mim(x,y);    
              norm = norm+im(x,y)*im(x,y);
            }
          }
        }
        else
        {
          for (int x = 0;  x < width; x++)
          {
            for (int y = 0; y < height; y++)
            {
              evecImage(vnum-1,x,y) = im(x,y);    
              norm = norm+im(x,y)*im(x,y);
            }
          }
        }
        norm = sqrt(norm);
        for (int x = 0;  x < width; x++)
        {
          for (int y = 0; y < height; y++)
          {
            evecImage(vnum-1,x,y) = evecImage(vnum-1,x,y)/norm;
          }
        }
      }

      if (egradFileName != NULL)
      {
        Util::ImageStack pbs;
        computeEigengradients(evecImage,pbs);
        Util::Message::debug("writing eigenvector-gradient file");
        std::ofstream out(egradFileName);
        if (!out.good()) {
          throw Util::Exception (Util::String ("Error opening '%s' for writing.", egradFileName));
        }
        out << pbs;
        out.close();
      }

      if ((graphFileName != NULL) || (segFileName != NULL))
      {
        Util::Message::debug("transposing eigenvectors");
        Util::Array2D<float> embed(width*height,numVectors);
        for (int x = 0; x < width; x++)
        {
          for (int y = 0; y < height; y++)
          {
            for (int vnum = 0; vnum < numVectors; vnum++)
            {
              embed(x*height+y,vnum) = evecImage(vnum,x,y);
            }
          }
        }

        affinities->undoNormalizedLaplacian(rowSum);

        Util::Message::startBlock(5,"finding superpixels");
        Util::Image means(numSuperpixels, numVectors);
        Util::Array1D<int> bestMembership(width*height);
        Util::kmeans(embed, numSuperpixels, 1e-4, 30, true, means, bestMembership);
        float bestNcut = affinities->computeNCut(rowSum,bestMembership,numSuperpixels);
        for (int iter = 0; iter < 10; iter++)
        {
          Util::Message::stepBlock();
          Util::Array1D<int> membership(width*height);
          Util::kmeans(embed, numSuperpixels, 1e-4, 30, true, means, membership);
          float ncut = affinities->computeNCut(rowSum,membership,numSuperpixels);
          if (ncut < bestNcut)
          {
            ncut = bestNcut;
            bestMembership = membership;
          }
        }
        Util::Message::endBlock();


        Util::Segmentation seg(width,height); 
        for (int x = 0; x < width; x++)
        {
          for (int y = 0; y < height; y++)
          {
            seg(x,y) = bestMembership(x*height+y);
          }
        }

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

        if (segFileName != NULL)
        {
          Util::Message::debug(Util::String("writing segmentation file %s",segFileName));
          std::ofstream out(segFileName);
          out << seg;
          out.close();
        }

        Util::Message::debug("coarsening graph G_1 -> G_2");
        Util::Array2D<float> Wsum(numSP,numSP);
        Wsum.init(0);
        Util::Array2D<float> D(numSP,numSP);
        D.init(0);
        for (int r = 0; r < affinities->n; r++)
        {
          int rx = r % width;
          int ry = (r - rx) / width;
          int rseg = seg(rx,ry);
          for (int i = 0; i < affinities->nz[r]; i++)
          {
            int c = affinities->col[r][i];       
            int cx = c % width;
            int cy = (c - cx) / width;
            int cseg = seg(cx,cy);
            Wsum(rseg,cseg) += affinities->values[r][i];
            D(rseg,cseg)++;
          }
        }

        Util::Message::debug("densifying G_2");
        Util::Array2D<float> W2(numSP,numSP);
        for (int i = 0; i < numSP; i++)
        {
          for (int j = 0; j < numSP; j++)
          {
            float dense = 0;
            if (D(i,j) != 0)
            {
              dense = segmentSizes(i)*(Wsum(i,j) / D(i,j))*segmentSizes(j);
            }
            W2(i,j) = dense;
          } 
        }

        if (graphFileName != NULL)
        {
          Util::Message::debug(Util::String("writing coarsened graph file %s",graphFileName));
          std::ofstream out(graphFileName);
          out << seg << Wsum << D << W2;
          out.close();
        }
     
     /*
        Util::Array2D<float> W3;
        Util::Segmentation seg3(width,height);
        densify(W2, seg, 80, seg3, W3);

        Util::Array2D<float> W4;
        Util::Segmentation seg4(width,height);
        densify(W3, seg3, 60, seg4, W4);

        Util::Array2D<float> W5;
        Util::Segmentation seg5(width,height);
        densify(W4, seg4, 30, seg5, W5);

        Util::Array2D<float> W6;
        Util::Segmentation seg6(width,height);
        densify(W5, seg5, 20, seg6, W6);

        Util::Array2D<float> W7;
        Util::Segmentation seg7(width,height);
        densify(W6, seg6, 10, seg7, W7);
      */  

        std::cerr << "done." << std::endl;
      }
    }

    return 0;
}

//////////////////////////////////////////////////////////////////////
//////////// Generic program wrapper below this point ////////////////
//////////////////////////////////////////////////////////////////////


// Print some useful info to stderr.
static void
infoBegin (const int argc, const char** argv)
{
    static char hostname[256];
    if (gethostname (hostname, 256) == 0) {
        Util::Message::debug(Util::String("host: %s",hostname),1);
    }

    static char cwd[1000];
    if (getcwd (cwd, 1000) != NULL) {
        Util::Message::debug(Util::String("pwd: %s",cwd),1);
    }

    Util::String args ("command: %s", argv[0]);
    for (int i = 1; i < argc; i++) {
        args.append (" %s", argv[i]);
    }
    Util::Message::debug(args,1);

    time_t t = time (NULL);
    Util::Message::debug(Util::String("start time: %s",ctime(&t)),1);
}

int
main (const int argc, const char** argv)
{
    int status = 0;
    try
    {
        // Register configuration options.
        registerConfig();

        // Initialize all modules.
        try
        {
            int index = Configure::init (argc, (const char**) argv);
            init (index, argc, argv);
            infoBegin (argc, argv);
            Configure::show();
        }
        catch (Util::Exception& e)
        {
            Util::Message::debug(Util::String("usage: %s  -image <image.jpg> [options]",argv[0]));
            Configure::usage();
            Util::Message::error(e.msg());
            exit (1);
        }
        // Run the main program.
        status = _main();
    }
    catch (Util::Exception& e)
    {
        // Some sort of uncaught runtime error.
        Util::Message::error(e.msg());
        Util::Message::debug("EXITING WITH ERROR!!!!");
        status = 0; //exit with status 0 even though there was a problem so that milrun still goes
    }
    return status;
}