11/* *
2- * codac2_SepPolarCart.cpp
2+ * \file codac2_SepPolarCart.cpp
33 * ----------------------------------------------------------------------------
44 * \date 2026
5- * \author Benoît Desrochers, (Simon Rohou)
5+ * \author Simon Rohou, from a former implementation of Benoît Desrochers
66 * \copyright Copyright 2026 Codac Team
77 * \license GNU Lesser General Public License (LGPL)
88 */
99
1010#include " codac2_SepPolarCart.h"
1111#include " codac2_CtcPolar.h"
12+ #include " codac2_cart_prod.h"
13+ #include " codac2_hull.h"
1214
1315using namespace codac2 ;
1416
15- BoxPair SepPolarCart::separate (const IntervalVector& x) const
17+ /* *
18+ * Contracts x with the Cartesian projection of the complement of a polar box.
19+ *
20+ * For a polar box [rho]x[theta], the complement is represented as the union
21+ * of four polar boxes:
22+ * - rho <= rho.lb(), with theta in [theta]
23+ * - rho >= rho.ub(), with theta in [theta]
24+ * - theta before theta.lb(), with rho in R+
25+ * - theta after theta.ub(), with rho in R+
26+ *
27+ * The angular complement is computed in the 2*pi-wide window centered on [theta].
28+ */
29+ IntervalVector contract_polar_box_complement (
30+ const CtcPolar& ctc_polar,
31+ const IntervalVector& x,
32+ const IntervalVector& x_pol)
1633{
1734 assert_release (x.size () == 2 );
35+ assert_release (x_pol.size () == 2 );
1836
19- CtcPolar ctc_polar;
20- IntervalVector x_cart ( 2 ), x_pol (x) ;
37+ if (x. is_empty ())
38+ return x ;
2139
22- ctc_polar.contract (x_cart[0 ], x_cart[1 ], x_pol[0 ], x_pol[1 ]);
23- auto x_cart_sep = _sep.front ()->separate (x_cart);
24- BoxPair x_polar_sep { x, x };
25- ctc_polar.contract (x_cart_sep.inner [0 ], x_cart_sep.inner [1 ], x_polar_sep.inner [0 ], x_polar_sep.inner [1 ]);
26- ctc_polar.contract (x_cart_sep.outer [0 ], x_cart_sep.outer [1 ], x_polar_sep.outer [0 ], x_polar_sep.outer [1 ]);
40+ // Complement of the empty polar box: no point is removed by the inner
41+ // contractor, so the Cartesian box is left unchanged.
42+ if (x_pol.is_empty ())
43+ return x;
2744
28- return x_polar_sep ;
29- }
45+ const Interval& rho = x_pol[ 0 ] ;
46+ const Interval& theta = x_pol[ 1 ];
3047
48+ std::list<IntervalVector> parts;
3149
32- SepPolarXY::SepPolarXY (Interval rho, Interval theta) : rho(rho), theta(theta), Sep(2 ) {
33- rho_m = Interval (0 , rho.lb ());
34- rho_p = Interval (rho.ub (), POS_INFINITY );
35- double limit = theta.mid () - M_PI ;
36- theta_m = Interval (limit, theta.lb ());
37- theta_p = Interval (theta.ub (), limit + 2 *M_PI );
38- cmpl = Interval (0 , 2 *M_PI );
39- }
50+ auto add_polar_part = [&ctc_polar,&parts](const IntervalVector& x,
51+ const Interval& rho,
52+ const Interval& theta)
53+ {
54+ if (rho.is_empty () || theta.is_empty ())
55+ return ;
4056
57+ Interval rho_ (rho), theta_ (theta);
58+ IntervalVector x_ (x);
59+ ctc_polar.contract (x_[0 ],x_[1 ],rho_,theta_);
60+ parts.push_back (x_);
61+ };
4162
42- void SepPolarXY::contractOut (IntervalVector &x_out){
43- Interval th = this -> theta ;
44- Interval r = this -> rho ;
45- this -> ctc . contract (x_out[ 0 ], x_out[ 1 ], r, th );
46- if (x_out[ 0 ]. is_empty () || x_out[ 1 ]. is_empty ())
47- x_out. set_empty ();
48- }
63+ // Radial complement, restricted to the current angular sector. The angular
64+ // outside parts below cover the remaining angles.
65+ if ( std::isfinite (rho. lb ()) && rho. lb () >= 0 )
66+ add_polar_part (x, { 0 ,rho. lb ()}, theta );
67+
68+ if ( std::isfinite (rho. ub ()))
69+ add_polar_part (x, {rho. ub (),oo}, theta);
4970
71+ // Angular complement. If [theta] already covers at least one full turn,
72+ // there is no angular complement in the 2*pi-periodic sense.
73+ if (std::isfinite (theta.lb ()) && std::isfinite (theta.ub ())
74+ && theta.ub () - theta.lb () < 2 *PI )
75+ {
76+ double limit = theta.mid ()-PI ;
77+ add_polar_part (x, {0 ,oo}, {limit,theta.lb ()});
78+ add_polar_part (x, {0 ,oo}, {theta.ub (),limit+2 *PI });
79+ }
5080
51- void SepPolarXY::contractIn (IntervalVector &x_in){
52- Interval x1 (x_in[0 ]); Interval y1 (x_in[1 ]);
53- Interval x2 (x_in[0 ]); Interval y2 (x_in[1 ]);
54- Interval x3 (x_in[0 ]); Interval y3 (x_in[1 ]);
55- Interval x4 (x_in[0 ]); Interval y4 (x_in[1 ]);
56-
57- Interval ALLREALS1 = Interval::POS_REALS ;
58- Interval ALLREALS2 = Interval::POS_REALS ;
59- Interval cmpl1 (cmpl);
60- Interval cmpl2 (cmpl);
61- Interval theta_m_tmp (theta_m);
62- Interval theta_p_tmp (theta_p);
63- Interval rho_m_tmp (rho_m);
64- Interval rho_p_tmp (rho_p);
65-
66- this ->ctc .contract (x1, y1, ALLREALS1 , theta_m_tmp);
67- this ->ctc .contract (x2, y2, ALLREALS2 , theta_p_tmp);
68- this ->ctc .contract (x3, y3, rho_m_tmp, cmpl1);
69- this ->ctc .contract (x4, y4, rho_p_tmp, cmpl2);
70- x_in[0 ] &= (x1 | x2 | x3 | x4);
71- x_in[1 ] &= (y1 | y2 | y3 | y4);
72- if (x_in[0 ].is_empty () || x_in[1 ].is_empty ())
73- x_in.set_empty ();
81+ return x & hull (parts);
7482}
83+
84+ BoxPair SepPolarCart::separate (const IntervalVector& x) const
85+ {
86+ assert_release (x.size () == 2 );
87+
88+ BoxPair x_cart_sep { x, x };
89+
90+ if (x.is_empty ())
91+ return x_cart_sep;
92+
93+ CtcPolar ctc_polar;
94+ IntervalVector x_cart (x), x_pol (2 );
95+
96+ // Cartesian input -> polar enclosure
97+ ctc_polar.contract (x_cart[0 ], x_cart[1 ], x_pol[0 ], x_pol[1 ]);
98+
99+ if (x_cart.is_empty () || x_pol.is_empty ())
100+ return x_cart_sep;
101+
102+ // Separation in the polar space
103+ const BoxPair x_pol_sep = _sep.front ()->separate (x_pol);
104+
105+ // Outer contraction: keep the Cartesian points whose polar coordinates may
106+ // belong to the polar separator output
107+ IntervalVector copy_x_pol_sep = x_pol_sep.outer .subvector (0 ,1 );
108+ ctc_polar.contract (x_cart_sep.outer [0 ], x_cart_sep.outer [1 ], copy_x_pol_sep[0 ], copy_x_pol_sep[1 ]);
109+
110+ // Inner contraction: remove points belonging to the polar set by contracting
111+ // with the complement of the polar outer box
112+ x_cart_sep.inner = contract_polar_box_complement (ctc_polar, x_cart_sep.inner , x_pol_sep.outer );
113+
114+ return x_cart_sep;
115+ }
0 commit comments