convex hull example script - ObjectVision/GeoDMS GitHub Wiki
The following script can be used for a simple convex hull with only one ring. With iterations the number of points can be simplified.
container Simplified := Simplify(iter_domain, 50); // domain and number of iterations
template Simplify
{
// begin case parameters
unit<uint32> iter_domain;
parameter<uint32> nrIterations;
// end case parameters
unit<uint32> iter := range(uint32, 0, nrIterations)
{
attribute<string> name := 'I' + string(id(.));
attribute<string> PrevName := MakeDefined(name[sub_or_null(id(.), 1)], 'StartingState');
}
container StartingState
{
unit<uint32> NextValue := sequence2points(iter_domain/geometry_rd);
}
container Iters := for_each_ne(iter/name, 'IterT('+iter/PrevName+', iter_domain)');
container LastIter := =last('Iters/'+iter/name);
template IterT
{
// begin case parameters
container PrevIterator;
unit<uint32> domain;
// end case parameters
unit<uint32> PrevValue := PrevIterator/NextValue;
container ConvexHullT := MakeConvexHullT(domain, PrevValue);
unit<uint32> NextValue := ConvexHullT/sequence/Convex_hull_points;
}
container MakeFinal
{
unit <uint32> domain_with_endpoints := union_unit(LastIter/NextValue, iter_domain)
{
attribute<geometries/rdc> point :=
union_data(., LastIter/NextValue/point, first(LastIter/NextValue/point, LastIter/NextValue/sequence_rel));
attribute<uint32> sequence_rel := union_data(., LastIter/NextValue/sequence_rel, id(iter_domain));
}
attribute<geometries/rdc> convex_hull (iter_domain, poly) :=
points2polygon(domain_with_endpoints/point, domain_with_endpoints/sequence_rel);
}
template MakeConvexHullT
{
// begin case parameters
unit<uint32> domain;
unit<uint32> seq;
// end case parameters
unit<uint32> sequence := seq
{
attribute<geometries/rdc> point := seq/point;
attribute<uint32> sequence_rel := seq/sequence_rel;
attribute<uint32> min_index (domain) := min_index(id(.), sequence_rel);
attribute<uint32> max_index (domain) := max_index(id(.), sequence_rel);
attribute<bool> IsFirst := id(.) == min_index[sequence_rel];
attribute<bool> IsLast := id(.) == max_index[sequence_rel];
attribute<uint32> prev_id := IsFirst ? rjoin(sequence_rel, id(domain), max_index) : sub_or_null(id(.), 1);
attribute<uint32> next_id := IsLast ? rjoin(sequence_rel, id(domain), min_index) : id(.)+1;
attribute<geometries/rdc> A := point[prev_id];
attribute<geometries/rdc> B := point;
attribute<geometries/rdc> C := point[next_id];
attribute<geometries/rdc> p := B - A;
attribute<geometries/rdc> q := C - B;
attribute<float32> Px := pointcol(p);
attribute<float32> Py := pointrow(p);
attribute<float32> Qx := pointcol(q);
attribute<float32> Qy := pointrow(q);
attribute<float32> det := Px*Qy - Py*Qx;
unit<uint32> Convex_hull_points :=
select_with_org_rel(det <= 0f && (point != rjoin(prev_id, id(.), point)))
{
attribute<geometries/rdc> point := ../point[org_rel];
attribute<uint32> sequence_rel := ../sequence_rel[org_rel];
}
}
}
}
This script computes the exact convex hull of a point set. pareto_optimal solves a related but different problem: it keeps the elements of a partition that are not dominated by any other element on a set of minimised criteria. The two are easy to confuse because both can be computed with one sort and one sweep, but the sweep tests something different at each step: pareto_optimal compares a point against the best seen so far on every criterion, while this script's det <= 0 test is an orientation test over three consecutive points of the hull built so far.
For two criteria the two results relate in one direction only. Every point on the side of the hull facing the origin, the side a model wants when both criteria are minimised, is also Pareto-optimal, but a Pareto-optimal point does not have to be on the hull: it can be undominated by any single other point and still lie inside the hull of the other Pareto-optimal points, on a concave dent of the front. Multi-objective optimisation theory calls this the difference between supported and unsupported Pareto-optimal points. So pareto_optimal cannot stand in for this script, nor the other way round, and a hull computed on an already Pareto-optimal point set still needs this script's orientation test, not pareto_optimal's dominance test.
GeoDMS#1331 proposes an actual operator for this hull, and a second mode restricted to the side facing the origin, for the two-criterion case where that distinction matters.