Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

LineData.cc

Go to the documentation of this file.
00001 //-*-c++-*-
00002 #include <iostream>
00003 #include <math.h>
00004 #include <vector>
00005 #include <list>
00006 
00007 #include "SketchSpace.h"
00008 #include "Sketch.h"
00009 #include "Region.h"
00010 #include "visops.h"
00011 
00012 #include "ShapeSpace.h"
00013 #include "ShapeRoot.h"
00014 
00015 #include "PointData.h"
00016 #include "LineData.h"
00017 #include "ShapeLine.h"
00018 
00019 #include "Crew/MapBuilder.h"
00020 #include "VRmixin.h"
00021 
00022 #ifdef PLATFORM_APERIOS
00023 //! this is normally defined in <math.h>, but the OPEN-R cross-compiler isn't configured right
00024 template<typename T> static inline int signbit(T x) { return x<0; }
00025 #endif
00026 
00027 using namespace std;
00028 
00029 namespace DualCoding {
00030 
00031 DATASTUFF_CC(LineData);
00032 
00033 const float LineData::DEFAULT_MIN_LENGTH = 0.025f;
00034 //const float LineData::DEFAULT_MIN_LENGTH = 8.0f;
00035 
00036 const Point LineData::origin_pt = Point(0,0);
00037 
00038 LineData::LineData(ShapeSpace& _space, const Point &p1, orientation_t orient)
00039   : BaseData(_space,getStaticType()), end1_pt(p1), end2_pt(), 
00040     rho_norm(0), theta_norm(0), orientation(0), length(0) {
00041   int const width = space->getDualSpace().getWidth();
00042   int const height = space->getDualSpace().getHeight();
00043   // Use a large offset from p1 to p2 because SketchGUI must calculate
00044   // the line slope from p1/p2 coords transmitted as strings; we don't
00045   // transmit the orientation.  But don't use an offset so large that the line
00046   // goes off-screen.
00047   float p2x=0, p2y=0;
00048   if ( fabs(orient-M_PI/2) < 0.001 ) {
00049     p2x = p1.coordX();
00050     p2y = p1.coordY() > height/2 ? 0 : height-1;
00051   } else {
00052     float slope = tan(orient);
00053     float intcpt = p1.coordY() - p1.coordX()*slope;
00054     p2x = p1.coordX() >= width/2 ? 0 : width-1;
00055     p2y = p2x * slope + intcpt;
00056     if ( p2y > height-1 ) {
00057       p2y = height-1;
00058       p2x = (p2y-intcpt) / slope;
00059     } else if ( p2y < 0 ) {
00060       p2y = 0;
00061       p2x = -intcpt/slope;
00062     }
00063   }
00064   end2_pt = Point(p2x,p2y);
00065   end1_pt.setValid(false);
00066   end1_pt.setActive(false);
00067   end2_pt.setValid(false);
00068   end2_pt.setActive(false);
00069   update_derived_properties();
00070 }
00071 
00072 LineData::LineData(const LineData& other)
00073   : BaseData(other),
00074     end1_pt(other.end1_pt), end2_pt(other.end2_pt),
00075     rho_norm(other.rho_norm), theta_norm(other.theta_norm),
00076     orientation(other.orientation), length(other.length) {}
00077 
00078 Point LineData::getCentroid() const { return (end1Pt()+end2Pt())*0.5f; }
00079 
00080 void LineData::setInfinite(bool value) {
00081   end1_pt.setActive(!value);
00082   end2_pt.setActive(!value);
00083 }
00084 
00085 #define NORMPOINT_MATCH_DISTSQ 400
00086 #define LINE_MATCH_OVERLAP -20
00087 #define LINE_ORIENTATION_MATCH M_PI/12
00088 #define LINE_PERPDIST 30
00089 //#define LINE_PERPDIST 20
00090 
00091 bool LineData::isMatchFor(const ShapeRoot& other) const {
00092   if (!(isSameTypeAs(other) && isSameColorAs(other)))
00093     return false;
00094   else {
00095     const Shape<LineData>& other_ln = ShapeRootTypeConst(other,LineData);
00096     return isMatchFor(other_ln.getData());
00097   }
00098 }
00099 
00100 bool LineData::isMatchFor(const LineData& other_line) const {
00101   //float const px = rho_norm*cos(theta_norm), 
00102   //  py = rho_norm*sin(theta_norm);
00103   //float const qx = other_line.rho_norm*cos(other_line.theta_norm),
00104   //  qy = other_line.rho_norm*sin(other_line.theta_norm);
00105   AngPi theta_diff = float(theta_norm) - float(other_line.theta_norm);
00106   if ( (orientation_t)theta_diff > (orientation_t)M_PI/2 )
00107     theta_diff = (orientation_t)M_PI - theta_diff;
00108   //float normpointdistsq = (px-qx)*(px-qx) + (py-qy)*(py-qy);
00109    //cout << "id=" << getId() << " id2=" << other_line.getId() << " px=" << px << "  py=" << py << "  qx=" << qx << "  qy=" << qy << "  normpointdistsq=" << normpointdistsq 
00110    //  << " th=" << theta_norm << " th2=" << other_line.theta_norm << " theta_diff=" << theta_diff
00111    //  << "  perpdist=" << perpendicularDistanceFrom(other_line.getCentroid())
00112    //  << "  isoverlapped=" << isOverlappedWith(other_line,LINE_MATCH_OVERLAP) << endl;
00113   float const linePerpDist =
00114     space->getRefFrameType() == camcentric ? 20 : 100;
00115   AngPi const lineOrientationMatch =
00116     space->getRefFrameType() == camcentric ? M_PI/6 : M_PI/4;
00117   //std::cout << "theta_diff: " << theta_diff 
00118   //  << ", perpendicularDistance1: " << perpendicularDistanceFrom(other_line.getCentroid()) 
00119   //  << ", perpendicularDistance2: " << other_line.perpendicularDistanceFrom(getCentroid()) 
00120   //  << ", isOverlappedWith: " << isOverlappedWith(other_line,LINE_MATCH_OVERLAP) << std::endl;
00121   //float perp1 = perpendicularDistanceFrom(other_line.getCentroid());
00122   //float perp2 = other_line.perpendicularDistanceFrom(getCentroid());
00123   return //normpointdistsq < NORMPOINT_MATCH_DISTSQ  // *** DST hack
00124     theta_diff < lineOrientationMatch
00125     && (perpendicularDistanceFrom(other_line.getCentroid()) < linePerpDist || other_line.perpendicularDistanceFrom(getCentroid()) < linePerpDist)  // was 100
00126     && isOverlappedWith(other_line,LINE_MATCH_OVERLAP);
00127 }
00128 
00129 bool LineData::isOverlappedWith(const LineData& otherline, int amount) const {
00130   fmat::Column<3> thisVec = end2_pt.coords - end1_pt.coords;
00131   fmat::Column<2> v = fmat::SubMatrix<2,1>(thisVec);
00132   float theta = atan2(v[1],v[0]);
00133   fmat::Matrix<2,2> rot = fmat::rotation2D(-theta);
00134   v = rot * v;
00135   fmat::Column<3> otherPt1 = otherline.end1_pt.coords - end1_pt.coords;
00136   fmat::Column<3> otherPt2 = otherline.end2_pt.coords - end1_pt.coords;
00137   fmat::Column<2> op1 = rot * fmat::SubMatrix<2,1>(otherPt1);
00138   fmat::Column<2> op2 = rot * fmat::SubMatrix<2,1>(otherPt2);
00139   float thisMin = min(0.f, v[0]);
00140   float thisMax = max(0.f, v[0]);
00141   float otherMin = min(op1[0], op2[0]);
00142   float otherMax = max(op1[0], op2[0]);
00143   return (thisMin < otherMin) ? (otherMin+amount <= thisMax) : (thisMin+amount <= otherMax);
00144 }
00145 
00146 void LineData::mergeWith(const ShapeRoot& other) {
00147   const Shape<LineData>& other_line = ShapeRootTypeConst(other,LineData);
00148   if (other_line->confidence <= 0)
00149     return;
00150   //const int other_conf = other_line->confidence;
00151   //confidence += other_conf;
00152   // Favor longer lines because they're likely to be more reliable
00153   float const totalLengthSq = length * length + other_line->length * other_line->length;
00154   EndPoint new_end1_pt = (end1_pt*length*length/totalLengthSq + other_line->end1Pt()*other_line->length*other_line->length/totalLengthSq);
00155   EndPoint new_end2_pt = (end2_pt*length*length/totalLengthSq + other_line->end2Pt()*other_line->length*other_line->length/totalLengthSq);
00156   new_end1_pt.valid = end1_pt.isValid() && other_line->end1Pt().isValid();
00157   new_end2_pt.valid = end2_pt.isValid() && other_line->end2Pt().isValid();
00158   end1_pt = new_end1_pt;
00159   end2_pt = new_end2_pt;
00160   update_derived_properties();
00161 }
00162 
00163 bool LineData::isValidUpdate(coordinate_t c1_cur, coordinate_t c2_cur, coordinate_t c1_new, coordinate_t c2_new) {
00164   const float c1_noise = 10 + std::abs(c1_cur+c1_new) / 20.f; // allow larger error for shapes far from the robot
00165   const float c2_noise = 10 + std::abs(c2_cur+c2_new) / 20.f;
00166   return (c1_new-c1_noise < c1_cur && c2_cur < c2_new+c2_noise);
00167 }
00168 
00169 
00170 //! Update a line in the local map with info from a matching line in the ground space.
00171 bool LineData::updateParams(const ShapeRoot& ground_root, bool force) {
00172   const Shape<LineData>& ground_line = ShapeRootTypeConst(ground_root,LineData);
00173   return updateParams(ground_line.getData(), force);
00174 }
00175 
00176 //! Update a line in the local map with info from a matching line in the ground space.
00177   bool LineData::updateParams(const LineData &ground_line, bool force) {
00178   //  cout << "updating local Line " << getId() << " with data from ground line " << ground_line->getId() << ":" << endl;
00179   //  ground_line->printEnds();
00180   //  cout << "Update from " << endl;
00181   //  printEnds();
00182   
00183   const coordinate_t c1_cur = firstPtCoord();
00184   const coordinate_t c2_cur = secondPtCoord();
00185 
00186   Point _end1_pt = firstPt();
00187   Point _end2_pt = secondPt();
00188   
00189   updateLinePt(firstPt(), firstPtCoord(), firstPt(ground_line), firstPtCoord(ground_line), -1);
00190   updateLinePt(secondPt(), secondPtCoord(), secondPt(ground_line), secondPtCoord(ground_line), +1);
00191   //  cout << "to" << endl;
00192   //  printEnds();
00193   
00194   const coordinate_t c1_new = firstPtCoord();
00195   const coordinate_t c2_new = secondPtCoord();
00196   
00197   if (isValidUpdate(c1_cur, c2_cur, c1_new, c2_new) || force){
00198     //    cout << "was accepted, line updated" << endl;
00199     //    ++nsamples;
00200     update_derived_properties();
00201     return true;
00202   }
00203   
00204   //  cout << "was denied, line not updated" << endl;
00205   setEndPts(_end1_pt, _end2_pt);
00206   return false;
00207 }
00208 
00209 void LineData::updateLinePt(EndPoint& localPt, coordinate_t local_coord,
00210                             const EndPoint& groundPt, coordinate_t ground_coord,
00211                             int sign) {
00212   if ( groundPt.isValid() ) {
00213     if ( localPt.isValid() )
00214       localPt.updateParams(groundPt);
00215     else
00216       localPt = groundPt;
00217   }
00218   else if ( (ground_coord - local_coord)*sign > 0 )
00219     localPt = groundPt;
00220 }
00221 
00222 bool LineData::isAdmissible() const {
00223   if (end1Pt().isValid() && end2Pt().isValid())
00224     return length >= VRmixin::mapBuilder->curReq->minLineLength;
00225   else
00226     return length >= VRmixin::mapBuilder->curReq->minRayLength;
00227 }
00228 
00229 //! Print information about this shape. (Virtual in BaseData.)
00230 void LineData::printParams() const {
00231   cout << "Type = " << getTypeName() << "  ID=" << getId() << "  ParentID=" << getParentId() << endl;
00232   cout << "  end1{" << end1Pt().coordX() << ", " << end1Pt().coordY()  << "}"
00233        << " active=" << end1Pt().isActive()
00234        << " valid=" << end1Pt().isValid() << endl;
00235   
00236   cout << "  end2{" << end2Pt().coordX() << ", " << end2Pt().coordY() <<  "}"
00237        << " active=" << end2Pt().isActive()
00238        << " valid=" << end2Pt().isValid() << std::endl;
00239   
00240   cout << "  rho_norm=" << rho_norm
00241        << ", theta_norm=" << theta_norm
00242        << ", orientation=" << getOrientation() 
00243        << ", length=" << getLength() << endl;
00244   
00245   printf("  color = %d %d %d\n",getColor().red,getColor().green,getColor().blue);
00246   
00247   cout << "  mobile=" << getMobile()
00248        << ", viewable=" << isViewable() << endl;
00249   
00250   vector<float> abc = lineEquation_abc();
00251   printf("  equ = %f %f %f\n",abc[0],abc[1],abc[2]);
00252 }
00253 
00254 void LineData::printEnds() const {
00255   cout << "  end1{" << end1Pt().coordX() << ", " << end1Pt().coordY()  << "}";
00256   cout << "  active=" << end1Pt().isActive() << ", valid=" << end1Pt().isValid() << endl;
00257   cout << "  end2{" << end2Pt().coordX() << ", " << end2Pt().coordY() <<  "}";
00258   cout << "  active=" << end2Pt().isActive() << ", valid=" << end2Pt().isValid() << endl;
00259   // cout << endl;
00260 }
00261 
00262 
00263 
00264 // *** Transformations. *** //
00265 
00266 //! Apply a transformation to this shape.
00267 void LineData::applyTransform(const fmat::Transform& Tmat, const ReferenceFrameType_t newref) {
00268   end1Pt().applyTransform(Tmat,newref);
00269   end2Pt().applyTransform(Tmat,newref);
00270   update_derived_properties();
00271 }
00272 
00273 void LineData::projectToGround(const fmat::Transform& camToBase, const PlaneEquation& groundplane) {
00274   end1Pt().projectToGround(camToBase,groundplane);
00275   end2Pt().projectToGround(camToBase,groundplane);
00276   update_derived_properties();
00277 }
00278 
00279 
00280 // *** Logical endpoints *** //
00281 
00282 EndPoint& LineData::leftPt() { return end1Pt().isLeftOf(end2Pt()) ? end1_pt : end2_pt; }
00283 const EndPoint& LineData::leftPt() const { return end1Pt().isLeftOf(end2Pt()) ? end1_pt : end2_pt; }
00284 EndPoint& LineData::rightPt() { return end1Pt().isLeftOf(end2Pt()) ? end2_pt : end1_pt; }
00285 const EndPoint& LineData::rightPt() const { return end1Pt().isLeftOf(end2Pt()) ? end2_pt : end1_pt; }
00286 EndPoint& LineData::topPt() { return end1Pt().isAbove(end2Pt()) ? end1_pt : end2_pt; }
00287 const EndPoint& LineData::topPt() const { return end1Pt().isAbove(end2Pt()) ? end1_pt : end2_pt; }
00288 EndPoint& LineData::bottomPt() { return end1Pt().isAbove(end2Pt()) ? end2_pt : end1_pt; }
00289 const EndPoint& LineData::bottomPt() const { return end1Pt().isAbove(end2Pt()) ? end2_pt : end1_pt; }
00290 
00291 Shape<PointData> LineData::leftPtShape() {
00292   Shape<PointData> result(new PointData(*space, leftPt()));
00293   result->setName("leftPt");
00294   result->inheritFrom(*this);
00295   result->setViewable(false);
00296   return result;
00297 }
00298 
00299 Shape<PointData> LineData::rightPtShape() {
00300   Shape<PointData> result(new PointData(*space, rightPt()));
00301   result->setName("rightPt");
00302   result->inheritFrom(*this);
00303   result->setViewable(false);
00304   return result;
00305 }
00306 
00307 Shape<PointData> LineData::topPtShape() {
00308   Shape<PointData> result(new PointData(*space, topPt()));
00309   result->setName("topPt");
00310   result->inheritFrom(*this);
00311   result->setViewable(false);
00312   return result;
00313 }
00314 
00315 Shape<PointData> LineData::bottomPtShape() {
00316   Shape<PointData> result(new PointData(*space, bottomPt()));
00317   result->setName("bottomPt");
00318   result->inheritFrom(*this);
00319   result->setViewable(false);
00320   return result;
00321 }
00322 
00323 EndPoint& LineData::firstPt() {
00324   if ( isNotVertical() )
00325     if ( end1Pt().coordX() < end2Pt().coordX() )
00326       return end1Pt();
00327     else return end2Pt();
00328   else
00329     if ( end1Pt().coordY() < end2Pt().coordY() )
00330       return end1Pt();
00331     else return end2Pt();
00332 }
00333     
00334 const EndPoint& LineData::firstPt() const {
00335   if ( isNotVertical() )
00336     if ( end1Pt().coordX() < end2Pt().coordX() )
00337       return end1Pt();
00338     else return end2Pt();
00339   else
00340     if ( end1Pt().coordY() < end2Pt().coordY() )
00341       return end1Pt();
00342     else return end2Pt();
00343 }
00344     
00345 const EndPoint& LineData::firstPt(const LineData &otherline) const {
00346   if ( isNotVertical() )
00347     if ( otherline.end1Pt().coordX() < otherline.end2Pt().coordX() )
00348       return otherline.end1Pt();
00349     else return otherline.end2Pt();
00350   else
00351     if ( otherline.end1Pt().coordY() < otherline.end2Pt().coordY() )
00352       return otherline.end1Pt();
00353     else return otherline.end2Pt();
00354 }
00355     
00356 EndPoint& LineData::secondPt() {
00357   if ( isNotVertical() )
00358     if ( end1Pt().coordX() > end2Pt().coordX() )
00359       return end1Pt();
00360     else return end2Pt();
00361   else
00362     if ( end1Pt().coordY() > end2Pt().coordY() )
00363       return end1Pt();
00364     else return end2Pt();
00365 }
00366     
00367 const EndPoint& LineData::secondPt() const {
00368   if ( isNotVertical() )
00369     if ( end1Pt().coordX() > end2Pt().coordX() )
00370       return end1Pt();
00371     else return end2Pt();
00372   else
00373     if ( end1Pt().coordY() > end2Pt().coordY() )
00374       return end1Pt();
00375     else return end2Pt();
00376 }
00377     
00378 const EndPoint& LineData::secondPt(const LineData &otherline) const {
00379   if ( isNotVertical() )
00380     if ( otherline.end1Pt().coordX() > otherline.end2Pt().coordX() )
00381       return otherline.end1Pt();
00382     else return otherline.end2Pt();
00383   else
00384     if ( otherline.end1Pt().coordY() > otherline.end2Pt().coordY() )
00385       return otherline.end1Pt();
00386     else return otherline.end2Pt();
00387 }
00388     
00389 Shape<PointData> LineData::firstPtShape() {
00390   Shape<PointData> result(new PointData(*space, firstPt()));
00391   result->setName("firstPt");
00392   result->inheritFrom(*this);
00393   result->setViewable(false);
00394   return result;
00395 }
00396 
00397 Shape<PointData> LineData::secondPtShape() {
00398   Shape<PointData> result(new PointData(*space, secondPt()));
00399   result->setName("secondPt");
00400   result->inheritFrom(*this);
00401   result->setViewable(false);
00402   return result;
00403 }
00404 
00405 coordinate_t LineData::firstPtCoord() const {
00406   return  isNotVertical() ? firstPt().coordX() : firstPt().coordY();
00407 }
00408 
00409 coordinate_t LineData::firstPtCoord(const LineData &otherline) const {
00410   return  isNotVertical() ?
00411     firstPt(otherline).coordX() :
00412     firstPt(otherline).coordY();
00413 }
00414 
00415 coordinate_t LineData::secondPtCoord() const {
00416   return  isNotVertical() ? secondPt().coordX() : secondPt().coordY();
00417 }
00418 
00419 coordinate_t LineData::secondPtCoord(const LineData &otherline) const {
00420   return  isNotVertical() ?
00421     secondPt(otherline).coordX() :
00422     secondPt(otherline).coordY();
00423 }
00424 
00425 
00426 // *** Functions to set endpoints. *** //
00427 
00428 void LineData::setEndPts(const EndPoint& _end1_pt, const EndPoint& _end2_pt) {
00429   end1_pt.setCoords(_end1_pt);
00430   end1_pt.setActive(_end1_pt.isActive());
00431   end1_pt.setValid(_end1_pt.isValid());
00432   end1_pt.setNumUpdates(_end1_pt.numUpdates());
00433 
00434   end2_pt.setCoords(_end2_pt);
00435   end2_pt.setActive(_end2_pt.isActive());
00436   end2_pt.setValid(_end2_pt.isValid());
00437   end2_pt.setNumUpdates(_end2_pt.numUpdates());
00438 
00439   update_derived_properties();
00440 }
00441 
00442 
00443 // *** Properties functions. *** //
00444 
00445 std::pair<float,float> LineData::lineEquation_mb() const {
00446   float m;
00447   if ((fabs(end2Pt().coordX() - end1Pt().coordX()) * BIG_SLOPE)
00448       <= fabs(end2Pt().coordY() - end2Pt().coordY()))
00449     m = BIG_SLOPE;
00450   else 
00451     m = (end2Pt().coordY() - end1Pt().coordY())/(end2Pt().coordX() - end1Pt().coordX());
00452   float b = end1Pt().coordY() - m * end1Pt().coordX();
00453   return pair<float,float>(m,b);
00454 }
00455 
00456 
00457 //! Determine parameters a, b, c, d satisfying the equation ax + bz = c.
00458 std::vector<float> LineData::lineEquation_abc_xz() const {
00459   float dx = end2Pt().coordX() - end1Pt().coordX();
00460   float dz = end2Pt().coordZ() - end1Pt().coordZ();
00461 
00462   std::vector<float> abc(3,1);
00463   float& a = abc[0];
00464   float& b = abc[1];
00465   float& c = abc[2];
00466 
00467   // If vertical...b = 0
00468   if((dx == 0)
00469      || (dz/dx > BIG_SLOPE)) {
00470     a = 1;
00471     b = 0;
00472     c = end1Pt().coordX();
00473   }
00474   
00475   // If horizontal...a = 0
00476   else if((dz == 0)
00477     || (dx/dz > BIG_SLOPE)) {
00478     a = 0;
00479     b = 1;
00480     c = end1Pt().coordZ();
00481   }
00482   
00483   // If not horizontal or vertical...a = 1.0
00484   else {
00485     a = 1;
00486     b = (end1Pt().coordX() - end2Pt().coordX()) 
00487       / (end2Pt().coordZ() - end1Pt().coordZ());
00488     c = end1Pt().coordX() + b*end1Pt().coordZ();
00489   }
00490   
00491   return(abc);
00492 
00493 }
00494 
00495 //! Determine parameters a, b, c satisfying the equation ax + by = c.
00496 std::vector<float> LineData::lineEquation_abc() const {
00497   float dx = end2Pt().coordX() - end1Pt().coordX();
00498   float dy = end2Pt().coordY() - end1Pt().coordY();
00499   
00500   std::vector<float> abc(3,1);
00501   float& a = abc[0];
00502   float& b = abc[1];
00503   float& c = abc[2];
00504 
00505   // If vertical...b = 0
00506   if( std::abs(dx) < 1.0e-6f || Slope(dy/dx) > BIG_SLOPE) {
00507     a = 1;
00508     b = 0;
00509     c = end1Pt().coordX();
00510   }
00511   
00512   // If horizontal...a = 0
00513  else if ( std::abs(dy) < 1.0e-6f || Slope(dx/dy) > BIG_SLOPE ) {
00514     a = 0;
00515     b = 1;
00516     c = end1Pt().coordY();
00517   }
00518   
00519   // If not horizontal or vertical...a = 1.0
00520   else {
00521     a = 1;
00522     //    b = (end1Pt().coordX() - end2Pt().coordX()) 
00523     // / (end2Pt().coordY() - end1Pt().coordY());
00524     b = -dx / dy;
00525     c = end1Pt().coordX() + b*end1Pt().coordY();
00526   }
00527   
00528   return(abc);
00529 }
00530 
00531 
00532 // *** Functions to set values dealing with orientation. *** //
00533 
00534 void LineData::update_derived_properties() {
00535   rho_norm = perpendicularDistanceFrom(origin_pt);
00536   const vector<float> abc = lineEquation_abc();
00537   const float& a1 = abc[0];
00538   const float& b1 = abc[1];
00539   const float& c1 = abc[2];
00540   const float c1sign = (c1 >= 0) ? 1 : -1;
00541   theta_norm = atan2(b1*c1sign, a1*c1sign);
00542   orientation = theta_norm + AngPi((orientation_t)M_PI/2);
00543   length = end1Pt().distanceFrom(end2Pt());
00544   const ReferenceFrameType_t ref = getRefFrameType();
00545   end1_pt.setRefFrameType(ref);
00546   end2_pt.setRefFrameType(ref);
00547   deleteRendering();
00548 }
00549 
00550 bool LineData::isNotVertical() const {
00551   const AngPi threshold = (orientation_t)M_PI / 3;
00552   const AngPi orient = getOrientation();
00553   return (orient <= threshold) || (orient >= (orientation_t)M_PI - threshold);
00554 }
00555 
00556 /*
00557 bool LineData::isRoughlyPerpendicularTo(Shape<LineData>& other) {
00558   AngPi threshold = M_PI_4;
00559   AngPi orientation_diff = getOrientation() - other->getOrientation();
00560   if((orientation_diff >= threshold) && (orientation_diff < (M_PI-threshold)))
00561     return true;
00562   else
00563     return false;
00564 }
00565 
00566 bool LineData::isExactlyPerpendicularTo(Shape<LineData>& other) {
00567   AngPi orientation_diff = getOrientation() - other->getOrientation();
00568   return (orientation_diff == M_PI_2);
00569 }
00570 
00571 */
00572 
00573 //================================================================
00574 
00575 /* There is a standard formula to determine whether a point is to the
00576 left of a ray.  However, we're using lines, not rays (so there is no
00577 direction), and we mean "left of" from the viewer's point of view, not
00578 the ray's.  So this code is a little more complex.  We must decide on
00579 canonical start and end points for the line (ray) based on the
00580 comparison we're trying to make.  (To test this, draw two lines
00581 forming an "X" and make sure the code works correctly for both lines.)
00582 Also, we want this to work in both camcentric and ego/allocentric
00583 reference frames, which adds further complexity. -- DST */
00584 
00585 bool LineData::pointIsLeftOf(const Point& pt) const {
00586   const EndPoint &p1 = 
00587     ( getRefFrameType() == camcentric )
00588     ? ((end1_pt.coordY() < end2_pt.coordY()) ? end1_pt : end2_pt)
00589     : ((end1_pt.coordX() < end2_pt.coordX()) ? end1_pt : end2_pt);
00590   const EndPoint &p2 = (&p1 == &end1_pt) ? end2_pt : end1_pt;
00591   if ( (getRefFrameType()==camcentric) ? (p1.coordY() == p2.coordY()) : (p1.coordX() == p2.coordX()) )
00592     return false;
00593   else
00594     return ( (p2.coordX() - p1.coordX()) * (pt.coordY() - p1.coordY()) -
00595        (p2.coordY() - p1.coordY()) * (pt.coordX() - p1.coordX()) ) > 0;
00596 }
00597 
00598 bool LineData::pointIsRightOf(const Point& pt) const {
00599   const EndPoint &p1 = 
00600     ( getRefFrameType() == camcentric )
00601     ? ((end1_pt.coordY() < end2_pt.coordY()) ? end1_pt : end2_pt)
00602     : ((end1_pt.coordX() < end2_pt.coordX()) ? end1_pt : end2_pt);
00603   const EndPoint &p2 = (&p1 == &end1_pt) ? end2_pt : end1_pt;
00604   if ( (getRefFrameType()==camcentric) ? (p1.coordY() == p2.coordY()) : (p1.coordX() == p2.coordX()) )
00605     return false;
00606   else
00607     return ( (p2.coordX() - p1.coordX()) * (pt.coordY() - p1.coordY()) -
00608        (p2.coordY() - p1.coordY()) * (pt.coordX() - p1.coordX()) ) < 0;
00609 }
00610 
00611 bool LineData::pointIsAbove(const Point& pt) const {
00612   const EndPoint &p1 = 
00613     ( getRefFrameType() == camcentric )
00614     ? ((end1_pt.coordX() < end2_pt.coordX()) ? end1_pt : end2_pt)
00615     : ((end1_pt.coordY() < end2_pt.coordY()) ? end1_pt : end2_pt);
00616   const EndPoint &p2 = (&p1 == &end1_pt) ? end2_pt : end1_pt;
00617   if ( (getRefFrameType()==camcentric) ? (p1.coordX() == p2.coordX()) : (p1.coordY() == p2.coordY()) )
00618     return false;
00619   else
00620     return ( (p2.coordX() - p1.coordX()) * (pt.coordY() - p1.coordY()) -
00621        (p2.coordY() - p1.coordY()) * (pt.coordX() - p1.coordX()) ) > 0;
00622 }
00623 
00624 bool LineData::pointIsBelow(const Point& pt) const {
00625   const EndPoint &p1 = 
00626     ( getRefFrameType() == camcentric )
00627     ? ((end1_pt.coordX() < end2_pt.coordX()) ? end1_pt : end2_pt)
00628     : ((end1_pt.coordY() < end2_pt.coordY()) ? end1_pt : end2_pt);
00629   const EndPoint &p2 = (&p1 == &end1_pt) ? end2_pt : end1_pt;
00630   if ( (getRefFrameType()==camcentric) ? (p1.coordX() == p2.coordX()) : (p1.coordY() == p2.coordY()) )
00631     return false;
00632   else
00633     return ( (p2.coordX() - p1.coordX()) * (pt.coordY() - p1.coordY()) -
00634        (p2.coordY() - p1.coordY()) * (pt.coordX() - p1.coordX()) ) < 0;
00635 }
00636 
00637 // *** These functions are true based on line length. *** //
00638 
00639 bool LineData::isLongerThan(const Shape<LineData>& other) const {
00640   return length > other->length; }
00641 
00642 bool LineData::isLongerThan(float ref_length) const {
00643   return length > ref_length; }
00644 
00645 bool LineData::isShorterThan(const Shape<LineData>& other) const {
00646   return length < other->length; }
00647 
00648 bool LineData::isShorterThan(float ref_length) const {
00649   return length < ref_length; }
00650 
00651 bool LineData::isBetween(const Point &p, const LineData &other) const {
00652   if (getOrientation() == other.getOrientation()) { // parallel lines
00653     float dl = perpendicularDistanceFrom(other.end1Pt()); // distance between the lines
00654     return (perpendicularDistanceFrom(p) <= dl && other.perpendicularDistanceFrom(p) <= dl);
00655   }
00656   else {
00657     bool b;
00658     const LineData p_line (*space,  p,  // line from intersection of this and other to p
00659          intersectionWithLine(other, b, b));
00660     const AngPi theta_pline = p_line.getOrientation();
00661     const AngPi theta_greater = 
00662       (getOrientation() > other.getOrientation()) ? getOrientation() : other.getOrientation();
00663     const AngPi theta_smaller = 
00664       (getOrientation() < other.getOrientation()) ? getOrientation() : other.getOrientation();
00665     if (theta_greater - theta_smaller > M_PI/2)
00666       return (theta_pline >= theta_greater || theta_pline <= theta_smaller);
00667     else
00668       return (theta_pline <= theta_greater && theta_pline >= theta_smaller);
00669   }
00670 }
00671 
00672   
00673 // ***  Check intersection. *** //
00674 
00675 bool
00676 LineData::intersectsLine(const Shape<LineData>& other) const {
00677   return intersectsLine(other.getData());
00678 }
00679 
00680 bool
00681 LineData::intersectsLine(const LineData& other) const {
00682   // Calculate F(x,y) = 0 for this line (x1,y1) to (x2,y2).
00683   pair<float,float> F = lineEquation_mb();
00684   
00685   // Calculate G(x,y) = 0 for L (x3,y3) to (x4,y4).
00686   pair<float,float> G = other.lineEquation_mb();
00687   
00688   // NOTE: These tests are assumed to be taking place in the space of
00689   // "this" line.  Therefore, the limits of line extent (for lines
00690   // with inactive endpoints) are calculated in the space of "this"
00691   // line.
00692   
00693   // JJW *** YOU NEED TO TAKE ACCOUNT OF END POINTS BEING TURNED OFF.
00694   
00695   //  float end1x, end1y, end2x, end2y, other_end1x, other_end1y, other_end2x, other_end2y;
00696   //  if(end1Pt().isActive()) {
00697   //    end1x = end1Pt().coordX();
00698   //    end1y = end1Pt().coordY();
00699   //  } else {
00700   //    end1x = 
00701   
00702   
00703   // TEST 1
00704   
00705   // Calculate r3 = F(x3,y3)
00706   float r3 = F.first * other.end1Pt().coordX() + F.second - other.end1Pt().coordY();
00707   
00708   // Calculate r4 = F(x4,y4)
00709   float r4 = F.first * other.end2Pt().coordX() + F.second - other.end2Pt().coordY();
00710   
00711   // If r3 != 0...
00712   // ...AND r4 != 0...
00713   // ...AND signbit(r3) == signbit(r4)...
00714   // ...THEN the lines do not intersect.
00715   
00716   if((r3 != 0)
00717      && (r4 != 0)
00718      && (signbit(r3) == signbit(r4))
00719      )
00720     return false;
00721   
00722   
00723   // TEST 2
00724   
00725   // Calculate r1 = G(x1,y1)
00726   float r1 = G.first * end1Pt().coordX() + G.second - end1Pt().coordY();
00727   
00728   // Calculate r2 = G(x2,y2)
00729   float r2 = G.first * end2Pt().coordX() + G.second - end2Pt().coordY(); 
00730   
00731   // If r1 != 0...
00732   // ...AND r2 != 0...
00733   // ...AND signbit(r1) == signbit(r2)...
00734   // ...THEN the lines do not intersect.
00735   
00736   if((r1 != 0)
00737      && (r2 != 0)
00738      && (signbit(r1) == signbit(r2))
00739      )
00740     return false;
00741   
00742   // Otherwise, the lines DO intersect.
00743   return true;
00744 }
00745 
00746 
00747 Point LineData::intersectionWithLine(const Shape<LineData>& other,
00748              bool& intersection_on_this,
00749              bool& intersection_on_other) const { 
00750   return intersectionWithLine(other.getData(), intersection_on_this,intersection_on_other); 
00751 }
00752 
00753 Point
00754 LineData::intersectionWithLine(const LineData& other, 
00755              bool& intersection_on_this, bool& intersection_on_other) const {
00756 
00757   // Based upon algorithm written by Paul Bourke, April 1989.
00758   // http://astronomy.swin.edu/~pbourke/geometry/lineline2d/  
00759   // Accessed July 20th 2004
00760   
00761   // Points 1 and 2 are on "this" line. Points 3 and 4 define the "other" line.
00762   float x1, x2, x3, x4, y1, y2, y3, y4;
00763   x1 = end1Pt().coordX();
00764   x2 = end2Pt().coordX();
00765   x3 = other.end1Pt().coordX();
00766   x4 = other.end2Pt().coordX();
00767   y1 = end1Pt().coordY();
00768   y2 = end2Pt().coordY();
00769   y3 = other.end1Pt().coordY();
00770   y4 = other.end2Pt().coordY();
00771   
00772   // x1 + u_a(x2-x1) = x3 + u_b(x4-x3)
00773   // y1 + u_a(y2-y1) = y3 + u_b(y4-y3)
00774   
00775   // u_a = (x4-x3)(y1-y3) - (y4-y3)(x1-x3)
00776   //       -------------------------------
00777   //       (y4-y3)(x2-x1) - (x4-x3)(y2-y1)
00778   
00779   // u_b = (x2-x1)(y1-y3) - (y2-y1)(x1-x3)
00780   //       -------------------------------
00781   //       (y4-y3)(x2-x1) - (x4-x3)(y2-y1)
00782   
00783   float denom = ((y4-y3)*(x2-x1) - (x4-x3)*(y2-y1));
00784   float u_a_numerator = ((x4-x3)*(y1-y3) - (y4-y3)*(x1-x3));
00785   float u_b_numerator = ((x2-x1)*(y1-y3) - (y2-y1)*(x1-x3));
00786   
00787   // If denominators of u_a and u_b are zero, then lines are parallel.
00788   if(denom == 0.0) {
00789     if (u_a_numerator == 0.0 && u_b_numerator == 0.0) {
00790       // cerr << "intersectionWithLine: lines are coincident!\n";
00791       // printParams();
00792       // other.printParams();
00793       return(end1Pt());
00794     }
00795     else {
00796        cout << x1 << " " << x2 << " " << x3 << " " << x4 << " "
00797       << y1 << " " << y2 << " " << y3 << " " << y4 << endl;
00798        cout << "this theta; " << getOrientation() << ", other theta: " << other.getOrientation() << endl;
00799        cerr << "ERROR in intersectionWithLine: lines are parallel!\n";
00800       return(Point(-9999.0f,-99999.0f));
00801     }
00802   }
00803   
00804   else {
00805     float u_a = u_a_numerator / denom;
00806     float u_b = u_b_numerator / denom;
00807     
00808     // If 0 <= u_a <=1  then intersection point is on segment a.
00809     if(0<=u_a && u_a<=1)
00810       intersection_on_this = true;
00811     else
00812       intersection_on_this = false;
00813     
00814     // If 0 <= u_b <=1  then intersection point is on segment b.
00815     if(0<=u_b && u_b<=1)
00816       intersection_on_other = true;
00817     else
00818       intersection_on_other = false;
00819     
00820     return(Point((x1+u_a*(x2-x1)),
00821      (y1+u_a*(y2-y1))));
00822   }
00823 }
00824 
00825 
00826 // *** Distance. *** //
00827 
00828 float LineData::perpendicularDistanceFrom(const Point& otherPt) const {
00829   // NOTE that this formula is rather slow, as it involves the
00830   // calculation of the line equation in addition to a square root here.
00831   
00832   // Using the formula from http://www.tpub.com/math2/8.htm
00833   vector<float> abc = lineEquation_abc();
00834   float& a = abc[0];
00835   float& b = abc[1];
00836   float& c = abc[2];
00837   
00838   // Distance...
00839   //    (x*A + y*B - C  )
00840   // abs(-------------  )
00841   //    (sqrt(a^2 + b^2)) 
00842   float d = fabs((a * otherPt.coordX() + b * otherPt.coordY() - c)/sqrt(a*a + b*b));
00843   //  cout << "abcd: " << a << ", " << b << ", " << c << ", " << d << endl;
00844   return(d);
00845 }
00846 
00847 
00848 
00849 // ==================================================
00850 // BEGIN LINE EXTRACTION CODE
00851 // ==================================================
00852 
00853 //#define EXTRACT_LINE_MIN_AREA 20
00854 #define EXTRACT_LINE_MIN_AREA 60
00855 
00856 //#define EXTRACT_LINE_MIN_SCORE 60
00857 //#define EXTRACT_LINE_MIN_SCORE 40
00858 #define EXTRACT_LINE_MIN_SCORE 15 // was 25
00859 #define BIG_SLOPE_CS 5000.0
00860 
00861 bool LineData::pointsOnPerimeterOfWindow(Sketch<bool> const& sketch, double t, double r, Point& pt1, Point& pt2) {
00862   const int width = sketch->getWidth(), height = sketch->getHeight();
00863   r -= RSIZE / 2;
00864   t *= M_PI/180;
00865   
00866   pt1.setCoords(cos(t), sin(t));
00867   pt1 *= r;
00868   t += M_PI/2.0;
00869   pt2.setCoords(cos(t), sin(t));
00870   pt2 *= 10;
00871   pt2 += pt1;
00872   
00873   //check if line is vertical, then we don't want to calculate slope
00874   if (fabs(pt2.coordX() - pt1.coordX()) <= 0.001) {
00875     pt2.setCoords(pt1.coordX(), height);
00876     return true;
00877   }
00878     
00879   double m = (pt2.coordY() - pt1.coordY()) / (pt2.coordX() - pt1.coordX());
00880   double b = pt1.coordY() - m*pt1.coordX();
00881   
00882   static vector<Point> accOfPoints;
00883   accOfPoints.clear();
00884   
00885   //If line intercepts y axis in the image frame
00886   if (b >= 0 && b <= height)
00887     accOfPoints.push_back(Point(0, b));
00888   
00889   //If the line intercepts the right of the image frame
00890   if ((m*width + b) >= 0 && (m*width + b) <= height)
00891     accOfPoints.push_back(Point( width, m*width + b));
00892     
00893   //if the line intercepts the bottom of the image frame
00894   if ((height-b)/m >= 0 && (height-b)/m <= width)
00895     accOfPoints.push_back(Point((height-b)/m, height));
00896     
00897   //if the line intercepts the top of the image frame
00898   if (-b/m >= 0 && -b/m <= width)
00899     accOfPoints.push_back(Point( -b/m, 0));
00900   
00901   //See if we have a valid intersection
00902   if (accOfPoints.size() != 2)
00903     return false;
00904     
00905   pt1 = accOfPoints[0];
00906   pt2 = accOfPoints[1];
00907   return true;
00908 }
00909 
00910 double LineData::getScore(int t, int r, int height, short hough[TSIZE][RSIZE], float ptsArray[TSIZE][RSIZE]) {
00911   int tally = hough[t][r];
00912   
00913   //scale by how long this line could be in the picture frame. Bias against horizontal lines.
00914   //Point pt1, pt2;
00915   //pointsOnPerimeterOfWindow(t, r, pt1, pt2);
00916   //double scale = min((height / (pt1-pt2).xyzNorm()), 1.0f);
00917   if (ptsArray[t][r] == 0)
00918     return 0;
00919   double scale = min((height / ptsArray[t][r]), 1.0f);
00920   if (tally * scale > 15)
00921     return tally * scale;
00922   return 0;
00923 }
00924 
00925 vector<Shape<LineData> > LineData::houghExtractLines(Sketch<bool> const& sketch, 
00926                                                      Sketch<bool> const& occluders,
00927                                                      const int num_lines) {
00928   const int width = sketch->getWidth(), height = sketch->getHeight();
00929   SketchSpace &SkS = sketch->getSpace();
00930   ShapeSpace &ShS = SkS.getDualSpace();
00931   NEW_SKETCH_N(skelDist,uint,visops::mdist(sketch));
00932   NEW_SKETCH_N(occ,bool,occluders)  
00933   short hough[TSIZE][RSIZE];    // bin counts for Hough table
00934   float ptsArray[TSIZE][RSIZE]; // max number of on-screen pixels for this bin (for score biasing)
00935   
00936   // Initialize matrices
00937   memset(hough, 0, TSIZE*RSIZE*sizeof(short));
00938   memset(ptsArray, 0, TSIZE*RSIZE*sizeof(float));
00939   Point pt1, pt2;
00940   for (int t = 0; t < TSIZE; t++)
00941     for (int r = 0; r < RSIZE; r++)
00942       if (pointsOnPerimeterOfWindow(sketch, t, r, pt1, pt2))
00943         ptsArray[t][r] = (pt1-pt2).xyNorm();
00944   
00945   // Clean out pixels on image edges
00946   NEW_SKETCH_N(edges, bool, visops::non_bounds(sketch, 2));
00947   
00948   // Populate the Hough table
00949   for (int x = 0; x < width; x++) {
00950     for (int y = 0; y < height; y++) {
00951       if ( edges(x,y) ) {
00952         for ( int t = 0; t < TSIZE; t++ ) {
00953           double theta = M_PI/180.0 * t;
00954           int r = int(x*cos(theta) + y*sin(theta) + RSIZE/2);
00955           hough[t][r]++;
00956         }
00957       }
00958     }
00959   }
00960   
00961   // Bias the scores for their chance of actually being a line
00962   for(int t = 0; t < TSIZE; t++) {
00963     for(int r = 0; r < RSIZE; r++) {
00964       hough[t][r] = (short)getScore(t, r, height, hough, ptsArray);
00965     }
00966   }
00967   
00968   vector<Shape<LineData> > lines_vec;
00969   while ( (int)lines_vec.size() < num_lines ) {
00970     int maxScore = 0;
00971     int maxT = 0;
00972     int maxR = 0;
00973     // Find highest score in Hough table
00974     for (int i = 0; i < TSIZE; i++)
00975       for(int j = 0; j < RSIZE; j++)
00976         if (hough[i][j] > maxScore) {
00977           maxScore = hough[i][j];
00978           maxT = i;
00979           maxR = j;
00980         }
00981     // std::cout << ". Score: " << maxScore << ", T: " << maxT << ", R: " << maxR << std::endl;
00982     if (maxScore < EXTRACT_LINE_MIN_SCORE) break;
00983     if (! pointsOnPerimeterOfWindow(sketch, maxT, maxR, pt1, pt2)) {
00984       std::cout << "Bad line " << endl;
00985       continue;
00986     }
00987     
00988     // Make some line segments
00989     float x_normpoint = (maxR - RSIZE/2)*cos(maxT * M_PI/180.0);
00990     float y_normpoint = (maxR - RSIZE/2)*sin(maxT * M_PI/180.0);
00991     // m is slope of this line (not the normal line)
00992     float m = std::max(std::min(std::tan((maxT * M_PI/180.0)+AngPi((orientation_t)M_PI/2)), BIG_SLOPE_CS), -BIG_SLOPE_CS);
00993     float b = y_normpoint - m*x_normpoint;
00994     Point pt(0, b);
00995     AngPi orientation = atan(m);
00996     int count = 0;
00997     for ( int startval = 0;; ) {
00998       if ( ++count > 6 ) break;  // safety check: can't have more than 6 colinear lines in an image
00999       Shape<LineData> extracted_line(ShS, pt, orientation);
01000       if ( std::abs(m) <= 1 )  // when |slope| <= 1, scan along x, else scan along y
01001         startval = extracted_line->scanHorizForEndPts(skelDist,occluders,m,b,startval);
01002       else
01003         startval = extracted_line->scanVertForEndPts(skelDist,occluders,m,b,startval);
01004       if ( startval < 0 ) {
01005         extracted_line.deleteShape();
01006         break;
01007       }
01008       extracted_line->setColor(sketch->getColor());
01009       extracted_line->setParentId(sketch->getViewableId());
01010       extracted_line->firstPt().checkValidity(width,height,beg_dist_thresh);
01011       extracted_line->secondPt().checkValidity(width,height,beg_dist_thresh);
01012       // Transform from SketchSpace coordinates to ShapeSpace coordinates
01013       SkS.applyTmatinv(extracted_line->firstPt().getCoords());
01014       SkS.applyTmatinv(extracted_line->secondPt().getCoords());
01015       extracted_line->update_derived_properties();
01016       // cout << "Extracted line from " << extracted_line->firstPt() << " to " << extracted_line->secondPt() << endl;
01017       bool matched = false;
01018       for ( vector<Shape<LineData> >::iterator ln_it = lines_vec.begin(); ln_it != lines_vec.end(); ln_it++ ) {
01019         Shape<LineData> &ln = *ln_it;
01020         if ( extracted_line->isMatchFor(ln) ) {
01021           matched = true;
01022           // earlier line will have a higher score, so only replace it if new one is at least 3 pixels longer
01023           if ( extracted_line->getLength() > 3 + ln->getLength() ) {
01024             //std::cout << "Extracted line " << extracted_line->firstPt() << "," << extracted_line->secondPt()
01025             //         << " len " << extracted_line->getLength()
01026             //         << " replacing " << ln->firstPt() << "," << ln->secondPt()
01027             //         << " len " << ln->getLength()
01028             //         << std::endl;
01029              ln_it->deleteShape();
01030              *ln_it = extracted_line;
01031           }
01032           else {
01033             extracted_line.deleteShape();
01034             //std::cout << "Extracted line " << extracted_line->firstPt() << "," << extracted_line->secondPt()
01035             //        << " len " << extracted_line->getLength()
01036             //        << " inferior to " << " " << ln->firstPt() << "," << ln->secondPt()
01037             //        << " len " << ln->getLength()
01038             //        << std::endl;
01039           }
01040           break;
01041         }
01042       }
01043       if (matched == false)
01044         lines_vec.push_back(extracted_line);
01045     }
01046     
01047     // Done extracting segments; clear this neighborhood of the Hough table
01048     const int tRange = 3; // was 10
01049     const int rRange = 3; // was 20
01050     for (int i = maxT-tRange; i <= maxT+tRange; i++ )
01051       for (int j = maxR-rRange; j <= maxR+rRange; j++ )
01052         if (i >=0 && i < TSIZE && j >=0 && j < RSIZE) {
01053           // std::cout << "clear: " << i << ", " << j << std::endl;
01054           hough[i][j] = 0;
01055         }
01056   } // end of loop to extract lines
01057   return lines_vec;
01058 }
01059 
01060 Shape<LineData> LineData::extractLine(Sketch<bool>& sketch) {
01061   //copy sketch as occluders
01062   NEW_SKETCH_N(occluders, bool, visops::copy(sketch));
01063   return houghExtractLines(sketch, occluders, 1)[0];
01064   
01065   //return houghExtractLines(sketch, sketch, 1)[0];
01066 }
01067 
01068 Shape<LineData> LineData::extractLine(Sketch<bool>& sketch, const Sketch<bool>& occlusions) {
01069   //copy sketch as occluders
01070   //NEW_SKETCH_N(occluders, bool, visops::copy(sketch));
01071   //return houghExtractLines(sketch, occluders, 1)[0];
01072   
01073   return houghExtractLines(sketch, occlusions, 1)[0];
01074 }
01075 
01076 vector<Shape<LineData> > LineData::extractLines(Sketch<bool> const& sketch, const int num_lines) {
01077   //copy sketch as occluders
01078   NEW_SKETCH_N(occluders, bool, visops::copy(sketch));
01079   return houghExtractLines(sketch, occluders, num_lines);
01080   
01081   //return houghExtractLines(sketch, sketch, num_lines);
01082 }
01083 
01084 vector<Shape<LineData> > LineData::extractLines(Sketch<bool> const& sketch, Sketch<bool> const& occluders, const int num_lines) {
01085   //copy sketch as occluders
01086   //NEW_SKETCH_N(occluders, bool, visops::copy(sketch));
01087   //return houghExtractLines(sketch, occluders, num_lines);
01088   
01089   return houghExtractLines(sketch, occluders, num_lines);
01090 }
01091 
01092 /*
01093 Shape<LineData> LineData::extractLine(Sketch<bool>& sketch) {
01094   NEW_SKETCH(fatmask,bool,visops::fillin(sketch,1,2,8));
01095   NEW_SKETCH(skel,bool,visops::skel(fatmask));
01096   return extractLine(skel, sketch);
01097 }
01098 
01099 
01100 Shape<LineData> LineData::oldExtractLine(Sketch<bool>& sketch) {
01101   NEW_SKETCH_N(fatmask,bool,visops::fillin(sketch,1,2,8));
01102   NEW_SKETCH_N(skel,bool,visops::skel(fatmask));
01103   return oldExtractLine(skel, sketch);
01104 }
01105 
01106 
01107 Shape<LineData> LineData::extractLine(Sketch<bool>& skeleton, 
01108               const Sketch<bool>& occlusions) {
01109   const int width = skeleton->getWidth(), height = skeleton->getHeight();
01110   SketchSpace &SkS = skeleton->getSpace();
01111   ShapeSpace &ShS = SkS.getDualSpace();
01112 
01113   // approximate largest skel region with line
01114   NEW_SKETCH_N(fatskel, bool, visops::fillin(skeleton,2,2,8));
01115   //NEW_SKETCH(fatskel, bool, visops::fillin(skeleton,1,2,8));
01116   NEW_SKETCH_N(labels,uint,visops::oldlabelcc(fatskel,visops::EightWayConnect));
01117   //NEW_SKETCH(labels,uint,visops::labelcc(skeleton, 10));
01118   list<Region> regions = Region::extractRegions(labels,EXTRACT_LINE_MIN_AREA);
01119   if ( regions.empty() ) {
01120     ShapeRoot invalid; // need to define a named variable to avoid warning on next line
01121     return ShapeRootType(invalid,LineData);
01122   }
01123   Region skelchunk = regions.front();  // regions not empty, so use the largest....
01124   AngPi orientation(skelchunk.findPrincipalAxisOrientation());
01125   Point centroid(skelchunk.findCentroid());
01126   if (! skelchunk.isContained(centroid, 3) ) {  //  region does not look like a straight line
01127     Shape<LineData> result(splitLine(ShS, skelchunk, skeleton, occlusions));
01128     if ( result.isValid() )
01129       return result;
01130     else
01131       cout << "*** LineData::splitLine returned an invalid shape!" << endl;
01132   }
01133 
01134   AngPi d_orientation = 0.005;
01135   int twists = 5;
01136   vector<LineData> lv_temp;
01137   NEW_SKETCH_N(skelDist,uint,visops::mdist(skeleton));
01138   for (int j = -twists; j <  twists; j++) {
01139     // Create a symbolic line object.
01140     LineData extracted_line(ShS, centroid, orientation + AngPi(j*d_orientation));
01141     // rho and theta describe the normal line.
01142     // normpoint is where the normal line intersects with this line.
01143     // normpoint may actually be off-image.
01144     float x_normpoint = extracted_line.getRhoNorm()*cos(extracted_line.getThetaNorm());
01145     float y_normpoint = extracted_line.getRhoNorm()*sin(extracted_line.getThetaNorm());
01146 
01147     // m is slope of this line (not the normal line)
01148     float m = std::max(std::min(std::tan(extracted_line.getThetaNorm()+AngPi((orientation_t)M_PI/2)), BIG_SLOPE), -BIG_SLOPE);
01149     float b = y_normpoint - m*x_normpoint;
01150     if ( std::abs(m) <= 1 )  // when |slope| <= 1, scan along x, else scan along y
01151       extracted_line.scanHorizForEndPts(skelDist,occlusions,m,b);
01152     else
01153       extracted_line.scanVertForEndPts(skelDist,occlusions,m,b);
01154     lv_temp.push_back(extracted_line);
01155     extracted_line.update_derived_properties();
01156   }
01157   stable_sort(lv_temp.begin(),lv_temp.end(),not2(LineData::LineDataLengthLessThan()));
01158   NEW_SHAPE(adjusted_line, LineData, new LineData(lv_temp[0]));
01159   adjusted_line->setName("LineData");
01160   adjusted_line->setColor(skeleton->getColor());
01161   adjusted_line->setParentId(skeleton->getViewableId());
01162   // Check to see if any endpoints are near any edge of the screen.
01163   // If they are, invalidate them, assuming that line continues beyond screen.
01164   adjusted_line->firstPt().checkValidity(width,height,beg_dist_thresh);
01165   adjusted_line->secondPt().checkValidity(width,height,beg_dist_thresh);
01166   // Transform from SketchSpace coordinates to ShapeSpace coordinates
01167   SkS.applyTmatinv(adjusted_line->firstPt().getCoords());
01168   SkS.applyTmatinv(adjusted_line->secondPt().getCoords());
01169   adjusted_line->update_derived_properties();
01170   return adjusted_line;
01171 }
01172 
01173 Shape<LineData> LineData::oldExtractLine(Sketch<bool>& skeleton, 
01174               const Sketch<bool>& occlusions) {
01175   const int width = skeleton->getWidth(), height = skeleton->getHeight();
01176   SketchSpace &SkS = skeleton->getSpace();
01177   ShapeSpace &ShS = SkS.getDualSpace();
01178 
01179   // approximate largest skel region with line
01180   NEW_SKETCH_N(labels,uint,visops::oldlabelcc(skeleton,visops::EightWayConnect));
01181   list<Region> regions = Region::extractRegions(labels,EXTRACT_LINE_MIN_AREA);
01182   if ( regions.empty() ) {
01183     ShapeRoot invalid; // need to define a named variable to avoid warning on next line
01184     return ShapeRootType(invalid,LineData);
01185   }
01186   Region skelchunk = regions.front();  // regions not empty, so use the largest....
01187   AngPi orientation(skelchunk.findPrincipalAxisOrientation());
01188   Point centroid(skelchunk.findCentroid());
01189   // cout << " orientation=" << orientation << "  centroid=" << centroid.getCoords() << endl;
01190   if (! skelchunk.isContained(centroid, 3) ) {   //  region does not look like a straight line
01191     Shape<LineData> result(oldSplitLine(ShS, skelchunk, skeleton, occlusions));
01192     if ( result.isValid() )
01193       return result;
01194   }
01195   
01196   // Create a symbolic line object.
01197   Shape<LineData> extracted_line(ShS, centroid, orientation);
01198   extracted_line->setColor(skeleton->getColor());
01199   extracted_line->setParentId(skeleton->getViewableId());
01200 
01201   // rho and theta describe the normal line.
01202   // normpoint is where the normal line intersects with this line.
01203   // normpoint may actually be off-image.
01204   float x_normpoint = extracted_line->rho_norm*cos(extracted_line->theta_norm);
01205   float y_normpoint = extracted_line->rho_norm*sin(extracted_line->theta_norm);
01206   
01207   // m is slope of this line (not the normal line)
01208   // float m = (y_normpoint != 0) ? -(x_normpoint/y_normpoint) : BIG_SLOPE;
01209   float m = std::max(std::min(std::tan(extracted_line->theta_norm+AngPi((orientation_t)M_PI/2)), BIG_SLOPE), -BIG_SLOPE);
01210   float b = y_normpoint - m*x_normpoint;
01211   // cout << "  x_normpoint=" << x_normpoint << "  y_normpoint=" << y_normpoint
01212   // << "  m=" << m << "  b=" << b <<std::endl;
01213   
01214   NEW_SKETCH_N(skelDist,uint,visops::mdist(skeleton));
01215   if ( std::abs(m) <= 1 )  // when |slope| <= 1, scan along x, else scan along y
01216     extracted_line->scanHorizForEndPts(skelDist,occlusions,m,b);
01217   else
01218     extracted_line->scanVertForEndPts(skelDist,occlusions,m,b);
01219 
01220   // Check to see if any endpoints are near any edge of the screen.
01221   // If they are, invalidate them, assuming that line continues beyond screen.
01222   extracted_line->end1_pt.checkValidity(width,height,beg_dist_thresh);
01223   extracted_line->end2_pt.checkValidity(width,height,beg_dist_thresh);
01224 
01225   // Transform from SketchSpace coordinates to ShapeSpace coordinates
01226   SkS.applyTmatinv(extracted_line->end1_pt.getCoords());
01227   SkS.applyTmatinv(extracted_line->end2_pt.getCoords());
01228   extracted_line->update_derived_properties();
01229   return extracted_line;
01230 }
01231 */
01232 
01233   // hack to split a non-straight region by kei
01234 Shape<LineData> LineData::splitLine(ShapeSpace &ShS, Region &skelchunk, 
01235         Sketch<bool> &skeleton, const Sketch<bool> &occlusions) {
01236   //cout << "this region is not straight, needs be split into two regions" << endl;
01237   Point bounds[4] = {skelchunk.findTopBoundPoint(), skelchunk.findLeftBoundPoint(),
01238          skelchunk.findRightBoundPoint(), skelchunk.findBotBoundPoint()};
01239   //  cout << "top(0): " << bounds[0] << endl;
01240   //  cout << "left(1): " << bounds[1] << endl;
01241   //  cout << "right(2): " << bounds[2] << endl;
01242   //  cout << "bottom(3): " << bounds[3] << endl;
01243   for (int i = 0; i < 4; i++) {
01244     for (int j = i+1; j < 4; j++) {
01245       if (bounds[i].distanceFrom(bounds[j]) > 20 && ! skelchunk.isContained((bounds[i]+bounds[j])/2.f, 3)) {
01246         //  cout << "[" << i << "," << j << "] form a line from which most distant point is where the region should split " << endl;
01247         LineData ln(ShS,bounds[i],bounds[j]);
01248         Point most_distant(skelchunk.mostDistantPtFrom(ln));
01249         ShS.getDualSpace().applyTmat(most_distant.getCoords());
01250         int const clear_dist = 15;  // was 10
01251         int x1 = max(0, int(most_distant.coordX() - clear_dist));
01252         int x2 = min(skeleton.width-1, int(most_distant.coordX() + clear_dist));
01253         int y1 = max(0, int(most_distant.coordY() - clear_dist));
01254         int y2 = min(skeleton.height-1, int(most_distant.coordY() + clear_dist));
01255         for (int x = x1; x <= x2; x++)
01256           for (int y = y1; y <= y2; y++)
01257             skeleton(x,y) = false;
01258         return extractLine(skeleton, occlusions);
01259       }
01260     }
01261   }
01262   ShapeRoot invalid; // need to define a named variable to avoid warning on next line
01263   return ShapeRootType(invalid,LineData);
01264 }
01265 
01266 /*
01267   // hack to split a non-straight region by kei
01268 Shape<LineData> LineData::oldSplitLine(ShapeSpace &ShS, Region &skelchunk, 
01269         Sketch<bool> &skeleton, const Sketch<bool> &occlusions) {
01270   //cout << "this region is not straight, needs be split into two regions" << endl;
01271   Point bounds[4] = {skelchunk.findTopBoundPoint(), skelchunk.findLeftBoundPoint(),
01272          skelchunk.findRightBoundPoint(), skelchunk.findBotBoundPoint()};
01273   //  cout << "top(0): " << bounds[0] << endl;
01274   //  cout << "left(1): " << bounds[1] << endl;
01275   //  cout << "right(2): " << bounds[2] << endl;
01276   //  cout << "bottom(3): " << bounds[3] << endl;
01277   for (int i = 0; i < 4; i++) {
01278     for (int j = i+1; j < 4; j++) {
01279       if (bounds[i].distanceFrom(bounds[j]) > 20 && ! skelchunk.isContained((bounds[i]+bounds[j])/2.f, 3)) {
01280   //  cout << "[" << i << "," << j << "] form a line from which most distant point is where the region should split " << endl;
01281   LineData ln(ShS,bounds[i],bounds[j]);
01282   PointData most_distant_pt(ShS, skelchunk.mostDistantPtFrom(ln));
01283   //  cout << "Point of split: " << most_distant_pt.getCentroid() << endl;
01284   const Sketch<bool>& point_rendering = most_distant_pt.getRendering();
01285   uint const clear_dist = 10;
01286   NEW_SKETCH_N(not_too_close, bool, visops::mdist(point_rendering) >= clear_dist);
01287   skeleton &= not_too_close;
01288   return oldExtractLine(skeleton, occlusions);
01289       }
01290     }
01291   }
01292   ShapeRoot invalid; // need to define a named variable to avoid warning on next line
01293   return ShapeRootType(invalid,LineData);
01294 }
01295 */
01296 int LineData::scanHorizForEndPts(const Sketch<uint>& skelDist,
01297                                  const Sketch<bool>& occlusions,
01298                                  float m, float b, int xstart) {
01299   // Scan along the infinite line, looking for a segment in the image.
01300   // Adjust the line's endpoints based on what we find.
01301   // Return stopping point so we can be called again for additional segments.
01302   bool on_line = false;  // true if tracing a line skeleton
01303   int curxstart = -1;    // start of current segment we're examining
01304   int possxstop = -1;    // end of current segment unless we bridge a gap
01305   int bestlength = -1;   // length of best segment seen so far
01306   int goodpixels;        // number of pixels that aren't gaps
01307   int const xstop = skelDist.height - line_min_length;
01308   while ( xstart < xstop ) {
01309     possxstop = -1;
01310     for (int x = xstart, y, dist=0; x < skelDist.width; x++) {
01311       y = int(m*x+b);    
01312       if (y < 0 || y >= skelDist.height) continue;  // make sure y is on-screen
01313       dist = skelDist(x,y);
01314       if (on_line == false) {   // not currently on a line segment
01315         if (dist <= beg_dist_thresh) { // start a new segment
01316           on_line = true;
01317           goodpixels = 0;
01318           curxstart = x - dist;
01319           // if new segment begins at an occluder, back up over the occluder
01320           int curystart;
01321           bool backupflag = false;
01322           while ( curxstart >= 0 && (curystart=int(m*curxstart+b)) >= 0 && curystart < skelDist.height ) {
01323             if ( occlusions(curxstart,curystart) || skelDist(curxstart,curystart) == 0 ) {
01324               --curxstart;
01325               backupflag = true;
01326             } else { // if we backed up past the occluder or segment, go one step forward
01327               curxstart += backupflag;
01328               break;
01329             }
01330           }
01331           if ( curxstart < 0) // occluder extended to left edge of image
01332             curxstart = 0;
01333         }
01334       }
01335       else {   // on_line == true:  currently on a line segment
01336         if ( dist <= end_dist_thresh || occlusions(x,y) ) {
01337           ++goodpixels;
01338           possxstop = x;
01339         }
01340         else if ( dist > extractorGapTolerance ) {
01341           // we're traversing a gap, and it just got too long
01342           on_line = false;
01343           bestlength = possxstop - curxstart;
01344           if ( bestlength >= line_min_length )
01345             break; // proceed if we've found a good segment
01346         }
01347       }
01348     }
01349     bestlength = possxstop - curxstart;
01350     // cout << "  vert bestlength=" << bestlength << "  goodpixels=" << goodpixels << endl;
01351     if ( bestlength <= 0 )  // could not find a candidate line; we're done
01352       return -1;
01353     if ( bestlength >= line_min_length && float(goodpixels)/bestlength >= 0.8 ) {  // found a good line; return it
01354       float y1 = m*curxstart + b;
01355       float y2 = m*possxstop + b;
01356       setEndPts(Point(curxstart,y1), Point(possxstop,y2));
01357       // cout << "m = " << m << "   b = " << b << endl;
01358       return possxstop + end_dist_thresh;
01359     }
01360     xstart = max(xstart+1,possxstop) + end_dist_thresh;
01361   }
01362   // ran off the end of the screen; we're done
01363   return -1;
01364 }
01365 
01366 int LineData::scanVertForEndPts(const Sketch<uint>& skelDist, 
01367                                 const Sketch<bool>& occlusions, 
01368                                 float m, float b, int ystart) {
01369   // Scan along the infinite line, looking for segments in the image.
01370   // Adjust the line's endpoints based on what we find.
01371   // Return stopping point so we can be called again for additional segments.
01372   bool on_line = false;  // true if currently tracing a line skeleton
01373   int curystart = -1;    // start of current segment we're examining
01374   int possystop = -1;    // end of current segment unless we bridge a gap
01375   int bestlength = -1;   // length of segment seen so far
01376   int goodpixels = 0;    // number of pixels that aren't gaps
01377   int const ystop = skelDist.height - line_min_length;
01378   while ( ystart < ystop ) {
01379     // scan for a candidate line
01380     possystop = -1;
01381     for (int x, y = ystart, dist=0; y < skelDist.height; y++) {
01382       x = int((y-b)/m);
01383       if (x < 0 || x >= skelDist.width) continue;  // make sure x is on-screen
01384       dist = skelDist(x,y);
01385       if (on_line == false) {   // not currently on a line segment
01386         if (dist <= beg_dist_thresh) { // start a new segment
01387           on_line = true;
01388           goodpixels = 0;
01389           curystart = y - dist;
01390           // if new segment begins at an occluder, back up over the occluder
01391           int curxstart;
01392           bool backupflag = false;
01393           while ( curystart >= 0 && (curxstart=int((curystart-b)/m)) >= 0 && curxstart < skelDist.width )
01394             if ( occlusions(curxstart,curystart) || skelDist(curxstart,curystart) == 0 ) {
01395               --curystart;
01396               backupflag = true;
01397             } else { // if we backed up past the occluder or segment, go one step forward
01398               curystart += backupflag;
01399               break;
01400             }
01401           if ( curystart < 0) // occluder extended to top edge of image
01402             curystart = 0;
01403         }
01404       }
01405       else {   // on_line == true:  currently on a line segment
01406         if ( dist <= end_dist_thresh || occlusions(x,y) ) {
01407           ++goodpixels;
01408           possystop = y;
01409         }
01410         else if ( dist > extractorGapTolerance ) {
01411           // we're traversing a gap, and it just got too long
01412           on_line = false;
01413           bestlength = possystop - curystart;
01414           if ( bestlength >= line_min_length )
01415             break; // proceed if we've found a good setment
01416         } 
01417       }
01418     }
01419     bestlength = possystop - curystart;
01420     // cout << "  vert bestlength=" << bestlength << "  goodpixels=" << goodpixels << endl;
01421     if ( bestlength <= 0 )  // could not find a candidate line; we're done
01422       return -1;
01423     if ( bestlength >= line_min_length && float(goodpixels)/bestlength >= 0.8 ) {  // found a good line; return it
01424       float x1 = (curystart-b)/m;
01425       float x2 = (possystop-b)/m;
01426       setEndPts(Point(x1,curystart), Point(x2,possystop));
01427       return possystop + end_dist_thresh;
01428     }
01429     // found a candidate but it didn't pan out; look some more
01430     // cout << " ystart=" << ystart << " curystart=" << curystart << " possystop=" << possystop;
01431     ystart = max(ystart+1,possystop) + end_dist_thresh;
01432     // cout << " new ystart=" << ystart << endl;
01433   }
01434   // ran off the end of the screen; we're done
01435   return -1;
01436 }
01437 
01438 //! Clear out pixels that are on or close to this line.
01439 void LineData::clearLine(Sketch<bool>& sketch) {
01440   const Sketch<bool>& line_rendering = getRendering();
01441   uint const clear_dist = 5;
01442   Sketch<bool> not_too_close = (visops::mdist(line_rendering) >= clear_dist);
01443   sketch &= not_too_close;
01444 } 
01445 
01446   
01447 
01448 // ==================================================
01449 // BEGIN LINE RENDERING CODE
01450 // ==================================================
01451 
01452 Sketch<bool>& LineData::getRendering() {
01453   if ( ! end1Pt().rendering_valid || ! end2Pt().rendering_valid )
01454     deleteRendering();
01455   if ( rendering_sketch == NULL )
01456     rendering_sketch = render();
01457   return *rendering_sketch;
01458 }
01459 
01460 //! Render line to SketchSpace and return a pointer to the sketch.
01461 //! This function does not link the Sketch<bool>* in the shape to the sketch returned.
01462 Sketch<bool>* LineData::render() const {
01463   SketchSpace &renderspace = space->getDualSpace();
01464   Sketch<bool>* draw_result = 
01465     new Sketch<bool>(renderspace, "render("+getName()+")");
01466   //  (*draw_result)->setParentId(getViewableId());
01467   //  (*draw_result)->setColor(getColor());
01468   (*draw_result)->inheritFrom(*this);
01469   *draw_result = 0;
01470   renderOnTo(draw_result);
01471   return draw_result;
01472 }  
01473 
01474 void LineData::renderOnTo(Sketch<bool>* draw_result) const {
01475   SketchSpace &renderspace = space->getDualSpace();
01476   int const width = renderspace.getWidth();
01477   int const height = renderspace.getHeight();
01478   float x1,y1,x2,y2;
01479   setDrawCoords(x1, y1, x2, y2, width, height);
01480 
01481   drawline2d(*draw_result, (int)x1, (int)y1, (int)x2, (int)y2);
01482   end1Pt().rendering_valid = true;
01483   end2Pt().rendering_valid = true;
01484 }   
01485   
01486 // This function will be called by both LineData and PolygonData renderers
01487 void LineData::setDrawCoords(float& x1,float& y1, float& x2, float& y2,
01488            const int width, const int height) const {
01489   EndPoint e1(end1Pt());
01490   EndPoint e2(end2Pt());
01491   //std::cout << "setDrawCoords: e1=" << e1 << " e2=" << e2 << std::endl;
01492   space->getDualSpace().applyTmat(e1.getCoords());
01493   space->getDualSpace().applyTmat(e2.getCoords());
01494   const EndPoint &left_point = e1.coordX() <= e2.coordX() ? e1 : e2;
01495   const EndPoint &right_point = e1.coordX() <= e2.coordX() ? e2 : e1;
01496   //std::cout << "setDrawCoords: left=" << left_point << " right=" << right_point << std::endl;
01497 
01498   // Check if horizontal
01499   if( left_point.coordY() == right_point.coordY()) { 
01500     if (!left_point.isActive())
01501       x1 = 0;
01502     else x1 = max(0.0f,min(width-1.0f,left_point.coordX()));
01503     
01504     if (!right_point.isActive())
01505       x2 = width-1;
01506     else x2 = max(0.0f,min(width-1.0f,right_point.coordX()));
01507     
01508     y1 = left_point.coordY();
01509     y2 = y1;
01510   } 
01511   else if ( left_point.coordX() == right_point.coordX()) {  // Check if vertical...
01512 
01513     const EndPoint &top_point =    e1.coordY() <= e2.coordY() ? e1 : e2;
01514     const EndPoint &bottom_point = e1.coordY() <= e2.coordY() ? e2 : e1;
01515 
01516   //std::cout << "top_point=" << top_point << " bottom_point=" << bottom_point << std::endl;
01517     if(!top_point.isActive())
01518       y1 = 0;
01519     else y1 = max(0.0f,min(height-1.0f,top_point.coordY()));
01520     
01521     if (!bottom_point.isActive())
01522       y2 = height-1;
01523     else y2 = max(0.f,min(height-1.0f,bottom_point.coordY()));
01524 
01525     x1 = left_point.coordX();
01526     x2 = x1;
01527   } 
01528 
01529   else {   // Neither horizontal nor vertical...
01530     float const m = (right_point.coordY()-left_point.coordY())/(right_point.coordX()-left_point.coordX());
01531     float const b = left_point.coordY() - m*left_point.coordX();
01532     
01533     // find image edge intersections
01534     int const i0x = (int)((0-b)/m);
01535     int const ihx = (int)(((height-1)-b)/m);
01536     int const i0y = (int)(m*0+b);
01537     int const iwy = (int)(m*(width-1)+b);
01538     
01539     // If left point is active, set starting x and y accordingly.
01540     if ( left_point.isActive() ) {
01541       x1 = left_point.coordX();
01542       y1 = left_point.coordY();
01543     } 
01544     
01545     // If endpoint 1 extends past image edge...
01546     else { 
01547       
01548       // intersects left edge
01549       if ( i0y >= 0 && i0y < height ) { 
01550   x1 = 0;
01551   y1 = i0y;
01552       } 
01553       
01554       // intersects top or bottom edge
01555       else { 
01556   
01557   // intersects top first
01558   if ( i0x < ihx ) { 
01559     x1 = i0x;
01560     y1 = 0;
01561   } 
01562   
01563   // intersects bottom first
01564   else { 
01565     x1 = ihx;
01566     y1 = height-1;
01567   }
01568       }
01569     }
01570     
01571     // If right point is active, set starting x and y accordingly.
01572     if ( right_point.isActive() ) {
01573       x2 = right_point.coordX();
01574       y2 = right_point.coordY();
01575     } 
01576     else { // endpoint extends to image edge
01577       if ( iwy >= 0 && iwy < height ) { // intersects right edge
01578   x2 = width-1;
01579   y2 = iwy;
01580       } 
01581       else { // intersects top or bottom edge
01582   if ( i0x > ihx ) { // intersects top last 
01583     x2 = i0x;
01584     y2 = 0;
01585   } 
01586   else { // intersects bottom last
01587     x2 = ihx;
01588     y2 = height-1;
01589   }
01590       }
01591     }
01592   }
01593 }
01594 
01595 void LineData::drawline2d(Sketch<bool>& canvas, int x0, int y0, int x1, int y1) {
01596   int width = canvas->getWidth();
01597   int height = canvas->getHeight();
01598   
01599   //std::cout << "LineData::drawline2D x0=" << x0 << " y0=" << y0 << " x1=" << x1 << " y1=" << y1 << std::endl;
01600 
01601   // code below from free Graphics Gems repository, graphicsgems.org
01602   int dx = x1-x0,  ax = abs(dx)<<1,  sx = signbit(dx) ? -1 : 1;
01603   int dy = y1-y0,  ay = abs(dy)<<1,  sy = signbit(dy) ? -1 : 1;
01604   int d, x=x0, y=y0;
01605   
01606   if ( ax > ay ) {    /* x dominant */
01607     //std::cout << "  x dominant\n";
01608     d = ay-(ax>>1);
01609     for (;;) {
01610       if (x >= 0 && y >= 0 && x < width && y < height)
01611   canvas(x,y) = true;
01612       
01613       if (x==x1)
01614   return;
01615       
01616       if ( d >= 0 ) {
01617   y += sy;
01618   d -= ax;
01619       }
01620       
01621       x += sx;
01622       d += ay;
01623     }
01624   }
01625   else {      /* y dominant */
01626     //std::cout << "  y dominant\n";
01627     d = ax-(ay>>1);
01628     for (;;) {
01629       if (x >= 0 && y >= 0 && x < width && y < height)
01630   canvas(x,y) = true;
01631       
01632       if ( y == y1 )
01633   return;
01634       
01635       if ( d >= 0 ) {
01636   x += sx;
01637   d -= ay;
01638       }
01639       
01640       y += sy;
01641       d += ax;
01642     }
01643   }
01644 }
01645 //}
01646 
01647 
01648 //****** NOT CURRENTLY USED *******
01649 std::vector<Shape<LineData> > LineData::houghTransform(const Sketch<bool>& fat, 
01650             const Sketch<bool>& skinny,
01651             const Sketch<bool>& occlusions,
01652             const size_t num_lines, int minLength)
01653 {
01654   if ( minLength == -1 )
01655     minLength = static_cast<int>(fat->getWidth() * DEFAULT_MIN_LENGTH);
01656   std::vector<Shape<LineData> > result;
01657   ShapeSpace& ShS = fat->getSpace().getDualSpace();
01658 
01659   const int width = fat->getWidth(), height = fat->getHeight();
01660   const int numR = 120, numTheta = 120;
01661   const int numRf = 40, numThetaf = 40;
01662   int hsize = numR*numTheta, hsizef = numRf * numThetaf;
01663   int htab[hsize], hfine[hsizef]; // Hough accumulator table
01664   float minTheta = 0, maxTheta = 2*(float)M_PI; //0.9999*M_PI;
01665   float minR = 0, maxR = std::sqrt((float)width*(float)width+(float)height*(float)height); //240.0;
01666   float thetaSpread = maxTheta - minTheta;
01667   float rSpread = maxR - minR;
01668   //float bestR = 0, bestTheta = 0;
01669 
01670   for (int i = 0; i < hsize; i++)
01671     htab[i] = 0; // zero accumulator table
01672 
01673   // build accumulator table
01674   float theta, r;
01675   int ridx;
01676   for (int x = 0; x < width; x++) {
01677     for (int y = 0; y < height; y++) {
01678       if (fat(x,y) == true) {
01679         for (int tidx = 0; tidx < numTheta ; tidx++) {
01680           theta = minTheta + tidx*thetaSpread/numTheta;
01681           r = (float)x*cos(theta)+(float)y*sin(theta);
01682           ridx = (int)((r-minR)*(float)numR/rSpread); // Is this right calc?
01683           if (ridx >= 0 && ridx < numR) // check ridx bounds
01684             htab[ridx+tidx*numR]++; 
01685         }
01686       }
01687     }
01688   }
01689 
01690 
01691   /*float lineLen = 500;
01692   for (int i = 0; i < numR*numTheta; i++) {
01693     if (htab[i] >= minCount)
01694       {
01695   const float curR = (i%numR)*rSpread/numR + minR;
01696   const float curTheta = (i/numR)*thetaSpread/numTheta + minTheta;
01697   const float x1 = curR*cos(curTheta), y1 = curR*sin(curTheta);
01698   const float x2 = x1+lineLen*cos(curTheta+M_PI/2), y2 = y1+lineLen*sin(curTheta+M_PI/2);
01699   const float x3 = x1-lineLen*cos(curTheta+M_PI/2), y3 = y1-lineLen*sin(curTheta+M_PI/2);
01700 
01701   Point e1(x3,y3), e2(x2,y2);
01702   Shape<LineData> line(ShS,e1,e2);
01703   line->end1Pt().setValid(false);
01704   line->end2Pt().setValid(false);
01705   result.push_back(line);
01706       }
01707       }*/
01708 
01709   // Finely explore the first n best lines in the rough table
01710   int max_val = -1, max_i = 0;
01711   while(result.size() < num_lines) {
01712 
01713     // find maximum value in accumulator table
01714     max_val = -1;
01715     max_i = 0;
01716     for (int i = 0; i < numR*numTheta; i++) {
01717       if (max_val < htab[i]) {
01718         max_val = htab[i];
01719         max_i = i;
01720       }
01721     }
01722 
01723     // Give up if the lines start getting too small
01724     if (max_val < minLength)
01725       break;
01726     
01727     float bestR = (max_i%numR)*rSpread/numR + minR;
01728     float bestTheta = (max_i/numR)*thetaSpread/numTheta + minTheta;
01729 
01730     // determine parameters for the fine hough iteration
01731     const float fthetaSpread = (float)M_PI/40.f;
01732     const float frSpread = maxR/40.f;
01733     float fminTheta = bestTheta-fthetaSpread/2.f;
01734     float fminR = bestR-frSpread/2.f;
01735 
01736     for (int i = 0; i < hsizef; i++)
01737       hfine[i] = 0; // zero fine table
01738     
01739     // build fine table
01740     float thetaf, rf;
01741     int ridxf;
01742     for (int x = 0; x < width; x++) {
01743       for (int y = 0; y < height; y++) {
01744         if (skinny(x,y) == true) {
01745           for (int tidx = 0; tidx < numThetaf ; tidx++) {
01746             thetaf = fminTheta + tidx*fthetaSpread/numThetaf;
01747             rf = (float)x*cos(thetaf)+(float)y*sin(thetaf);
01748             ridxf = (int)((rf-fminR)*(float)numRf/frSpread); // Is this right calc?
01749             if (ridxf >= 0 && ridxf < numRf) // check ridx bounds
01750               hfine[ridxf+tidx*numRf]++;  
01751           }
01752         }
01753       }
01754     }
01755 
01756     // Pull the best line out of the fine table
01757     int max_val_f = -1, max_i_f = 0;
01758     for (int i = 0; i < numRf*numThetaf; i++) {
01759       if (max_val_f < hfine[i]) {
01760         max_val_f = hfine[i];
01761         max_i_f = i;
01762       }
01763     }
01764     float bestRf = (max_i_f%numRf)*frSpread/numRf + fminR;
01765     float bestThetaf = (max_i_f/numRf)*fthetaSpread/numThetaf + fminTheta;
01766 
01767     // Add viable segments to the vector
01768     // Step along the line, looking for segments. 
01769     float lineLen = 500;
01770     // Can be simpler
01771     const float x1 = bestRf*cos(bestThetaf), y1 = bestRf*sin(bestThetaf);
01772     const float x2 = x1+lineLen*std::cos(bestThetaf+(float)M_PI_2), y2 = y1+lineLen*std::sin(bestThetaf+(float)M_PI_2);
01773     const float x3 = x1-lineLen*std::cos(bestThetaf+(float)M_PI_2), y3 = y1-lineLen*std::sin(bestThetaf+(float)M_PI_2);
01774     
01775     Point e1(x3,y3), e2(x2,y2);
01776     Shape<LineData> line(ShS,e1,e2);
01777     line->setColor(skinny->getColor());
01778     line->setParentId(skinny->getViewableId());
01779 
01780     // m is slope of this line (not the normal line)
01781     float m = (y1 != 0) ? -(x1/y1) : BIG_SLOPE;
01782     float b = y1 - m*x1;
01783     
01784     NEW_SKETCH_N(skelDist,uint,visops::mdist(skinny));
01785     if ( abs(m) <= 1 )  // when |slope| <= 1, scan along x, else scan along y
01786       line->scanHorizForEndPts(skelDist,occlusions,m,b,0);
01787     else
01788       line->scanVertForEndPts(skelDist,occlusions,m,b,0);
01789     
01790     //! Check to see if any endpoints are near any edge of the screen.
01791     //! If they are, invalidate them, assuming that line continues beyond screen.
01792     line->end1_pt.checkValidity(width,height,beg_dist_thresh);
01793     line->end2_pt.checkValidity(width,height,beg_dist_thresh);
01794 
01795     // Experimental: check how much of the line actually lies along edge pixels
01796     // only include lines that are at least 50% above edge pixels. 
01797     
01798     float len = line->getLength();
01799     EndPoint start = line->end1_pt;
01800     EndPoint dx = (line->end2_pt - start)/len;
01801     float onCount = 0;
01802     for (float i=0; i<len; i++) {
01803       EndPoint cur = start+dx*i;
01804       if (skinny((int)(cur.coordX()), (int)(cur.coordY()))) {
01805         onCount++;
01806       }
01807     }
01808 
01809     std::cout << "Line with " << onCount << " / " << len << " pixels";
01810     if (line->isLongerThan(minLength) && (onCount / len) > .5 ) {
01811         std::cout<<" accepted";
01812         result.push_back(line);
01813       }
01814     std::cout<<std::endl;
01815 
01816     // Empty the corresponding bins in the rough hough table
01817     // std::cout<<max_i%numR<<","<<max_i/numR<<std::endl;
01818     for(int row = -5; row <= 5; row++) {
01819       for (int col = -(5-abs(row)); col <=5-abs(row); col++) {
01820         int index = (max_i - (max_i%numR) + ((max_i+col)%numR)+row*numR + numR*numTheta)%(numR*numTheta);
01821         std::cout << "Empty hough " << (index%numR) << "," << (index/numR) << " " << std::endl;
01822         htab[(max_i - (max_i%numR) + ((max_i+col)%numR)+row*numR + numR*numTheta)%(numR*numTheta)] = 0;
01823       }
01824       //std::cout<<std::endl;
01825     }
01826     //std::cout<<std::endl;
01827   }
01828   
01829   return result;
01830 }
01831 
01832 bool LineData::linesParallel(Shape<LineData> l1, Shape<LineData>l2)
01833 {
01834   const float maxDTheta = .15f, minDThetaR = 10.0f;
01835   float dTheta = l1->getOrientation() - l2->getOrientation(); 
01836   if (dTheta > (float)M_PI_2)
01837     dTheta = dTheta - (float)M_PI;
01838   if (std::abs(dTheta) < maxDTheta)
01839     {
01840       if (std::abs (100*dTheta + (l1->rho_norm - l2->rho_norm)) > minDThetaR)
01841   return true;
01842     }
01843   return false;
01844 }
01845 
01846 LineData& LineData::operator=(const LineData& other) {
01847   if (&other == this)
01848     return *this;
01849   BaseData::operator=(other);
01850   end1_pt = other.end1_pt;
01851   end2_pt = other.end2_pt;
01852   rho_norm = other.rho_norm;
01853   theta_norm = other.theta_norm;
01854   orientation = other.orientation;
01855   length = other.length;
01856   return *this;
01857 }
01858 
01859 
01860 // Compute if two points are on the same side of the line
01861 bool LineData::pointsOnSameSide(const Point& p1, const Point& p2)
01862 {
01863   float dx = end2_pt.coordX() - end1_pt.coordX();
01864   float dy = end2_pt.coordY() - end1_pt.coordY();
01865   
01866   float p1val = (p1.coordY() - end1_pt.coordY())*dx - (p1.coordX() - end1_pt.coordX())*dy;
01867   float p2val = (p2.coordY() - end1_pt.coordY())*dx - (p2.coordX() - end1_pt.coordX())*dy;
01868 
01869   return (p1val>0) == (p2val>0);
01870 }
01871 
01872 // Checks if the distance from the line to the point is less than 1 pixel,
01873 // and checks that the point is within the end points of the line.
01874 bool LineData::pointOnLine(const Point& p)
01875 {
01876   const float BOUNDS_EXTEND = 1;
01877   float dx = end2_pt.coordX() - end1_pt.coordX();
01878   float dy = end2_pt.coordY() - end1_pt.coordY();
01879   float val = (p.coordY() - end1_pt.coordY())*dx - (p.coordX() - end1_pt.coordX())*dy;
01880   val /= length;
01881   bool inBounds = (p.coordX() >= (leftPt().coordX() - BOUNDS_EXTEND)) &&
01882     (p.coordX() <= (rightPt().coordX() + BOUNDS_EXTEND)) &&
01883     (p.coordY() >= (topPt().coordY() - BOUNDS_EXTEND)) &&
01884     (p.coordY() <= (bottomPt().coordY() + BOUNDS_EXTEND));
01885   return (val > -1) && (val < 1) && inBounds;
01886 }
01887 
01888 
01889 // Comparison predicates
01890 
01891 bool LineData::LineDataLengthLessThan::operator() (const LineData &line1, const LineData &line2) const {
01892       return line1.getLength() < line2.getLength();
01893 }
01894 
01895 bool LineData::LengthLessThan::operator() (const Shape<LineData> &line1, const Shape<LineData> &line2) const {
01896       return line1->getLength() < line2->getLength();
01897 }
01898 
01899 bool LineData::ParallelTest::operator() (const Shape<LineData> &line1, const Shape<LineData> &line2) const {
01900   return angdist(line1->getOrientation(),line2->getOrientation()) <= tolerance;
01901 }
01902 
01903 bool LineData::PerpendicularTest::operator() (const Shape<LineData> &line1, const Shape<LineData> &line2) const {
01904   return angdist(angdist(line1->getOrientation(),line2->getOrientation()), (orientation_t)M_PI/2) <= tolerance;
01905 }
01906 
01907 bool LineData::ColinearTest::operator() (const Shape<LineData> &line1, const Shape<LineData> &line2) const {
01908   return ParallelTest(ang_tol)(line1,line2) && 
01909     abs( line1->getRhoNorm() - line2->getRhoNorm() ) <= dist_tol;
01910 }
01911 
01912 bool LineData::IsHorizontal::operator() (const Shape<LineData> &line) const {
01913       const AngPi orient = line->getOrientation();
01914       return  (orient <= threshold) || (orient >= M_PI - threshold);
01915     }
01916 
01917 bool LineData::IsVertical::operator() (const Shape<LineData> &line) const {
01918       const AngPi orient = line->getOrientation();
01919       return  (orient >= threshold) && (orient <= M_PI - threshold);
01920     }
01921 
01922 } // namespace

DualCoding 5.1CVS
Generated Mon May 9 04:56:26 2016 by Doxygen 1.6.3