# CONTAIN - Editorial

Author: Vikas Pandey
Tester: Felipe Mota
Editorialist: Rajarshi Basu

Medium

# PROBLEM:

You are given N \leq 5000 points (X_1, Y_1), (X_2, Y_2), \ldots, (X_N, Y_N) in a 2D plane. We also given Q \leq 2000 queries, where in each query we are given a point p = (x, y), where we place a candle.

We have to divide them into a non-negative number of layers. The deliciousness of the cake to be equal to the number of layers in the cake. A layer is a simple polygon whose vertices are some of N points. It is not required for each point to belong to a layer. The layers must satisfy the following conditions:

• For each layer, the candle must lie strictly inside it.
• Each of Chef’s N points that lies strictly inside some layer must belong to some other layer or be the candle point.
• No two layers may touch or intersect.
• The layers must be chosen in such a way that the deliciousness of the cake is maximum possible.

For every candle point given in the queries, find the maximum deliciousness of the cake.

# QUICK EXPLANATION:

The tiny little explanation
• Keep forming Convex hulls iteratively, one inside the other
• For every query, binary search the first convex hull inside which the point lies.

# EXPLANATION:

Simplifcation

Layer is effectively a convex hull.

Observation 1

Every layer is effectively completely contained inside another outer layer.

Observation 2

Nesting of layers is better than having two non-nested layers. Eg, Here, Fig1 is never better than Fig2, in that sense that no point in Fig1 is “covered/enclosed” by more convex hulls than the same corresponding point in Fig2.

Hint 1

Instead of trying to build the solution from inside out, try to build it from outside in.

Full Solution

We are essentially using Observation 2 as the basis of our solution. Here is what we will do:

1. Build a convex hull from the set of points in your Active Set. Initially, all points are in the Active Set.
2. Remove the points taken in the hull from your Active Set.
3. Repeat Step 1, as long as Active Set is not empty.

This overall algorithm takes O(N^2logn) time, and is our preprocesing step. After this,

After this process ends, we will have a series of nested Convex hulls. Now, how to deal with the queries?

Observation 3

if it is in a certain hull, it will also be inside all of the outer hulls.

Query Details

For every query, just check which all convex hulls contain the give point, or alternatively, binary search on the first hull that contains it. Checking whether a point lies inside a convex polygon can be done in O(N) time, or O(logN) time where the polygon has N points. (check links below for more details). Further, depending on searching all convex hulls, total complexity of answering a query is either O(N) (it might seem at first that it is O(N^2) since we can have O(N) convex hulls, but since the overall number of points is bounded by N, hence the total cost amortises), or O(log^2n)

But... isn't Rule 2 violated?

Well… no. It’s more of the Simplification being tweaked a little. Let all the points that belonged to convex hulls which do not include the candle point, be called Floating Points. Turns out, the first/smallest convex hull that covers the candle point can be made non-convex to pass over all the Floating Points, without excluding the Candle Point.

Time Complexity

O(N^2logn) preprocessing, O(N) or O(log^2n) per query. Hence overall O(QN + N^2logn) or O(Qlog^2n + N^2logn) depending on implementation.

# MORE LINKS FOR CONVEX HULL STUFF:

What you will find in these links are how to implement convex hulls, and how to check whether a point lies inside a polygon, both simple or convex, their time complexities, and their codes ( and much more) . Have Fun!

# QUICK REMINDERS:

Remind Me

Don’t forget to take care of colinear points!

# SOLUTIONS:

Setter's Code
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define pb push_back
#define ff first
#define ss second
// #define endl "\n"

struct point
{
int x, y;
point(int xx, int yy):x(xx), y(yy){}
point operator+(const point &p) const {return point(x+p.x, y+p.y);}
point operator-(const point &p) const {return point(x-p.x, y-p.y);}
bool operator<(const point &p) const {return x<p.x || (x==p.x and y<p.y);}
bool operator==(const point &p) const {return x==p.x and y==p.y;}
int cross(const point &p) const {return x*p.y-y*p.x;}
int dot(const point &p) const {return x*p.x+y*p.y;}
int cross(const point &a, const point &b) const {return (a-*this).cross(b-*this);}
int dot(const point &a, const point &b) const {return (a-(*this)).dot(b-(*this));}
int sqrlen() const {return this->dot(*this);}
};

bool cw(point a, point b, point c){return a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y)<0;}
bool ccw(point a, point b, point c){return a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y)>0;}

struct HULL
{
vector<point> vertices;
HULL(vector<point> v):vertices(v)
{
vertices.pb(vertices[0]);
reverse(vertices.begin(), vertices.end());
vertices.pop_back();
int pos=0, n=(int)vertices.size();
for(int i=1; i<n; i++)
if(vertices[i].y<vertices[pos].y or (vertices[i].y==vertices[pos].y and vertices[i].x<vertices[pos].x))
pos=i;
rotate(vertices.begin(), vertices.begin()+pos, vertices.end());
v.clear();
v.pb(vertices[0]);
v.pb(vertices[1]);
for(int i=2; i<=n; i++)
{
if((v[(int)v.size()-1]-v[(int)v.size()-2]).cross(vertices[i%n]-v[(int)v.size()-2])==0)
v.pop_back();
v.pb(vertices[i]);
}
v.pop_back();
vertices=v;
}
bool PIP(point p)
{
int n=(int)vertices.size();
bool ans=true;
for(int i=0; i<n; i++)
ans&=ccw(vertices[i], vertices[(i+1)%n], p);
return ans;
}
};

vector<HULL> hulls;

void convex_hull(vector<point> &a)
{
set<point> allnow;
for(auto z:a)
allnow.insert(z);
while(!((int)allnow.size()<3))
{
point p1=*allnow.begin(), p2=*allnow.rbegin();
vector<point> up, down;
up.pb(p1);
down.pb(p1);
for(auto it=++allnow.begin(); it!=allnow.end(); it++)
{
if(it==prev(allnow.end()) || !ccw(p1, *it, p2))
{
while((int)up.size()>1 and ccw(up[(int)up.size()-2], up[(int)up.size()-1], *it))
up.pop_back();
up.pb(*it);
}
if(it==prev(allnow.end()) || !cw(p1, *it, p2))
{
while((int)down.size()>1 and cw(down[(int)down.size()-2], down[(int)down.size()-1], *it))
down.pop_back();
down.pb(*it);
}
}
for(int i=(int)down.size()-2; i>0; i--)
up.pb(down[i]);
if((up[1]-up[0]).cross(up[(int)up.size()-1]-up[0])==0)
break;
for(auto z:up)
allnow.erase(z);
hulls.pb(HULL(up));
}
}

void solve()
{
int n, q;
cin>>n>>q;
vector<point> arr;
while(n--)
{
int xx, yy;
cin>>xx>>yy;
arr.pb(point(xx, yy));
}
hulls.clear();
convex_hull(arr);
while(q--)
{
int u, v;
cin>>u>>v;
point p(u, v);
int st=0, en=(int)hulls.size()-1;
while(st<=en)
{
int mid=(st+en)>>1;
if(!hulls[mid].PIP(p))
en=mid-1;
else
st=mid+1;
}
cout<<en+1<<endl;
}
}

int32_t main()
{
ios_base::sync_with_stdio(false);
cin.tie(NULL);
cout.tie(NULL);
int t=1;
cin>>t;
while(t--)
solve();
return 0;
}

Tester's Code
#include <bits/stdc++.h>
using namespace std;
const double E = 1e-9, pi = 2 * acos(0);
struct Point {
int x, y, id;
};
Point operator - (Point a, Point b){
return (Point){a.x - b.x, a.y - b.y, -1};
}
long long cross(Point a, Point b){
return 1ll * a.x * b.y - 1ll * a.y * b.x;
}
long long cross(Point a, Point b, Point c){
return cross(b - a, c - a);
}
vector<Point> convex_hull(vector<Point> points){
sort(points.begin(), points.end(), [&](Point a, Point b){
if(a.y != b.y) return a.y < b.y;
return a.x < b.x;
});
vector<Point> upper, lower;
for(int i = 0; i < (int)points.size(); i++){
while(upper.size() >= 2 && cross(upper.end()[-2], upper.end()[-1], points[i]) > 0)
upper.pop_back();
upper.push_back(points[i]);
}
for(int i = (int)points.size() - 1; i >= 0; i--){
while(lower.size() >= 2 && cross(lower.end()[-2], lower.end()[-1], points[i]) > 0)
lower.pop_back();
lower.push_back(points[i]);
}
if(lower.size() > 2) upper.insert(upper.end(), lower.begin() + 1, lower.end() - 1);
return upper;
}
// 1 => Strictly inside; -1 => Border; 0 => Outside
int point_in_poly(const vector<Point> & poly, Point p){
int many = 0;
for(int i = 0; i < (int)poly.size(); i++){
Point a = poly[i], b = poly[i + 1 < (int) poly.size() ? i + 1 : 0];
if(a.x > b.x) swap(a, b);
if(a.x <= p.x && p.x <= b.x){
if(abs(a.x - b.x) == 0){
if(min(a.y, b.y) <= p.y && p.y <= max(a.y, b.y)) return -1;
} else {
double y = a.y + 1. * (b.y - a.y) / (b.x - a.x) * (p.x - a.x);
if(abs(y - p.y) <= E) return -1;
if(y >= p.y && p.x < b.x) many++;
}
}
}
return many % 2;
}
int main(){
int t;
cin >> t;
while(t--){
int n, q;
cin >> n >> q;
vector<Point> poly(n);
vector<int> not_used(n, 1);
for(int i = 0; i < n; i++){
cin >> poly[i].x >> poly[i].y;
poly[i].id = i;
}
vector<vector<Point>> rings;
while(accumulate(not_used.begin(), not_used.end(), 0) >= 3){
vector<Point> can;
for(int i = 0; i < n; i++){
if(not_used[i]){
can.push_back(poly[i]);
}
}
vector<Point> hull = convex_hull(can);
for(auto & p : hull)
not_used[p.id] = 0;
rings.push_back(hull);
}
for(int i = 0; i < q; i++){
Point p;
cin >> p.x >> p.y;
int ans = 0;
for(auto & ring : rings){
if(point_in_poly(ring, p) == 1)
ans += 1;
}
cout << ans << '\n';
}
}
return 0;
}


Please give me suggestions if anything is unclear so that I can improve. Thanks

14 Likes

please tell me why my submission giving WA for 15 point
https://www.codechef.com/viewsolution/34415553

It would be great if someone make a video solution.

4 Likes

It gives a wrong answer on this test case:

1
4 1
0 0
5 0
3 0
0 5
2 0

Can someone please say why my submission gives WA?
My approach:
Consider all n points. Take their convex hull and store them in a shapes vector. Now delete all those point from the n points. Be doing this until you have less than 3 points. Ignore the leftover points.
Now I reversed my shapes vector. I ran a loop from i = 0 to i = shapes.size()-1. If pnpoly(shapes[i], p) == true (pnpoly(vert, test) returns true if point test is strictly in a polygon of vertices vert), answer is shapes.size()-i.

Your code gives a WA on this case:

1
8 1
0 0
1 0
2 0
2 1
2 2
1 2
0 2
0 1
1 1

Expected output is 1.

1 Like

Your Code is not working when the query point lies on the boundary of the Polygon.

2 Likes

Anybody else that had a problem in these specific test cases please tell me how to fix them. Here is the link to my code.

Thanks. This was my favorite problem in this contest. Very nice question

3 Likes

live discussion

Can anyone please provide me a test case where my code fails for subtask 2,45 points https://www.codechef.com/viewsolution/34312825

I came across that issue while solving the problem . My initial code was:

int pnpoly(vector <pt> vert, pt test)
{
int i, j, c = 0;
for (i = 0, j = vert.size()-1; i < vert.size(); j = i++) {
if ( ((vert[i].y>test.y) != (vert[j].y>test.y)) &&
(test.x < (vert[j].x-vert[i].x) * (test.y-vert[i].y) / (vert[j].y-vert[i].y) + vert[i].x))
c = !c;
}
return c;
}


This failed when the query point was one of the points in the convex hull. So I fixed that by adding a loop to check if the query point was one of the points in the convex hull.

int pnpoly(vector <pt> vert, pt test)
{
for (int i = 0; i < vert.size(); ++i) if (vert[i].x == test.x && vert[i].y == test.y) return 0;
int i, j, c = 0;
for (i = 0, j = vert.size()-1; i < vert.size(); j = i++) {
if ( ((vert[i].y>test.y) != (vert[j].y>test.y)) &&
(test.x < (vert[j].x-vert[i].x) * (test.y-vert[i].y) / (vert[j].y-vert[i].y) + vert[i].x))
c = !c;
}
return c;
}


This fixed that issue.
But can you please tell me how to detect if the query point is on the boundary?

https://www.codechef.com/viewsolution/34036993
Above solution i have coded in java giving me TLE for only 2 cases
https://www.codechef.com/viewsolution/34070969
Now this exactly same solution in cpp giving me correct answer
Please explain me why???

The main idea is to compute the dot-product and cross-product of the Query point and any two vertices and Check

1. If the cross-product is 0 this tells you if all the three points are aligned.

2. If the dot-product is positive and is less than the square of the distance between a and b.

You have to check the above conditions for all the vertices.

This is my code for checking if the point lies on the boundary of the polygon or not.

def pointOnPolygon(polygonVertices, point):
def isBetween(a, b, c):
crossproduct = (c[1] - a[1]) * (b[0] - a[0]) - (c[0] - a[0]) * (b[1] - a[1])
if abs(crossproduct) > 2.220446049250313e-16:
return False
dotproduct = (c[0] - a[0]) * (b[0] - a[0]) + (c[1] - a[1])*(b[1] - a[1])
if dotproduct < 0:
return False
squaredlength = pow(b[0] - a[0], 2) + pow(b[1] - a[1], 2)
if dotproduct > squaredlength:
return False
return True
n = len(polygonVertices)
for i in range(n):
p1 = polygonVertices[i]
p2 = polygonVertices[-n+i+1]
if isBetween(p1, p2, point):
return True
return False


If somebody would have coded in R language would have got in 10-15 lines because of this function that calulates convex hull peeling depth !

5 Likes

sir I fixed that issue
still geting WA
https://www.codechef.com/viewsolution/34448980

Do seeing onion peeling generally inspire you to write a solution?
https://en.wikipedia.org/wiki/Convex_layers#:~:text=In%20this%20context%2C%20the%20number,this%20notion%20of%20data%20depth.&text=convex%20layers%2C%20as%20do%20the,points%20within%20any%20convex%20shape.
cheers!

June long challenge editorial beginner friendly video explanation with animation and code

Delicious cake (CONTAIN) : https://youtu.be/tXD2yEVA0pM
Tom and jerry (EOEO) : https://youtu.be/ZWI5n0Ogir0
Even matrix (EVENM) : https://youtu.be/KA8WoO7jCg8

2 Likes

Can someone please tell me why my submission is giving WA?

How to optimally check for colinear points on the polygon?
I noticed that my convex hull function was only taking the two farthest points if there are colinear points on the polygon. In order to fix this i had to sacrifice the efficiency of my code and could not get more than 15 points.