construct_graph.cpp

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

   if (d3 <= d1 && d3 <= d2 && d3 <= d4){return d3;}
   if (d4 <= d1 && d4 <= d2 && d4 <= d3){return d4;}
   return 0;
}

int sign(double d) {
   if (d == 0) return 0;
   return ((d>0)?1:-1);
}


int dround(double d) {
   return ((int)floor(d + 0.5));
}



double area(double P1x, double P1y, double P2x, double P2y) {

   return ( fabs(P1x-P2x)*( (P1y>P2y ? P1y : P2y) - fabs(P1y-P2y)/2) );

}

double triangle_area(double Ax, double Ay, double Bx, double By, double Cx, double Cy) {

   return ( fabs(-Bx*Ay + Cx*Ay + Ax*By - Cx*By - Ax*Cy + Bx*Cy)/2 );

}



bool quad_inside(double p1x,double p1y,double p2x,double p2y,double p3x,double p3y,double p4x,double p4y,double px,double py) {

   int s1, s2, s3;

   // first, check if p inside triangle p1p2p3

   s1 = sign(-p2x*p1y + px*p1y + p1x*p2y - px*p2y - p1x*py + p2x*py);
   s2 = sign(-p3x*p2y + px*p2y + p2x*p3y - px*p3y - p2x*py + p3x*py);
   s3 = sign(-p1x*p3y + px*p3y + p3x*p1y - px*p1y - p3x*py + p1x*py);

   if ((s1 == s2) && (s1 == s3)) { // point is inside
      return true;
   } else { // check second triangle, p3p4p1

      s1 = sign(-p3x*p1y + px*p1y + p1x*p3y - px*p3y - p1x*py + p3x*py);
      s2 = sign(-p4x*p3y + px*p3y + p3x*p4y - px*p4y - p3x*py + p4x*py);
      s3 = sign(-p1x*p4y + px*p4y + p4x*p1y - px*p1y - p4x*py + p1x*py);

      if ((s1 == s2) && (s1 == s3))
         return true;
      else
         return false;

   }

}

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

   double sum = 0, *line, m, b;
   int minx, maxx, maxy;

   if (P1x < P2x) {
      minx = dround(P1x);
      maxx = dround(P2x);
   } else {
      minx = dround(P2x);
      maxx = dround(P1x);
   }

   if (minx == maxx) return sum;

   m = (P1y-P2y)/(P1x-P2x);
   b = P1y - m*P1x;

   for (int i=minx; i<=maxx; i++) {

      maxy = dround(m*i + b);

      for (int j=0; j <= maxy; j++) {
         line = G[j];
         sum += line[i];
      }

   }

   return sum;

}

double* find_intersection(double P1x, double P1y, double P2x, double P2y, double P3x, double P3y, double P4x, double P4y) {

   double *I = new double[2];
   double m1,m2,b1,b2;

   if (P1x == P2x)
      m1 = 0;
   else {
      m1 = (P1y - P2y)/(P1x - P2x);
      b1 = P1y - m1*P1x;
   }
   if (P3x == P4x)
      m2 = 0;
   else {
      m2 = (P3y - P4y)/(P3x - P4x);
      b2 = P3y - m2*P3x;
   }

   if (P1x == P2x) {
      if (P3x == P4x) { //parallel, no intersection
         I[0] = -1;
         I[1] = -1;
      } else {
         I[0] = P1x;
         I[1] = m2*I[0] + b2;
      }
   } else {
      if (P3x == P4x) {
         I[0] = P3x;
         I[1] = m1*I[0] + b1;
      } else {
         if (m1 == m2) { //parallel, no intersection
            I[0] = -1;
            I[1] = -1;
         } else {
            I[0] = (b2-b1)/(m1-m2);
            I[1] = m1*I[0] + b1;
         }
      }
   }

   // discard I if P1,P2 and I are aligned
   if (((I[0] - P1x) < 0.01) && ((I[1] - P1y) < 0.01)) {
      I[0] = -1; I[1] = -1;
   } else if (((I[0] - P3x) < 0.01) && ((I[1] - P3y) < 0.01)) {
      I[0] = -1; I[1] = -1;
   }

   return I;
}

double* find_corner(double P1x, double P1y, double P2x, double P2y, double Ix, double Iy, vector<double*> Ctable) {

   double *C = new double[2], *p;
   double c1[2],c2[2],c3[2],c4[2];
   double m,b1,b2, d = 0;

   C[0] = -1;
   C[1] = -1;

   if (Ix == -1) {
      return C;
   }

   if (distance(P1x,P1y,Ix,Iy) > distance(P2x,P2y,Ix,Iy)) {
      c1[0] = P1x;
      c1[1] = P1y;
   } else {
      c1[0] = P2x;
      c1[1] = P2y;
   }

   double *inctr = triangle_incenter(P1x,P1y,Ix,Iy,P2x,P2y);

   if (Ix == inctr[0]) { //angle bisector is vertical
      c2[0] = 2*Ix - c1[0];
      c2[1] = c1[1];
   } else {
      m = (Iy-inctr[1])/(Ix-inctr[0]);
      if (m == 0) {
         c2[0] = c1[0];
         c2[1] = 2*Iy - c1[1];
      } else {
         b1 = c1[1] + c1[0]/m;
         b2 = Iy - m*Ix;
         c2[0] = (b2-b1)/(-1/m - m);
         c2[1] = m*c2[0] + b2;
      }
   }

   c3[0] = 2*Ix - c1[0];
   c3[1] = 2*Iy - c1[1];
   c4[0] = 2*Ix - c2[0];
   c4[1] = 2*Iy - c2[1];

/*if (P1x == 84 && P2x == 84) {
cout << "quad" << endl;
cout << c1[0] << " " << c1[1] << endl;
cout << c2[0] << " " << c2[1] << endl;
cout << c3[0] << " " << c3[1] << endl;
cout << c4[0] << " " << c4[1] << endl;
}*/

   //Search area is {c1,c2,c3,c4} 

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

      p = Ctable[i];

      if (p[2] > d) {
         if (quad_inside(c1[0],c1[1],c2[0],c2[1],c3[0],c3[1],c4[0],c4[1],p[0],p[1])) {
            C[0] = p[0];
            C[1] = p[1];
            d = p[2];
         }
      }
   }

   return C;

}


double* triangle_incenter(double Ax,double Ay,double Bx,double By,double Cx,double Cy) {

   double a,b,c, perimeter;
   double *center = new double[2];

   a = distance(Ax,Ay,Bx,By);
   b = distance(Bx,By,Cx,Cy);
   c = distance(Cx,Cy,Ax,Ay);

   perimeter = a+b+c;

   if ( perimeter == 0) {
      center[0] = Ax;
      center[1] = Ay;
   } else {
      center[0] = (b*Ax + c*Bx + a*Cx)/perimeter;
      center[1] = (b*Ay + c*By + a*Cy)/perimeter;
   }

   return center;
}


double curvature(double P0x, double P0y, double P1x, double P1y, double P2x, double P2y, double P3x, double P3y) {

    double total = 0, delta, x, xx, y, yy, K;
    double interm[4][2];
    int n = 100;
    double length, u, x1, y1, x2, y2;

    delta = 1.0 / n;

    interm[0][0] = -(P0x) + (3 * P1x) - (3 * P2x) + P3x;
    interm[0][1] = -(P0y) + (3 * P1y) - (3 * P2y) + P3y;
    interm[1][0] = (3 * P0x) - (6 * P1x) + (3 * P2x);
    interm[1][1] = (3 * P0y) - (6 * P1y) + (3 * P2y);
    interm[2][0] = -(3 * P0x) + (3 * P1x);
    interm[2][1] = -(3 * P0y) + (3 * P1y);
    interm[3][0] = P0x;
    interm[3][1] = P0y;

    for (int i = 0; i < n; i++)
    {
        x = 3 * (delta * delta * i * i) * interm[0][0] + 2 * delta * i * interm[1][0] + interm[2][0];
        xx = 6 * delta * i * interm[0][0] + 2 * interm[1][0];
        y = 3 * (delta * delta * i * i) * interm[0][1] + 2 * delta * i * interm[1][1] + interm[2][1];
        yy = 6 * delta * i * interm[0][1] + 2 * interm[1][1];

        if ( (x == 0) && (y == 0) )
            K = 0;
        else
        {
            K = ( fabs(x * yy - xx * y) / pow((fabs(x * x + y * y)), 1.5) );
            K = K * K;
        }

        u = delta*i;
        x1 = interm[0][0] * u * u * u + interm[1][0] * u * u + interm[2][0] * u + interm[3][0];
        y1 = interm[0][1] * u * u * u + interm[1][1] * u * u + interm[2][1] * u + interm[3][1];
        u = delta*(i+1);
        x2 = interm[0][0] * u * u * u + interm[1][0] * u * u + interm[2][0] * u + interm[3][0];
        y2 = interm[0][1] * u * u * u + interm[1][1] * u * u + interm[2][1] * u + interm[3][1];

        length = distance(x1,y1,x2,y2);

        total += K * length;
    }

    return total;

}

double angle(double p1x, double p1y, double p2x, double p2y, double p3x, double p3y){
   double p1,p2, value;

   p1 = (p3x-p2x)*(p1x-p2x) + (p3y-p2y)*(p1y-p2y);
   p2 = (p3x-p2x)*(p1y-p2y) - (p3y-p2y)*(p1x-p2x);

   if (p1 == 0 && p2 == 0) return 0;

   value = atan2(p2,p1);

   return value;
}