construct_graph.cpp

A salient-boundary extraction software package based on the paper: S. Wang, T. Kubota, J. M. Siskind

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>

using namespace std;

const double MAX_DISTANCE = 50;
const double PI = 3.141592653589793;

int dround(double);
double distance(double x1, double y1, double x2, double y2);
int minimal(double, double, double, double);
double minimum(double, double, double, double);
int sign(double);
double area(double P1x, double P1y, double P2x, double P2y);
double triangle_area(double Ax, double Ay, double Bx, double By, double Cx, double Cy);
bool quad_inside(double p1x,double p1y,double p2x,double p2y,double p3x,double p3y,double p4x,double p4y,double px,double py);

double curvature(double P0x, double P0y, double P1x, double P1y, double P2x, double P2y, double P3x, double P3y);
double* triangle_incenter(double Ax,double Ay,double Bx,double By,double Cx,double Cy);
double* find_intersection(double P1x, double P1y, double P2x, double P2y, double P3x, double P3y, double P4x, double P4y);
double* find_corner(double P1x, double P1y, double P2x, double P2y, double Ix, double Iy, vector<double*> Ctable);
double angle(double p1x, double p1y, double p2x, double p2y, double p3x, double p3y);

double gradient(double P1x, double P1y, double P2x, double P2y, vector<double*> G);

//int main(int argc, char **argv) {
vector<double*>* construct_graph(vector<double*> Ltable, vector<double*> Ctable, vector<double*> grad_img, double WIDTH, double LAMBDA, bool homogeneity, bool use_corners, string datafile, double ALPHA) {

   double *dataread1, *dataread2;
   double *dataout;
   vector<double*> *Gtable = new vector<double*>;
   double MAX_W1 = 0, MAX_W2 = 0;
   ofstream output;
   int PSI = 0;
   if (homogeneity) PSI = 1;

   if (use_corners) output.open(datafile.c_str());
   

   // PART 1 - SOLID EDGE CONSTRUCTION //

   for (int i=0; i<(int)Ltable.size(); i++) {

      dataout = new double[5];
      dataout[0] = 4*i;
      dataout[1] = 4*i + 2;
      dataout[2] = 0; 
      dataout[3] = 0;
      dataout[4] = 1; // solid
      (*Gtable).push_back(dataout);

      dataout = new double[5];
      dataout[0] = 4*i + 1;
      dataout[1] = 4*i + 3;
      dataout[2] = 0;
      dataout[3] = 0;
      dataout[4] = 1; // solid
      (*Gtable).push_back(dataout);
   }

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

   // PART 2 - DASHED EDGE CONSTRUCTION //

   double P1[2], back1[2], P2[2], back2[2];
   double M1[2], M2[2], *C, *I, curv;
   double D, w1, area1, area2, carea, gaparea, alpha, beta;
   double G1 = 0, G2 = 0, gapG = 0;
   int k, n1, n2, count_dashed = 0;

   for (int i=0; i<(int)Ltable.size()-1; i++) {
      dataread1 = Ltable[i];

      for (int j=i+1; j<(int)Ltable.size(); j++) {
         dataread2 = Ltable[j];

         // find closest axis endpoints.
         D = minimum(distance(dataread1[0],dataread1[1],dataread2[0],dataread2[1]),distance(dataread1[0],dataread1[1],dataread2[2],dataread2[3]),distance(dataread1[2],dataread1[3],dataread2[0],dataread2[1]),distance(dataread1[2],dataread1[3],dataread2[2],dataread2[3]));
         k = minimal(distance(dataread1[0],dataread1[1],dataread2[0],dataread2[1]),distance(dataread1[0],dataread1[1],dataread2[2],dataread2[3]),distance(dataread1[2],dataread1[3],dataread2[0],dataread2[1]),distance(dataread1[2],dataread1[3],dataread2[2],dataread2[3]));
         // D is min distance, and k is index 1-4, e.g. 1 means O_2i and O_2j

         if (D > MAX_DISTANCE) continue;


         if (k < 3) {
            P1[0] = dataread1[0];
            P1[1] = dataread1[1];
            back1[0] = dataread1[2];
            back1[1] = dataread1[3];
         } else {
            P1[0] = dataread1[2];
            P1[1] = dataread1[3];
            back1[0] = dataread1[0];
            back1[1] = dataread1[1];
         }
         n1 = 2*i;
         area1 = area(P1[0],P1[1],back1[0],back1[1]);
         if (homogeneity)
            G1 = gradient(P1[0],P1[1],back1[0],back1[1], grad_img);
         if (P1[0] > back1[0])
            n1++;

         if (k % 2 == 1) {
            P2[0] = dataread2[0];
            P2[1] = dataread2[1];
            back2[0] = dataread2[2];
            back2[1] = dataread2[3];
         } else {
            P2[0] = dataread2[2];
            P2[1] = dataread2[3];
            back2[0] = dataread2[0];
            back2[1] = dataread2[1];
         }
         n2 = 2*j;
         area2 = area(P2[0],P2[1],back2[0],back2[1]);
         if (homogeneity)
            G2 = gradient(P2[0],P2[1],back2[0],back2[1],grad_img);
         if (P2[0] > back2[0])
            n2++;

         if (use_corners) {
            I = find_intersection(P1[0],P1[1],back1[0],back1[1],P2[0],P2[1],back2[0],back2[1]);
            C = find_corner(P1[0],P1[1],P2[0],P2[1],I[0],I[1], Ctable);
         } else {
            C = new double[2];
            C[0] = -1; C[1] = -1;
         }
         
         // calculate W1
         if (C[0] < 0) {
            M1[0] = (P1[0]+back1[0])/2;
            M2[0] = (P2[0]+back2[0])/2;
            M1[1] = (P1[1]+back1[1])/2;
            M2[1] = (P2[1]+back2[1])/2;
            if (LAMBDA > 0) {
               alpha = angle(back1[0], back1[1], P1[0], P1[1], P2[0], P2[1]);
               beta = angle(P1[0], P1[1], P2[0], P2[1], back2[0], back2[1]);
               if ((fabs(alpha) < (PI/2)) || (fabs(beta) < (PI/2))) continue;
               curv = curvature(M1[0],M1[1],P1[0],P1[1],P2[0],P2[1],M2[0],M2[1]);
               //cout << curv << endl;
               if (curv > 40) continue; //curv = 10; // avoid spikes in curvature
            } else curv = 0;
            w1 = pow(distance(P1[0],P1[1],P2[0],P2[1]),ALPHA) + LAMBDA*curv;
         } else {
            w1 = pow((distance(P1[0],P1[1],C[0],C[1]) + distance(C[0],C[1],P2[0],P2[1])), ALPHA);
         }

         // calculate W2
         gaparea = area(P1[0],P1[1],P2[0],P2[1]);
         if (C[0] < 0)
            carea = 0;
         else
            carea = triangle_area(P1[0],P1[1],C[0],C[1],P2[0],P2[1]);
         if (homogeneity)
            gapG = gradient(P1[0],P1[1],P2[0],P2[1], grad_img);
         if (P1[0] > P2[0]) {
            gaparea = -gaparea;
            carea = -carea;
            gapG = -gapG;
         }
         if (use_corners) delete I;

         // create dashed edge
         if (P1[0] > back1[0]) {

            dataout = new double[5];
            dataout[0] = 2*n1;

            if (P2[0] > back2[0]) {
               dataout[1] = 2*n2 + 1;
               dataout[2] = (area1 - area2)/2 + gaparea + carea - PSI*((G1 - G2)/2 + gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);

               dataout = new double[5];
               dataout[0] = 2*n1 + 1;
               dataout[1] = 2*n2;
               dataout[2] = (-area1 + area2)/2 - gaparea - carea - PSI*((-G1 + G2)/2 - gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);
               if (C[0] > 0) {
                  output << 2*n1 << " " << 2*n2+1 << " " << C[0] << " " << C[1] << endl;
                  output << 2*n1+1 << " " << 2*n2 << " " << C[0] << " " << C[1] << endl;
               }

            }
            else {
               dataout[1] = 2*n2;
               dataout[2] = (area1 + area2)/2 + gaparea + carea - PSI*((G1 + G2)/2 + gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);

               dataout = new double[5];
               dataout[0] = 2*n1 + 1;
               dataout[1] = 2*n2 + 1;
               dataout[2] = (-area1 - area2)/2 - gaparea - carea - PSI*((-G1 - G2)/2 - gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);
               if (C[0] > 0) {
                  output << 2*n1 << " " << 2*n2 << " " << C[0] << " " << C[1] << endl;
                  output << 2*n1+1 << " " << 2*n2+1 << " " << C[0] << " " << C[1] << endl;
               }

            }

         } else {

            dataout = new double[5];
            dataout[0] = 2*n1 + 1;

            if (P2[0] > back2[0]) {
               dataout[1] = 2*n2 + 1;
               dataout[2] = (-area1 - area2)/2 + gaparea + carea - PSI*((-G1 - G2)/2 + gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);

               dataout = new double[5];
               dataout[0] = 2*n1;
               dataout[1] = 2*n2;
               dataout[2] = (area1 + area2)/2 - gaparea - carea - PSI*((G1 + G2)/2 - gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);
               if (C[0] > 0) {
                  output << 2*n1+1 << " " << 2*n2+1 << " " << C[0] << " " << C[1] << endl;
                  output << 2*n1 << " " << 2*n2 << " " << C[0] << " " << C[1] << endl;
               }

            }
            else {
               dataout[1] = 2*n2;
               dataout[2] = (-area1 + area2)/2 + gaparea + carea - PSI*((-G1 + G2)/2 + gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);

               dataout = new double[5];
               dataout[0] = 2*n1;
               dataout[1] = 2*n2 + 1;
               dataout[2] = (area1 - area2)/2 - gaparea - carea - PSI*((G1 - G2)/2 - gapG);
               dataout[3] = w1;
               dataout[4] = 0; // dashed
               (*Gtable).push_back(dataout);
               if (C[0] > 0) {
                  output << 2*n1+1 << " " << 2*n2 << " " << C[0] << " " << C[1] << endl;
                  output << 2*n1 << " " << 2*n2+1 << " " << C[0] << " " << C[1] << endl;
               }

            }

         }

         delete[] C;

         count_dashed += 2;
         if (w1 > MAX_W1) MAX_W1 = w1;
         if (fabs(dataout[2]) > MAX_W2) MAX_W2 = fabs(dataout[2]);
      }
   }

   if (use_corners) output.close();

   cout << "Dashed edges: " << count_dashed << endl;

   dataout = new double[5];
   dataout[0] = MAX_W2;
   dataout[1] = MAX_W1;
   dataout[2] = 0;
   dataout[3] = 0;
   dataout[4] = 1;
   (*Gtable).push_back(dataout);

   cout << "MAX_W1: " << MAX_W1 << ", MAX_W2:" << MAX_W2 << endl;


   return Gtable;
}


double distance(double x1, double y1, double x2, double y2) {
   return sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
}

int minimal(double d1, double d2, double d3, double d4) {
   if (d1 <= d2 && d1 <= d3 && d1 <= d4){return 1;}
   if (d2 <= d1 && d2 <= d3 && d2 <= d4){return 2;}
   if (d3 <= d1 && d3 <= d2 && d3 <= d4){return 3;}
   if (d4 <= d1 && d4 <= d2 && d4 <= d3){return 4;}
   return 0;
}

double minimum(double d1, double d2, double d3, double d4) {
   if (d1 <= d2 && d1 <= d3 && d1 <= d4){return d1;}
   if (d2 <= d1 && d2 <= d3 && d2 <= d4){return d2;}