Skip to content
GitLab
Commit c0edcb8e authored by LEFEBVREJP email's avatar LEFEBVREJP email
Browse files

Adding Grid::inside(const Ray&) and format cleanup.

parent e4e7f588
Pipeline #7028 skipped
Hide whitespace changes
Inline Side-by-side
radixgeometry/grid.cc
View file @ c0edcb8e
......@@ -49,189 +49,195 @@ Grid::Grid(const Grid& orig)
}
Grid::Grid(Point3D lowerLeftCorner
, Real deltaX, Real deltaY, Real deltaZ
, int numX, int numY, int numZ)
: mTranslation(lowerLeftCorner)
, mXPlanes(numX+1)
, mYPlanes(numY+1)
, mZPlanes(numZ+1)
, Real deltaX, Real deltaY, Real deltaZ
, int numX, int numY, int numZ)
: mTranslation(lowerLeftCorner)
, mXPlanes(numX+1)
, mYPlanes(numY+1)
, mZPlanes(numZ+1)
{
radix_check(mXPlanes.size() >= 2);
radix_check(mYPlanes.size() >= 2);
radix_check(mZPlanes.size() >= 2);
//
// Initialize first planes with lower left corner
mXPlanes[0] = 0;
mYPlanes[0] = 0;
mZPlanes[0] = 0;
//
// auto populate grid with delta#
for(size_t i = 1; i < mXPlanes.size(); ++i)
{
mXPlanes[i] = mXPlanes[i-1] + deltaX;
}
for(size_t i = 1; i < mYPlanes.size(); ++i)
{
mYPlanes[i] = mYPlanes[i-1] + deltaY;
}
for(size_t i = 1; i < mZPlanes.size(); ++i)
{
radix_check(mXPlanes.size() >= 2);
radix_check(mYPlanes.size() >= 2);
radix_check(mZPlanes.size() >= 2);
//
// Initialize first planes with lower left corner
mXPlanes[0] = 0;
mYPlanes[0] = 0;
mZPlanes[0] = 0;
//
// auto populate grid with delta#
for(size_t i = 1; i < mXPlanes.size(); ++i)
{
mXPlanes[i] = mXPlanes[i-1] + deltaX;
}
for(size_t i = 1; i < mYPlanes.size(); ++i)
{
mYPlanes[i] = mYPlanes[i-1] + deltaY;
}
for(size_t i = 1; i < mZPlanes.size(); ++i)
{
mZPlanes[i] = mZPlanes[i-1] + deltaZ;
}
mZPlanes[i] = mZPlanes[i-1] + deltaZ;
}
}
Grid::Grid(const std::vector<Real>&xplanes
, const std::vector<Real>&yplanes
, const std::vector<Real>&zplanes)
Grid::Grid(const std::vector<Real>&xplanes
, const std::vector<Real>&yplanes
, const std::vector<Real>&zplanes)
: mXPlanes(xplanes)
, mYPlanes(yplanes)
, mZPlanes(zplanes) {
radix_check(xplanes.size() >= 2);
radix_check(yplanes.size() >= 2);
radix_check(zplanes.size() >= 2);
}
radix_check(xplanes.size() >= 2);
radix_check(yplanes.size() >= 2);
radix_check(zplanes.size() >= 2);
}
Grid::~Grid() {
}
Grid::~Grid() {
}
void Grid::arrayCoordinates(Real &xo, Real &yo, Real &zo) const
{
xo+=mTranslation.x;
yo+=mTranslation.y;
zo+=mTranslation.z;
}
void Grid::arrayCoordinates(Real &xo, Real &yo, Real &zo) const
{
xo+=mTranslation.x;
yo+=mTranslation.y;
zo+=mTranslation.z;
}
void Grid::origin(int xi, int yi, int zi, Point3D &o) const
{
radix_require(!mXPlanes.empty());
radix_require(!mYPlanes.empty());
radix_require(!mZPlanes.empty());
radix_check(xi < int(mXPlanes.size())-1);
radix_check(yi < int(mYPlanes.size())-1);
radix_check(zi < int(mZPlanes.size())-1);
void Grid::origin(int xi, int yi, int zi, Point3D &o) const
{
radix_require(!mXPlanes.empty());
radix_require(!mYPlanes.empty());
radix_require(!mZPlanes.empty());
radix_check(xi < int(mXPlanes.size())-1);
radix_check(yi < int(mYPlanes.size())-1);
radix_check(zi < int(mZPlanes.size())-1);
Real x1 = mXPlanes[xi+1];
Real x0 = mXPlanes[xi];
Real y1 = mYPlanes[yi+1];
Real y0 = mYPlanes[yi];
Real z1 = mZPlanes[zi+1];
Real z0 = mZPlanes[zi];
o.x = x0 + (0.5 * (x1 - x0));
o.y = y0 + (0.5 * (y1 - y0));
o.z = z0 + (0.5 * (z1 - z0));
}
std::string Grid::toString() const {
std::stringstream stream;
stream << "Grid x0="
<< mXPlanes.front()
<< ",x1=" << mXPlanes.back()
<< " (" << mXPlanes.size() << ")"
<< " y0=" << mYPlanes.front()
<< ",y1=" << mYPlanes.back()
<< " (" << mYPlanes.size() << ")"
<< " z0=" << mZPlanes.front()
<< ",z1=" << mZPlanes.back()
<< " (" << mZPlanes.size() << ")";
return stream.str();
}
Real x1 = mXPlanes[xi+1];
Real x0 = mXPlanes[xi];
Real y1 = mYPlanes[yi+1];
Real y0 = mYPlanes[yi];
Real z1 = mZPlanes[zi+1];
Real z0 = mZPlanes[zi];
o.x = x0 + (0.5 * (x1 - x0));
o.y = y0 + (0.5 * (y1 - y0));
o.z = z0 + (0.5 * (z1 - z0));
}
std::string Grid::toString() const {
std::stringstream stream;
stream << "Grid x0="
<< mXPlanes.front()
<< ",x1=" << mXPlanes.back()
<< " (" << mXPlanes.size() << ")"
<< " y0=" << mYPlanes.front()
<< ",y1=" << mYPlanes.back()
<< " (" << mYPlanes.size() << ")"
<< " z0=" << mZPlanes.front()
<< ",z1=" << mZPlanes.back()
<< " (" << mZPlanes.size() << ")";
return stream.str();
}
bool Grid::inside(Real x, Real y, Real z)const {
radix_require(mXPlanes.size() >= 2);
radix_require(mYPlanes.size() >= 2);
radix_require(mZPlanes.size() >= 2);
arrayCoordinates(x,y,z);
bool Grid::inside(const Ray &r) const
{
return inside(r.o.x, r.o.y, r.o.z);
}
bool in = x + halfKEpsilon >= mXPlanes.front()
&& x - halfKEpsilon <= mXPlanes.back()
&& y + halfKEpsilon >= mYPlanes.front()
&& y - halfKEpsilon <= mYPlanes.back()
&& z + halfKEpsilon >= mZPlanes.front()
&& z - halfKEpsilon <= mZPlanes.back();
// radix_tagged_line(toString()
// <<", tested against x="<<x<<",y="<<y<<",z="<<z
// <<", inside ? "<<(in?"yes":"no"));
return in;
}
bool Grid::inside(Real x, Real y, Real z) const
{
radix_require(mXPlanes.size() >= 2);
radix_require(mYPlanes.size() >= 2);
radix_require(mZPlanes.size() >= 2);
arrayCoordinates(x,y,z);
bool Grid::start(const Ray &r, int &xi, int &yi, int &zi) const
{
return start(r.o.x, r.o.y, r.o.z, r.d.x, r.d.y, r.d.z, xi, yi, zi);
}
/**
bool in = x + halfKEpsilon >= mXPlanes.front()
&& x - halfKEpsilon <= mXPlanes.back()
&& y + halfKEpsilon >= mYPlanes.front()
&& y - halfKEpsilon <= mYPlanes.back()
&& z + halfKEpsilon >= mZPlanes.front()
&& z - halfKEpsilon <= mZPlanes.back();
// radix_tagged_line(toString()
// <<", tested against x="<<x<<",y="<<y<<",z="<<z
// <<", inside ? "<<(in?"yes":"no"));
return in;
}
bool Grid::start(const Ray &r, int &xi, int &yi, int &zi) const
{
return start(r.o.x, r.o.y, r.o.z, r.d.x, r.d.y, r.d.z, xi, yi, zi);
}
/**
* Determine the starting array index for the given point
* @return true, if and only if the point starts inside the array
*/
bool Grid::start(Real xo, Real yo, Real zo,
Real xd, Real yd, Real zd,
int& xi, int& yi, int&zi
) const {
radix_tagged_line(Point3D(xo,yo,zo)<<" translated into by "<<translation<<"...");
arrayCoordinates(xo,yo,zo); // put world coordinates into array
radix_tagged_line(Point3D(xo,yo,zo));
Real tx = xd < 0 ? xo - halfKEpsilon : (xd > 0 ? xo + halfKEpsilon : xo);
Real ty = yd < 0 ? yo - halfKEpsilon : (yd > 0 ? yo + halfKEpsilon : yo);
Real tz = zd < 0 ? zo - halfKEpsilon : (zd > 0 ? zo + halfKEpsilon : zo);
std::vector<Real>::const_iterator xitr = std::lower_bound(mXPlanes.begin(), mXPlanes.end(), tx);
std::vector<Real>::const_iterator yitr = std::lower_bound(mYPlanes.begin(), mYPlanes.end(), ty);
std::vector<Real>::const_iterator zitr = std::lower_bound(mZPlanes.begin(), mZPlanes.end(), tz);
xi = xitr - mXPlanes.begin();
yi = yitr - mYPlanes.begin();
zi = zitr - mZPlanes.begin();
radix_tagged_line("Start lower_bound indices are "<<xi<<","<<yi<<","<<zi
<<", xo="<<xo<<",yo="<<yo<<",zo="<<zo
<<", xd="<<xd<<",yd="<<yd<<",zd="<<zd
<<", "<<toString());
bool foundx = xi != int(mXPlanes.size());
bool foundy = yi != int(mYPlanes.size());
bool foundz = zi != int(mZPlanes.size());
bool Grid::start(Real xo, Real yo, Real zo,
Real xd, Real yd, Real zd,
int& xi, int& yi, int&zi
) const {
radix_tagged_line(Point3D(xo,yo,zo)<<" translated into by "<<translation<<"...");
arrayCoordinates(xo,yo,zo); // put world coordinates into array
radix_tagged_line(Point3D(xo,yo,zo));
Real tx = xd < 0 ? xo - halfKEpsilon : (xd > 0 ? xo + halfKEpsilon : xo);
Real ty = yd < 0 ? yo - halfKEpsilon : (yd > 0 ? yo + halfKEpsilon : yo);
Real tz = zd < 0 ? zo - halfKEpsilon : (zd > 0 ? zo + halfKEpsilon : zo);
std::vector<Real>::const_iterator xitr = std::lower_bound(mXPlanes.begin(), mXPlanes.end(), tx);
std::vector<Real>::const_iterator yitr = std::lower_bound(mYPlanes.begin(), mYPlanes.end(), ty);
std::vector<Real>::const_iterator zitr = std::lower_bound(mZPlanes.begin(), mZPlanes.end(), tz);
xi = xitr - mXPlanes.begin();
yi = yitr - mYPlanes.begin();
zi = zitr - mZPlanes.begin();
radix_tagged_line("Start lower_bound indices are "<<xi<<","<<yi<<","<<zi
<<", xo="<<xo<<",yo="<<yo<<",zo="<<zo
<<", xd="<<xd<<",yd="<<yd<<",zd="<<zd
<<", "<<toString());
bool foundx = xi != int(mXPlanes.size());
bool foundy = yi != int(mYPlanes.size());
bool foundz = zi != int(mZPlanes.size());
bool isxplane = foundx && std::abs(mXPlanes[xi] - xo) < kEpsilon;
bool isyplane = foundy && std::abs(mYPlanes[yi] - yo) < kEpsilon;
bool iszplane = foundz && std::abs(mZPlanes[zi] - zo) < kEpsilon;
// radix_tagged_line("isx ? "<<(isxplane?"yes":"no")
// <<", isy ? "<<(isyplane?"yes":"no")
// <<", isz ? "<<(iszplane?"yes":"no"));
//
// determine if any of the axis are prior to their respective first plane
//
bool prior = (!isxplane && xi == 0)
|| (!isyplane && yi == 0)
|| (!iszplane && zi == 0);
bool isxplane = foundx && std::abs(mXPlanes[xi] - xo) < kEpsilon;
bool isyplane = foundy && std::abs(mYPlanes[yi] - yo) < kEpsilon;
bool iszplane = foundz && std::abs(mZPlanes[zi] - zo) < kEpsilon;
// radix_tagged_line("isx ? "<<(isxplane?"yes":"no")
// <<", isy ? "<<(isyplane?"yes":"no")
// <<", isz ? "<<(iszplane?"yes":"no"));
//
// determine if any of the axis are prior to their respective first plane
//
bool prior = (!isxplane && xi == 0)
|| (!isyplane && yi == 0)
|| (!iszplane && zi == 0);
bool isfirstxplane = isxplane && xi == 0;
bool isfirstyplane = isyplane && yi == 0;
bool isfirstzplane = iszplane && zi == 0;
bool islastxplane = isxplane && xi == int(mXPlanes.size()) - 1;
bool islastyplane = isyplane && yi == int(mYPlanes.size()) - 1;
bool islastzplane = iszplane && zi == int(mZPlanes.size()) - 1;
bool isfirstxplane = isxplane && xi == 0;
bool isfirstyplane = isyplane && yi == 0;
bool isfirstzplane = iszplane && zi == 0;
bool islastxplane = isxplane && xi == int(mXPlanes.size()) - 1;
bool islastyplane = isyplane && yi == int(mYPlanes.size()) - 1;
bool islastzplane = iszplane && zi == int(mZPlanes.size()) - 1;
bool exiting = (isfirstxplane && xd < 0)// on first xplane heading away
|| (isfirstyplane && yd < 0)// on first yplane heading away
|| (isfirstzplane && zd < 0)// on first zplane heading away
|| (islastxplane && xd > 0)// on last xplane heading away
|| (islastyplane && yd > 0)// on last yplane heading away
|| (islastzplane && zd > 0)// on last zplane heading away
;
bool exiting = (isfirstxplane && xd < 0)// on first xplane heading away
|| (isfirstyplane && yd < 0)// on first yplane heading away
|| (isfirstzplane && zd < 0)// on first zplane heading away
|| (islastxplane && xd > 0)// on last xplane heading away
|| (islastyplane && yd > 0)// on last yplane heading away
|| (islastzplane && zd > 0)// on last zplane heading away
;
if (!isxplane || xd) {
xi -= 1;
}
if (!isyplane || yd < 0) {
yi -= 1;
}
if (!iszplane || zd < 0) {
zi -= 1;
}
// radix_tagged_line("start - xi="<<xi<<",yi="<<yi<<",zi="<<zi
// <<" out of xcount="<<xplanes.size()<<",ycount="<<yplanes.size()
// <<",zcount="<<zplanes.size());
if (!isxplane || xd) {
xi -= 1;
}
if (!isyplane || yd < 0) {
yi -= 1;
}
if (!iszplane || zd < 0) {
zi -= 1;
}
// radix_tagged_line("start - xi="<<xi<<",yi="<<yi<<",zi="<<zi
// <<" out of xcount="<<xplanes.size()<<",ycount="<<yplanes.size()
// <<",zcount="<<zplanes.size());
return foundx && foundy && foundz && !prior && !exiting;
} // end of start
return foundx && foundy && foundz && !prior && !exiting;
} // end of start
/**
/**
* Determine the next array index that the given ray will move into
* @param xo - the x origin
* @param yo - the y origin
......@@ -245,63 +251,63 @@ Grid::Grid(Point3D lowerLeftCorner
* @param t - the track length to populate with the distance to the next index (if any)
* NOTE: It is assumed that the boundary conditions are managed by the caller
*/
void Grid::next(Real xo, Real yo, Real zo,
Real xd, Real yd, Real zd,
int& xi, int& yi, int& zi,
Real& t
) const {
// assume boundary conditions managed outside
arrayCoordinates(xo,yo,zo); // put world coordinates into array
Real tx = maxReal;
Real ty = maxReal;
Real tz = maxReal;
int cx = xi;
if (xd > 0) {
cx = xi + 1;
tx = (mXPlanes[cx] - xo) / xd;
} else if (xd < 0) {
cx = xi - 1;
tx = (mXPlanes[(mXPlanes[xi]<xo?xi:cx)] - xo) / xd;
}
int cy = yi;
if (yd > 0) {
cy = yi + 1;
ty = (mYPlanes[cy] - yo) / yd;
} else if (yd < 0) {
cy = yi - 1;
ty = (mYPlanes[(mYPlanes[yi]<yo?yi:cy)] - yo) / yd;
}
int cz = zi;
if (zd > 0) {
cz = zi + 1;
tz = (mZPlanes[cz] - zo) / zd;
} else if (zd < 0) {
cz = zi - 1;
tz = (mZPlanes[(mZPlanes[zi]<zo?zi:cz)] - zo) / zd;
}
// radix_tagged_line("xo=" << xo << ",yo=" << yo << ",zo=" << zo
// << ",xd=" << xd << ",yd=" << yd << ",zd=" << zd
// << ",tx=" << tx << ",ty=" << ty << ",tz=" << tz);
bool txety = std::abs(tx - ty) < halfKEpsilon;
bool txetz = std::abs(tx - tz) < halfKEpsilon;
bool tyetz = std::abs(ty - tz) < halfKEpsilon;
void Grid::next(Real xo, Real yo, Real zo,
Real xd, Real yd, Real zd,
int& xi, int& yi, int& zi,
Real& t
) const {
// assume boundary conditions managed outside
arrayCoordinates(xo,yo,zo); // put world coordinates into array
Real tx = maxReal;
Real ty = maxReal;
Real tz = maxReal;
int cx = xi;
if (xd > 0) {
cx = xi + 1;
tx = (mXPlanes[cx] - xo) / xd;
} else if (xd < 0) {
cx = xi - 1;
tx = (mXPlanes[(mXPlanes[xi]<xo?xi:cx)] - xo) / xd;
}
int cy = yi;
if (yd > 0) {
cy = yi + 1;
ty = (mYPlanes[cy] - yo) / yd;
} else if (yd < 0) {
cy = yi - 1;
ty = (mYPlanes[(mYPlanes[yi]<yo?yi:cy)] - yo) / yd;
}
int cz = zi;
if (zd > 0) {
cz = zi + 1;
tz = (mZPlanes[cz] - zo) / zd;
} else if (zd < 0) {
cz = zi - 1;
tz = (mZPlanes[(mZPlanes[zi]<zo?zi:cz)] - zo) / zd;
}
// radix_tagged_line("xo=" << xo << ",yo=" << yo << ",zo=" << zo
// << ",xd=" << xd << ",yd=" << yd << ",zd=" << zd
// << ",tx=" << tx << ",ty=" << ty << ",tz=" << tz);
bool txety = std::abs(tx - ty) < halfKEpsilon;
bool txetz = std::abs(tx - tz) < halfKEpsilon;
bool tyetz = std::abs(ty - tz) < halfKEpsilon;
bool txlty = tx < ty;
bool txltz = tx < tz;
bool tyltz = ty < tz;
t = maxReal; // make sure t is prepared for our algorithm below
bool txlty = tx < ty;
bool txltz = tx < tz;
bool tyltz = ty < tz;
t = maxReal; // make sure t is prepared for our algorithm below
if ((txlty || txety) && (txltz || txetz)) {
xi = cx;
t = tx;
}
if ((txety || !txlty) && (tyltz || tyetz)) {
yi = cy;
t = std::min(ty, t);
}
if (!txltz && !tyltz) {
zi = cz;
t = std::min(tz, t);
}
}// end of next
if ((txlty || txety) && (txltz || txetz)) {
xi = cx;
t = tx;
}
if ((txety || !txlty) && (tyltz || tyetz)) {
yi = cy;
t = std::min(ty, t);
}
if (!txltz && !tyltz) {
zi = cz;
t = std::min(tz, t);
}
}// end of next
} // end of namespace
radixgeometry/grid.hh
View file @ c0edcb8e
......@@ -51,6 +51,13 @@ public:
return mZPlanes.size();
}
std::string toString() const;
/**
* @brief inside convenience method for inside(Real, Real, Real.
* @param r const Ray&
* @return true if the ray is inside or on the grid
*/
bool inside(const Ray& r) const;
/**
* Determine if the given x,y,z is inside (or on) the array
*/
......
radixgeometry/tests/tstGrid.cc
View file @ c0edcb8e
......@@ -39,6 +39,11 @@ TEST(Grid, Constructor)
zplanes[5] = 4500;
EXPECT_EQ(zplanes, grid.zplanes());
Ray ray(Point3D(-116.15,37.16667,1000)
, Vector3D(-116.14, 37.165, 1000));
ray.d.normalize();
std::cout << ray.toString() << std::endl;
}
TEST(Grid, Uniform)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment