a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
This means that the merging routine doesn't need to be aware of direction,
but the distance routine does.
-*/
-/* If we start including very thin buildings, numerical accuracy errors can
- arise. Therefore, we ignore all buildings that are less than epsilon wide. */
-#define EPS 1e-5
+ From 2007 through 2012, buildings of width less than EPS were discarded,
+ citing numerical accuracy concerns. We remember that floating point
+ comparisons of nearly-equal values can be affected by rounding error.
+ Also, some target machines use the x87 floating point unit, which provides
+ extended precision for intermediate results held in registers. On this type
+ of hardware comparisons such as
+ double c = 1.0/3.0; boolean compare = (c == 1.0/3.0)
+ could go either way because the 1.0/3.0 is allowed to be kept
+ higher precision than the variable 'c'.
+ Alert to these considerations, we now accept buildings of zero-width.
+*/
static void
print_buildings (list<Building> const &b)
assert (!isinf (slope_) && !isnan (slope_));
if (isinf (start))
- {
- assert (start_height == end_height);
- y_intercept_ = start_height;
- }
- else
- y_intercept_ = start_height - slope_ * start;
+ assert (start_height == end_height);
+
+ start_height_ = start_height;
}
+// Because of the way slope is calculated as a ratio of usually small
+// differences, its precision is not sufficient for extrapolating
+// significantly from the original interval. For that reason, we
+// don't store the y0 value that would allow more precalculation by using
+// y = x * slope + y0
+// rather than
+// y = (x - xstart) * slope + ystart
+
inline Real
Building::height (Real x) const
{
- return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
+ return (isinf (x) || isinf (start_)) ? start_height_
+ : slope_ * (x - start_) + start_height_;
}
void
Building::print () const
{
- printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
+ printf ("%f (x - x0) + %f from %f to %f\n", slope_, start_height_, start_, end_);
}
inline Real
Building::intersection_x (Building const &other) const
{
- Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
+ Real ret = start_
+ - (other.start_height_ - start_height_
+ - other.slope_ * (other.start_ - start_))
+ /(other.slope_ - slope_);
+ /* A numerically symmetric version would be
+ x = (x1+x0)/2 - (y1-y0 - (s1+s0)(x1-x0)/2)/(s1-s0)
+ */
return isnan (ret) ? -infinity_f : ret;
}
end_ = chop;
}
-// Returns a shift s such that (x + s, y) intersects the roof of
-// this building. If no such shift exists, returns infinity_f.
-Real
-Building::shift_to_intersect (Real x, Real y) const
-{
- // Solve for s: y = (x + s)*m + b
- Real ret = (y - y_intercept_ - slope_ * x) / slope_;
-
- if (ret >= start_ && ret <= end_ && !isinf (ret))
- return ret;
- return infinity_f;
-}
-
static Real
first_intersection (Building const &b, list<Building> *const s, Real start_x)
{
// conceals and intersection_x involve multiplication and
// division. Avoid that, if we can.
- if (c.y_intercept_ == -infinity_f)
+ if (c.start_height_ == -infinity_f)
{
if (c.end_ > b.end_)
return b.end_;
return b.end_;
}
+// conceals returns true if *this is strictly above other at x
+
bool
Building::conceals (Building const &other, Real x) const
{
- if (slope_ == other.slope_)
- return y_intercept_ > other.y_intercept_;
-
- /* their slopes were not equal, so there is an intersection point */
- Real i = intersection_x (other);
- return (i <= x && slope_ > other.slope_)
- || (i > x && slope_ < other.slope_);
+ return height (x) > other.height (x);
}
// Remove redundant empty buildings from the skyline.
for (i = buildings_.begin (); i != buildings_.end (); i++)
{
- if (last_empty && i->y_intercept_ == -infinity_f)
+ if (last_empty && i->start_height_ == -infinity_f)
{
list<Building>::iterator last = i;
last--;
buildings_.erase (i);
i = last;
}
- last_empty = (i->y_intercept_ == -infinity_f);
+ last_empty = (i->start_height_ == -infinity_f);
}
assert (buildings_.front ().start_ == -infinity_f);
for (right = buildings_.begin (); right != buildings_.end (); right++)
{
- if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
+ if (center != buildings_.begin () && center->start_height_ == -infinity_f)
{
Real p1 = left->height (left->end_);
- Real p2 = right->height (right->start_);
+ Real p2 = right->start_height_;
*center = Building (center->start_, p1, p2, center->end_);
left = center;
// Optimization: if the other skyline is empty at this point,
// we can avoid testing some intersections. Just grab as many
// buildings from s1 as we can, and shove them onto the output.
- if (c.y_intercept_ == -infinity_f
+ if (c.start_height_ == -infinity_f
&& c.end_ >= b.end_)
{
list<Building>::iterator i = s1->begin ();
while (i != s1->end () && i->end_ <= c.end_)
i++;
+ s1->front ().start_height_ = s1->front ().height (x);
s1->front ().start_ = x;
result->splice (result->end (), *s1, s1->begin (), i);
x = result->back ().end_;
Real end = first_intersection (b, s2, x);
if (s2->empty ())
{
+ b.start_height_ = b.height (last_end);
b.start_ = last_end;
result->push_back (b);
break;
}
- /* only include buildings wider than epsilon */
- if (end > x + EPS)
+ if (end >= x)
{
b.leading_part (end);
+ b.start_height_ = b.height (last_end);
b.start_ = last_end;
last_end = b.end_;
result->push_back (b);
static void
single_skyline (Building b, list<Building> *const ret)
{
- if (b.end_ > b.start_ + EPS)
- {
- if (b.start_ != -infinity_f)
- ret->push_back (Building (-infinity_f, -infinity_f,
- -infinity_f, b.start_));
- ret->push_back (b);
- if (b.end_ != infinity_f)
- ret->push_back (Building (b.end_, -infinity_f,
- -infinity_f, infinity_f));
- }
- else
- {
- empty_skyline (ret);
- }
+ assert (b.end_ >= b.start_);
+
+ if (b.start_ != -infinity_f)
+ ret->push_back (Building (-infinity_f, -infinity_f,
+ -infinity_f, b.start_));
+ ret->push_back (b);
+ if (b.end_ != infinity_f)
+ ret->push_back (Building (b.end_, -infinity_f,
+ -infinity_f, infinity_f));
}
/* remove a non-overlapping set of boxes from BOXES and build a skyline
while (i != buildings->end ())
{
Real x1 = i->start_;
- Real y1 = i->height (i->start_);
+ Real y1 = i->start_height_;
Real x2 = i->end_;
Real y2 = i->height (i->end_);
continue;
}
- if (x1 > last_end + EPS)
+ // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
+ if (x1 > last_end)
result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
result.push_back (*i);
bool operator () (Building const &b1, Building const &b2)
{
return b1.start_ < b2.start_
- || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
+ || (b1.start_ == b2.start_ && b1.start_height_ > b2.start_height_);
}
};
/*
Build skyline from a set of boxes.
- Boxes should have fatness in the horizon_axis, otherwise they are ignored.
+ Boxes should be non-empty on both axes. Otherwise, they will be ignored
*/
Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
{
list<Building> buildings;
sky_ = sky;
- Axis vert_axis = other_axis (horizon_axis);
for (vsize i = 0; i < boxes.size (); i++)
- {
- Interval iv = boxes[i][horizon_axis];
- if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
- buildings.push_front (Building (boxes[i], horizon_axis, sky));
- }
+ if (!boxes[i].is_empty (X_AXIS)
+ && !boxes[i].is_empty (Y_AXIS))
+ buildings.push_front (Building (boxes[i], horizon_axis, sky));
buildings_ = internal_build_skyline (&buildings);
normalize ();
/*
build skyline from a set of line segments.
- Buildings should have fatness in the horizon_axis, otherwise they are ignored.
+ Segments can be articulated from left to right or right to left.
+ In the case of the latter, they will be stored internally as left to right.
*/
Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
{
Real y1 = left[other_axis (horizon_axis)] * sky;
Real y2 = right[other_axis (horizon_axis)] * sky;
- if (x1 + EPS < x2)
+ if (x1 <= x2)
buildings.push_back (Building (x1, y1, y2, x2));
}
Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
{
sky_ = sky;
- Building front (b, horizon_axis, sky);
- single_skyline (front, &buildings_);
- normalize ();
+ if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
+ {
+ Building front (b, horizon_axis, sky);
+ single_skyline (front, &buildings_);
+ normalize ();
+ }
}
void
}
/* do the same filtering as in Skyline (vector<Box> const&, etc.) */
- Interval iv = b[a];
- if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
+ if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
return;
my_bld.splice (my_bld.begin (), buildings_);
{
list<Building>::iterator end = buildings_.end ();
for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
- i->y_intercept_ += sky_ * r;
+ i->start_height_ += sky_ * r;
}
void
{
i->start_ += s;
i->end_ += s;
- i->y_intercept_ -= s * i->slope_;
}
}
Skyline
Skyline::padded (Real horizon_padding) const
{
+ if (horizon_padding < 0.0)
+ warning ("Cannot have negative horizon padding. Junking.");
+
+ if (horizon_padding <= 0.0)
+ return *this;
+
list<Building> pad_buildings;
for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
{
if (i->start_ > -infinity_f)
{
- Real height = i->height (i->start_);
+ Real height = i->start_height_;
if (height > -infinity_f)
{
// Add the sloped building that pads the left side of the current building.
list<Building>::const_iterator i;
for (i = buildings_.begin (); i != buildings_.end (); ++i)
{
- ret = max (ret, i->height (i->start_));
+ ret = max (ret, i->start_height_);
ret = max (ret, i->height (i->end_));
}
{
for (list<Building>::const_iterator i (buildings_.begin ());
i != buildings_.end (); i++)
- if (i->y_intercept_ > -infinity_f)
+ if (i->start_height_ > -infinity_f)
return i->start_;
return infinity_f;
{
for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
i != buildings_.rend (); ++i)
- if (i->y_intercept_ > -infinity_f)
+ if (i->start_height_ > -infinity_f)
return i->end_;
return -infinity_f;
Skyline::set_minimum_height (Real h)
{
Skyline s (sky_);
- s.buildings_.front ().y_intercept_ = h * sky_;
+ s.buildings_.front ().start_height_ = h * sky_;
merge (s);
}
if (!buildings_.size ())
return true;
Building b = buildings_.front ();
- return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
+ return b.end_ == infinity_f && b.start_height_ == -infinity_f;
}
void
projects / lilypond.git / blobdiff