I guys, I am trying to implement a ray tracing algorithm to check if a point lies inside a polygon. I got this python code and I am trying to implement the following part in CP-SAT
int pnpoly(int nvert, double *vertx, double *verty, double testx, double testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
here you can see that there are a lot of conditions to check inside if().
What I have tried
I tried to implement the same code for 2 points and 1 test point using CP_SAT and it looks something like this:
from ortools.sat.python import cp_model
model = cp_model.CpModel()
# 2 points
vert_x1 = model.NewIntVar(1, 1, 'x1')
vert_y1 = model.NewIntVar(4, 4, 'y1')
vert_x2 = model.NewIntVar(6, 6, 'x2')
vert_y2 = model.NewIntVar(0, 0, 'y2')
# 1 test point
testx = model.NewIntVar(2, 2, 'xc')
testy = model.NewIntVar(2, 2, 'yc')
result = model.NewBoolVar("result")
# booleans for y-axis conditions
b1 = model.NewBoolVar('b1')
b2 = model.NewBoolVar('b2')
b3 = model.NewBoolVar('b3')
#booleans for x-axis conditions
c1 = model.NewBoolVar('c1')
# booleans for final decision
d1 = model.NewBoolVar('d1')
d2 = model.NewBoolVar('d2')
#temp_variable for x_conditions
diff_1 = model.NewIntVar(-100, 100, 'diff_1')
diff_2 = model.NewIntVar(-100, 100, 'diff_2')
diff_3 = model.NewIntVarFromDomain(cp_model.Domain.FromIntervals([[-100, -1]]), 'diff_3')
add_1 = model.NewIntVar(-100, 100, 'add_1')
div_1 = model.NewIntVar(-100, 100, 'div_1')
prod_1 = model.NewIntVar(-1000, 1000, 'prod_1')
"""
(verty[i] > testy) && (verty[j] > testy)
"""
#check conditions for y_axis
model.Add(vert_y1 > testy).OnlyEnforceIf(b1)
model.Add(vert_y1 <= testy).OnlyEnforceIf(b1.Not())
model.Add(vert_y2 > testy).OnlyEnforceIf(b2)
model.Add(vert_y2 <= testy).OnlyEnforceIf(b2.Not())
model.Add(b1 != b2).OnlyEnforceIf(b3)
model.Add(b1 == b2).OnlyEnforceIf(b3.Not())
"""
(testx < (vertx[j] - vertx[i]) * (testy - vert[i]) / (verty[j] - verty[i]) + vertx[i] )
"""
#check conditions for x_axis
model.Add(diff_1 == vert_x2 - vert_x1)
model.Add(diff_2 == testy - vert_y1)
model.Add(diff_3 == vert_y2 - vert_y1)
model.AddMultiplicationEquality(prod_1, diff_1, diff_2)
model.AddDivisionEquality(div_1, prod_1, diff_3 )
model.Add(add_1 == div_1 + vert_x1)
model.Add(testx < add_1).OnlyEnforceIf(c1)
model.Add(testx >= add_1).OnlyEnforceIf(c1.Not())
"""
check if both conditions are True)
"""
model.Add(c1 == 1).OnlyEnforceIf(d1)
model.Add(b3 == 1).OnlyEnforceIf(d2)
"""toggle result
(I dont know how to implement toggle)
"""
model.Add(d1 != d2).OnlyEnforceIf(result.Not())
solver = cp_model.CpSolver()
solver.parameters.log_search_progress = True
status = solver.Solve(model)
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print("b1 = ", solver.Value(b1))
print("b2 = ", solver.Value(b2))
print("b3 = ", solver.Value(b3))
print("c1 = ", solver.Value(c1))
print("d1 = ", solver.Value(d1))
print("d2 = ", solver.Value(d2))
print("result : ", solver.Value(result))
print("")
print("diff_1", solver.Value(diff_1))
print("diff_2", solver.Value(diff_2))
print("diff_3", solver.Value(diff_3))
print("prod_1 = ", solver.Value(prod_1))
print("div_1", solver.Value(div_1))
print("addd_1", solver.Value(add_1))
this gives results for 2 points but I want to use for multiple 2d points.
**Issue: **
You can see that there is a division in the calculation having division as diff_3 and when I write the domain of diff_3 as:
diff_3 = model.NewIntVarFromDomain(cp_model.Domain.FromIntervals([[-100, -1], [1, 100]]), 'diff_3')
then it says the The Divisior cannot span accross zero in constraint
But how is it possible as the domain does not have a "0" in it.
**Also, **
How one can toggle a boolean variable in CP - SAT?
Any help will be appreciated :)
EDIT:
Thanks to to Laurent's answer. I hereby have implemented his code for multiple points. The instance of the code is written below:
denum_pos = [model.NewIntVar(1, 100, 'denum_pos_{i}') for i in range(num_points)]
denum_neg = [model.NewIntVar(-100, -1, 'denum_neg_{i}') for i in range(num_points)]
denum_is_positive = [model.NewBoolVar('denum_is_positive_{i}') for i in range(num_points)]
target_pos = [model.NewIntVar(1, 100, 'target_pos_{i}') for i in range(num_points)]
target_neg = [model.NewIntVar(-100, -1, 'target_neg_{i}') for i in range(num_points)]
for i in range(0, num_points):
j = i - 1 if i != 0 else num_points - 1
model.Add(diff_1_arr[i] == xp[j] - xp[i])
model.Add(diff_2_arr[i] == yc - yp[i])
model.Add(diff_3_arr[i] == yp[j] - yp[i])
model.AddMultiplicationEquality(prod_1_arr[i], diff_1_arr[i], diff_2_arr[i])
#perform division
model.Add(diff_3_arr[i] > 0).OnlyEnforceIf(denum_is_positive[i])
model.Add(diff_3_arr[i] < 0).OnlyEnforceIf(denum_is_positive[i].Not())
model.Add(denum_pos[i] == diff_3_arr[i]).OnlyEnforceIf(denum_is_positive[i])
model.Add(denum_pos[i] == 1).OnlyEnforceIf(denum_is_positive[i].Not())
model.Add(denum_neg[i] == 1).OnlyEnforceIf(denum_is_positive[i])
model.Add(denum_neg[i] == diff_3_arr[i]).OnlyEnforceIf(denum_is_positive[i].Not())
model.AddDivisionEquality(target_pos[i], prod_1_arr[i], denum_pos[i])
model.AddDivisionEquality(target_neg[i], prod_1_arr[i], denum_neg[i])
model.Add(div_1_arr[i] == target_pos[i]).OnlyEnforceIf(denum_is_positive[i])
model.Add(div_1_arr[i] == target_neg[i]).OnlyEnforceIf(denum_is_positive[i].Not())
model.Add(add_1_arr[i] == div_1_arr[i] + xp[i])
And the domain of denum is as such:
diff_3_arr = [model.NewIntVar(-1000, 1000, 'diff_3_arr_{i}') for i in range(num_points)]
When I comment out this whole part I get some output hence there should be no problem in the other constriants.
For the division with domains spanning across 0, you need to do a little coding
model.AddDivisionEquality(target, numerator, denum)
can be encoded as follows:
create 2 variables denum_neg and denum_pos with each half domain (positive and negative).
create a Boolean variable denum_is_positive.
link denum_is_positive, denum, and denum_neg, denum_pos
model.Add(denum > 0).OnlyEnforceIf(denum_is_positive)
model.Add(denum < 0).OnlyEnforceIf(denum_is_positive.Not())
model.Add(denum_pos == denum).OnlyEnforceIf(denum_is_positive)
model.Add(denum_pos == 1).OnlyEnforceIf(denum_is_positive.Not())
model.Add(denum_neg == -1).OnlyEnforceIf(denum_is_positive)
model.Add(denum_neg == denum).OnlyEnforceIf(denum_is_positive.Not())
now duplicate the division
model.AddDivisionEquality(target_pos, numerator, denum_pos)
model.AddDivisionEquality(target_neg, numerator, denum_neg)
channel back with `denum_is_positive)
model.Add(target == target_pos).OnlyEnforceIf(denum_is_positive)
model.Add(target == target_neg).OnlyEnforceIf(denum_is_positive.Not())