Skip to content

CSG: phi-wedge support for ZSphere/Cylinder (bake startPhi/deltaPhi) #345

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ggalgoczi wants to merge 2 commits into
base: main
Choose a base branch
from
Open

CSG: phi-wedge support for ZSphere/Cylinder (bake startPhi/deltaPhi) #345

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7dba2ca
CSG: phi-wedge support for ZSphere/Cylinder (bake startPhi/deltaPhi)
ggalgoczi May 29, 2026
8c2414b
CSG: add phi-wedge leaf-intersection regression test (ZSphere + Cylin…
ggalgoczi May 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion CSG/CSGNode.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -422,7 +422,10 @@ void CSGNode::setAABBLocal()
{
float px, py, a, radius, z1, z2 ;
getParam( px, py, a, radius, z1, z2) ;
setAABB( px-radius, py-radius, z1, px+radius, py+radius, z2 );
// F6: px,py (cx,cy) are repurposed as startPhi,deltaPhi for phi-wedge
// cylinders (sn::Cylinder 5-arg), so they are NOT the centre. Cylinder is
// origin-centred -> AABB is simply +/-radius in x/y over [z1,z2].
setAABB(-radius, -radius, z1, radius, radius, z2);
}
else if( tc == CSG_DISC )
{
Expand Down
202 changes: 155 additions & 47 deletions CSG/csg_intersect_leaf_cylinder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,51 +33,143 @@ into a float4 and then pick from that ?
LEAF_FUNC
void intersect_leaf_cylinder( bool& valid_isect, float4& isect, const quad& q0, const quad& q1, const float t_min, const float3& ray_origin, const float3& ray_direction )
{
const float& r = q0.f.w ;
const float& z1 = q1.f.x ;
const float& z2 = q1.f.y ;
const float& ox = ray_origin.x ;
const float& oy = ray_origin.y ;
const float& oz = ray_origin.z ;
const float& vx = ray_direction.x ;
const float& vy = ray_direction.y ;
const float& vz = ray_direction.z ;

const float r2 = r*r ;
const float a = vx*vx + vy*vy ; // see CSG/sympy_cylinder.py
const float b = ox*vx + oy*vy ;
const float c = ox*ox + oy*oy - r2 ;

float t_near, t_far, disc, sdisc ;
robust_quadratic_roots_disqualifying(t_min, t_near, t_far, disc, sdisc, a, b, c); // Solving: a t^2 + 2 b t + c = 0
float z_near = oz+t_near*vz ;
float z_far = oz+t_far*vz ;

const float t_z1cap = (z1 - oz)/vz ;
const float r2_z1cap = (ox+t_z1cap*vx)*(ox+t_z1cap*vx) + (oy+t_z1cap*vy)*(oy+t_z1cap*vy) ;

const float t_z2cap = (z2 - oz)/vz ;
const float r2_z2cap = (ox+t_z2cap*vx)*(ox+t_z2cap*vx) + (oy+t_z2cap*vy)*(oy+t_z2cap*vy) ;
const float& r = q0.f.w;
const float& z1 = q1.f.x;
const float& z2 = q1.f.y;
// Optional phi-wedge clip: q0.f.x=startPhi, q0.f.y=deltaPhi (rad).
// deltaPhi==0 (or >=2π) means no clip — original full-cylinder behaviour.
const float startPhi = q0.f.x;
const float deltaPhi = q0.f.y;
const bool has_phi_clip = (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf);
const float& ox = ray_origin.x;
const float& oy = ray_origin.y;
const float& oz = ray_origin.z;
const float& vx = ray_direction.x;
const float& vy = ray_direction.y;
const float& vz = ray_direction.z;

const float r2 = r * r;
const float a = vx * vx + vy * vy; // see CSG/sympy_cylinder.py
const float b = ox * vx + oy * vy;
const float c = ox * ox + oy * oy - r2;

float t_near, t_far, disc, sdisc;
robust_quadratic_roots_disqualifying(t_min, t_near, t_far, disc, sdisc, a, b, c); // Solving: a t^2 + 2 b t + c = 0
float z_near = oz + t_near * vz;
float z_far = oz + t_far * vz;

const float t_z1cap = (z1 - oz) / vz;
const float r2_z1cap = (ox + t_z1cap * vx) * (ox + t_z1cap * vx) + (oy + t_z1cap * vy) * (oy + t_z1cap * vy);

const float t_z2cap = (z2 - oz) / vz;
const float r2_z2cap = (ox + t_z2cap * vx) * (ox + t_z2cap * vx) + (oy + t_z2cap * vy) * (oy + t_z2cap * vy);

#ifdef DEBUG
//printf("// t_z1cap %10.4f r2_z1cap %10.4f \n", t_z1cap, r2_z1cap );
//printf("// t_z2cap %10.4f r2_z2cap %10.4f \n", t_z2cap, r2_z2cap );
// printf("// t_z1cap %10.4f r2_z1cap %10.4f \n", t_z1cap, r2_z1cap );
// printf("// t_z2cap %10.4f r2_z2cap %10.4f \n", t_z2cap, r2_z2cap );
#endif

float t_cand = CUDART_INF_F ;
if( t_near > t_min && z_near > z1 && z_near < z2 && t_near < t_cand ) t_cand = t_near ;
if( t_far > t_min && z_far > z1 && z_far < z2 && t_far < t_cand ) t_cand = t_far ;
if( t_z1cap > t_min && r2_z1cap <= r2 && t_z1cap < t_cand ) t_cand = t_z1cap ;
if( t_z2cap > t_min && r2_z2cap <= r2 && t_z2cap < t_cand ) t_cand = t_z2cap ;

valid_isect = t_cand > t_min && t_cand < CUDART_INF_F ;
auto phi_in_wedge = [&](float xh, float yh) -> bool {
if (!has_phi_clip)
return true;
float phi = atan2f(yh, xh);
const float twoPi = 2.f * M_PIf;
while (phi < startPhi) phi += twoPi;
return (phi - startPhi) <= deltaPhi;
};
// Side-surface candidates also need to be inside the wedge.
bool t_near_ok = (z_near > z1 && z_near < z2) && phi_in_wedge(ox + t_near * vx, oy + t_near * vy);
bool t_far_ok = (z_far > z1 && z_far < z2) && phi_in_wedge(ox + t_far * vx, oy + t_far * vy);
// Cap candidates need their hit point's r AND phi inside.
bool t_z1cap_ok = (r2_z1cap <= r2) && phi_in_wedge(ox + t_z1cap * vx, oy + t_z1cap * vy);
bool t_z2cap_ok = (r2_z2cap <= r2) && phi_in_wedge(ox + t_z2cap * vx, oy + t_z2cap * vy);

// Flat radial wedge-wall candidates at phi=startPhi and phi=startPhi+deltaPhi.
// The two half-planes contain the z-axis. A hit must land within
// r_h^2 <= r^2, z1 < z_h < z2, and on the +r side of the wall
// (i.e. not the antipodal half-plane through the same line).
// Outward normal at startPhi: n1 = ( sin(startPhi), -cos(startPhi), 0)
// Outward normal at startPhi+deltaPhi: n2 = (-sin(endPhi), cos(endPhi), 0)
float t_phi1 = CUDART_INF_F, t_phi2 = CUDART_INF_F;
bool t_phi1_ok = false, t_phi2_ok = false;
float n1x = 0.f, n1y = 0.f, n2x = 0.f, n2y = 0.f;
#ifndef CSG_CYL_NO_PHI_WALLS
if (has_phi_clip)
{
const float endPhi = startPhi + deltaPhi;
const float s1 = sinf(startPhi), c1 = cosf(startPhi);
const float s2 = sinf(endPhi), c2 = cosf(endPhi);
n1x = s1;
n1y = -c1;
n2x = -s2;
n2y = c2;

// phi=startPhi plane: { x*s1 - y*c1 = 0 } solved for t
const float denom1 = vx * s1 - vy * c1;
if (fabsf(denom1) > 1.e-12f)
{
t_phi1 = -(ox * s1 - oy * c1) / denom1;
const float xh = ox + t_phi1 * vx;
const float yh = oy + t_phi1 * vy;
const float zh = oz + t_phi1 * vz;
const float r2_h = xh * xh + yh * yh;
const bool on_pos_r = (xh * c1 + yh * s1) >= 0.f; // not antipodal half
t_phi1_ok = on_pos_r && (r2_h <= r2) && (zh > z1) && (zh < z2);
}
// phi=endPhi plane: { x*s2 - y*c2 = 0 } solved for t
const float denom2 = vx * s2 - vy * c2;
if (fabsf(denom2) > 1.e-12f)
{
t_phi2 = -(ox * s2 - oy * c2) / denom2;
const float xh = ox + t_phi2 * vx;
const float yh = oy + t_phi2 * vy;
const float zh = oz + t_phi2 * vz;
const float r2_h = xh * xh + yh * yh;
const bool on_pos_r = (xh * c2 + yh * s2) >= 0.f;
t_phi2_ok = on_pos_r && (r2_h <= r2) && (zh > z1) && (zh < z2);
}
}
#endif // CSG_CYL_NO_PHI_WALLS

float t_cand = CUDART_INF_F;
if (t_near > t_min && t_near_ok && t_near < t_cand)
t_cand = t_near;
if (t_far > t_min && t_far_ok && t_far < t_cand)
t_cand = t_far;
if (t_z1cap > t_min && t_z1cap_ok && t_z1cap < t_cand)
t_cand = t_z1cap;
if (t_z2cap > t_min && t_z2cap_ok && t_z2cap < t_cand)
t_cand = t_z2cap;
if (t_phi1 > t_min && t_phi1_ok && t_phi1 < t_cand)
t_cand = t_phi1;
if (t_phi2 > t_min && t_phi2_ok && t_phi2 < t_cand)
t_cand = t_phi2;

valid_isect = t_cand > t_min && t_cand < CUDART_INF_F;
if(valid_isect)
{
bool sheet = ( t_cand == t_near || t_cand == t_far ) ;
isect.x = sheet ? (ox + t_cand*vx)/r : 0.f ;
isect.y = sheet ? (oy + t_cand*vy)/r : 0.f ;
isect.z = sheet ? 0.f : ( t_cand == t_z1cap ? -1.f : 1.f) ;
isect.w = t_cand ;
bool sheet = (t_cand == t_near || t_cand == t_far);
bool cap = (t_cand == t_z1cap || t_cand == t_z2cap);
bool phi1 = (t_cand == t_phi1);
if (sheet)
{
isect.x = (ox + t_cand * vx) / r;
isect.y = (oy + t_cand * vy) / r;
isect.z = 0.f;
}
else if (cap)
{
isect.x = 0.f;
isect.y = 0.f;
isect.z = (t_cand == t_z1cap ? -1.f : 1.f);
}
else // phi-wall
{
isect.x = phi1 ? n1x : n2x;
isect.y = phi1 ? n1y : n2y;
isect.z = 0.f;
}
isect.w = t_cand;
}
}

Expand Down Expand Up @@ -136,18 +228,34 @@ The SDF rules for CSG combinations::
LEAF_FUNC
float distance_leaf_cylinder( const float3& pos, const quad& q0, const quad& q1 )
{
const float radius = q0.f.w ;
const float z1 = q1.f.x ;
const float z2 = q1.f.y ; // z2 > z1
const float radius = q0.f.w;
const float z1 = q1.f.x;
const float z2 = q1.f.y; // z2 > z1
const float startPhi = q0.f.x;
const float deltaPhi = q0.f.y;

float sd_capslab = fmaxf( pos.z - z2 , z1 - pos.z );
float sd_infcyl = sqrtf( pos.x*pos.x + pos.y*pos.y ) - radius ;
float sd = fmaxf( sd_capslab, sd_infcyl );
float sd_capslab = fmaxf(pos.z - z2, z1 - pos.z);
float sd_infcyl = sqrtf(pos.x * pos.x + pos.y * pos.y) - radius;
float sd = fmaxf(sd_capslab, sd_infcyl);

if (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf)
{
float phi = atan2f(pos.y, pos.x);
const float twoPi = 2.f * M_PIf;
while (phi < startPhi) phi += twoPi;
float dphi = phi - startPhi;
if (dphi > deltaPhi)
{
float r_xy = sqrtf(pos.x * pos.x + pos.y * pos.y);
float d_to_edge = r_xy * fminf(dphi - deltaPhi, twoPi - dphi);
sd = fmaxf(sd, d_to_edge);
}
}

#ifdef DEBUG
printf("//distance_leaf_cylinder sd %10.4f \n", sd );
printf("//distance_leaf_cylinder sd %10.4f \n", sd);
#endif
return sd ;
return sd;
}


87 changes: 72 additions & 15 deletions CSG/csg_intersect_leaf_zsphere.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,18 +3,38 @@
LEAF_FUNC
float distance_leaf_zsphere(const float3& pos, const quad& q0, const quad& q1 )
{
float3 center = make_float3(q0.f);
float radius = q0.f.w;
// ZSphere is always centred at the local origin (G4Sphere doesn't carry
// a translation in its primitive — that's done by the parent transform).
// q0.f.x and q0.f.y therefore carry an OPTIONAL phi-wedge clip:
// q0.f.x = startPhi (rad), q0.f.y = deltaPhi (rad).
// deltaPhi == 0 (or ≥ 2π) means no clip — original ZSphere behaviour.
const float startPhi = q0.f.x;
const float deltaPhi = q0.f.y;
const float radius = q0.f.w;
const float2 zdelta = make_float2(q1.f);
const float z2 = center.z + zdelta.y ;
const float z1 = center.z + zdelta.x ;
const float z2 = zdelta.y;
const float z1 = zdelta.x;

float3 p = pos - center;
float sd_sphere = length(p) - radius ;
float sd_capslab = fmaxf( pos.z - z2 , z1 - pos.z );
const float3 p = pos;
float sd_sphere = length(p) - radius;
float sd_capslab = fmaxf(pos.z - z2, z1 - pos.z);

float sd = fmaxf( sd_capslab, sd_sphere ); // CSG intersect
return sd ;

if (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf)
{
float phi = atan2f(p.y, p.x);
const float twoPi = 2.f * M_PIf;
while (phi < startPhi) phi += twoPi;
float dphi = phi - startPhi;
if (dphi > deltaPhi)
{
float r_xy = sqrtf(p.x * p.x + p.y * p.y);
float d_to_edge = r_xy * fminf(dphi - deltaPhi, twoPi - dphi);
sd = fmaxf(sd, d_to_edge);
}
}
return sd;
}

/**
Expand Down Expand Up @@ -44,10 +64,24 @@ See : notes/issues/unexpected_zsphere_miss_from_inside_for_rays_that_would_be_ex
LEAF_FUNC
void intersect_leaf_zsphere(bool& valid_isect, float4& isect, const quad& q0, const quad& q1, const float& t_min, const float3& ray_origin, const float3& ray_direction )
{
const float3 center = make_float3(q0.f);
float3 O = ray_origin - center;
// ZSphere is always centred at the local origin; q0.f.x/q0.f.y carry an
// OPTIONAL phi-wedge clip (see distance_leaf_zsphere comment).
const float3 center = make_float3(0.f, 0.f, 0.f);
float3 O = ray_origin - center;
float3 D = ray_direction;
const float radius = q0.f.w;
#ifdef ZSPHERE_PHI_DEBUG
if (q0.f.y > 0.f && q0.f.y < 6.28f)
{
static int cnt = 0;
if (cnt < 5)
{
printf("// ZSphere phi: startPhi=%.4f deltaPhi=%.4f radius=%.2f z1=%.2f z2=%.2f\n",
q0.f.x, q0.f.y, q0.f.w, q1.f.x, q1.f.y);
cnt++;
}
}
#endif

float b = dot(O, D); // t of closest approach to sphere center
float c = dot(O, O)-radius*radius; // < 0. indicates ray_origin inside sphere
Expand Down Expand Up @@ -109,18 +143,41 @@ void intersect_leaf_zsphere(bool& valid_isect, float4& isect, const quad& q0, co

// hmm somehow is seems unclean to have to use both z and t language

float t_cand = t_min ;
// Optional phi-wedge clip on the sphere-surface candidates.
// q0.f.x = startPhi (rad), q0.f.y = deltaPhi (rad). deltaPhi == 0
// (or ≥ 2π) means no clip — the original two-cap ZSphere behaviour.
const float startPhi = q0.f.x;
const float deltaPhi = q0.f.y;
const bool has_phi_clip = (deltaPhi > 0.f && deltaPhi < 2.f * M_PIf);
auto phi_in_wedge = [&](float t) -> bool {
if (!has_phi_clip)
return true;
float x = O.x + t * D.x; // local-frame hit coords (O = origin − center)
float y = O.y + t * D.y;
float phi = atan2f(y, x);
const float twoPi = 2.f * M_PIf;
while (phi < startPhi) phi += twoPi;
float dphi = phi - startPhi;
return dphi <= deltaPhi;
};
bool t1sph_ok = (z1sph > zmin && z1sph <= zmax) && phi_in_wedge(t1sph);
bool t2sph_ok = (z2sph > zmin && z2sph <= zmax) && phi_in_wedge(t2sph);

float t_cand = t_min;
if(sdisc > 0.f)
{

#ifdef DEBUG_RECORD
//std::raise(SIGINT);
// std::raise(SIGINT);
#endif

if( t1sph > t_min && z1sph > zmin && z1sph <= zmax ) t_cand = t1sph ; // t1sph qualified and t1cap disabled or disqualified -> t1sph
else if( t1cap > t_min ) t_cand = t1cap ; // t1cap qualifies -> t1cap
if (t1sph > t_min && t1sph_ok)
t_cand = t1sph; // t1sph qualified and t1cap disabled or disqualified -> t1sph
else if (t1cap > t_min)
t_cand = t1cap; // t1cap qualifies -> t1cap
else if( t2cap > t_min ) t_cand = t2cap ; // t2cap qualifies -> t2cap
else if( t2sph > t_min && z2sph > zmin && z2sph <= zmax) t_cand = t2sph ; // t2sph qualifies and t2cap disabled or disqialified -> t2sph
else if (t2sph > t_min && t2sph_ok)
t_cand = t2sph; // t2sph qualifies and t2cap disabled or disqialified -> t2sph
}

valid_isect = t_cand > t_min ;
Expand Down
3 changes: 2 additions & 1 deletion CSG/tests/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,8 @@ set(TEST_SOURCES

CSGCopyTest.cc

intersect_leaf_phicut_test.cc
intersect_leaf_phicut_test.cc
intersect_leaf_phi_wedge_test.cc
intersect_leaf_thetacut_test.cc
intersect_leaf_box3_test.cc
intersect_leaf_cylinder_test.cc
Expand Down
Loading
Loading