// Author: Colin Andrews
// Institution: Maxis Emeryville
// email: candrews@maxis.com
namespace
ColinAndrews_ElectronicArtsMaxis_sub7
{
typedef unsigned int size_t;
/// quick_sort
stuff - All massaged from Paul Pedriana's EASTL implementations
static const int
kQuickSortLimit = 28;
// TODO: template
based on value_type unnecessary - Always PointC
// it would be interesting to compare the two
if there's a performance
// difference in the generated code
template
<typename T, typename
Compare>
inline const T& median(const
T& a, const T& b, const T& c, Compare compare)
{
if(compare(a,
b))
{
if(compare(b,
c))
return
b;
else
if(compare(a, c))
return
c;
else
return
a;
}
else if(compare(a, c))
return
a;
else if(compare(b, c))
return
c;
return
b;
}
inline
size_t Log2(size_t n)
{
int i;
for(i =
0; n; ++i)
n >>= 1;
return
i - 1;
}
template
<typename value_type>
inline void simple_swap(value_type& a, value_type&
b)
{
value_type tmp(a);
a = b;
b = tmp;
}
template
<typename value_type, typename Compare>
inline
value_type* get_partition(value_type* first, value_type* last, value_type
pivotValue, Compare compare)
{
for(; ; ++first)
{
while(compare(*first,
pivotValue))
++first;
--last;
while(compare(pivotValue,
*last))
--last;
if(first
>= last)
return
first;
simple_swap(*first, *last);
}
}
template
<typename value_type, typename StrictWeakOrdering>
void
insertion_sort(value_type* first, value_type* last, StrictWeakOrdering compare)
{
if(first
!= last)
{
value_type* iCurrent;
value_type* iNext;
value_type* iSorted = first;
for(++iSorted;
iSorted != last; ++iSorted)
{
const
value_type temp(*iSorted);
iNext = iCurrent = iSorted;
for(--iCurrent;
(iNext != first) && compare(temp, *iCurrent); --iNext, --iCurrent)
*iNext = *iCurrent;
*iNext = temp;
}
}
} //
insertion_sort
template
<typename value_type, typename Compare>
inline void insertion_sort_simple(value_type* first,
value_type* last, Compare compare)
{
for(value_type*
current = first; current != last; ++current)
{
value_type* end(current);
value_type* prev(current);
const
value_type value(*current);
for(--prev;
compare(value, *prev); --end, --prev)
*end = *prev;
*end = value;
}
}
template
<typename value_type, typename Compare>
inline void quick_sort_impl(value_type* first, value_type*
last, size_t kRecursionCount, Compare compare)
{
while(((last
- first) > kQuickSortLimit) && (kRecursionCount > 0))
{
value_type*
position(get_partition<value_type, Compare>(first, last,
median<value_type, Compare>(*first, *(first + (last - first) / 2), *(last
- 1), compare), compare));
quick_sort_impl<value_type,
Compare>(position, last, --kRecursionCount, compare);
last = position;
}
if(kRecursionCount
== 0)
insertion_sort_simple<value_type, Compare>(first, last, compare); // replaced partial_sort
}
template
<typename value_type, typename Compare>
void
quick_sort(value_type* first, value_type* last, Compare compare)
{
if(first
!= last)
{
quick_sort_impl<value_type,
Compare>(first, last, 2 * Log2(size_t(last - first)), compare);
if((last
- first) > (size_t)kQuickSortLimit)
{
insertion_sort<value_type,
Compare>(first, first + kQuickSortLimit, compare);
insertion_sort_simple<value_type, Compare>(first +
kQuickSortLimit, last, compare);
}
else
insertion_sort<value_type,
Compare>(first, last, compare);
}
}
// my distance
functions
inline float xzdist2(PointC *pA, PointC *pB)
{
register
float d1(pA->x - pB->x);
register
float d2(pA->z - pB->z);
return
d2*d2 + d1*d1;
}
inline float yd2(PointC *pA, PointC *pB)
{
register
float d(pA->y - pB->y);
return
d*d;
}
inline float xydist2(PointC *pA, PointC *pB)
{
register
float d1(pA->x - pB->x);
register
float d2(pA->y - pB->y);
return
d2*d2 + d1*d1;
}
inline float zd2(PointC *pA, PointC *pB)
{
register
float d(pA->z - pB->z);
return
d*d;
}
inline float yzdist2(PointC *pA, PointC *pB)
{
register
float d1(pA->y - pB->y);
register
float d2(pA->z - pB->z);
return
d2*d2 + d1*d1;
}
inline float xd2(PointC *pA, PointC *pB)
{
register
float d(pA->x - pB->x);
return
d*d;
}
struct
compareX {
inline bool operator()(const PointC& a, const
PointC& b)
{
return
a.x < b.x;
}
};
struct
compareY {
inline bool operator()(const PointC& a, const
PointC& b)
{
return
a.y < b.y;
}
};
struct
compareZ {
inline bool operator()(const PointC& a, const
PointC& b)
{
return
b.z < a.z;
}
};
///////////////////////////////////////////////////////////////////////////////
// BoundaryCheck
// Checks for a closer pair that spans the
boundary between two regions
// of a y sorted cloud where the closest
distance in either region has
// already been determined to be fMin
inline float BoundaryCheckX(float
minDist, PointC *pBase, PointC *pMidPoint, PointC *pEnd) {
PointC *pUpper(pMidPoint);
PointC *pLower(pMidPoint-1);
register
float
y2(xd2(pLower,pUpper));
while
(y2<minDist) { // outer loop bumps lower
while
(y2<minDist) { // inner loop bumps upper
register
float dist(y2+yzdist2(pLower,pUpper));
if
(dist < minDist) {
minDist=dist;
}
++pUpper;
if
(pUpper==pEnd) {
break;
}
y2 = xd2(pLower, pUpper);
}
//
decrement lower, reset upper
--pLower;
if
(pLower<pBase) {
break;
}
pUpper = pMidPoint;
y2 = xd2(pLower, pUpper);
}
return
minDist;
}
inline float BoundaryCheckY(float
minDist, PointC *pBase, PointC *pMidPoint, PointC *pEnd) {
PointC *pUpper(pMidPoint);
PointC *pLower(pMidPoint-1);
register
float
y2(yd2(pLower,pUpper));
while
(y2<minDist) { // outer loop bumps lower
while
(y2<minDist) { // inner loop bumps upper
register
float dist(y2+xzdist2(pLower,pUpper));
if
(dist < minDist) {
minDist=dist;
}
++pUpper;
if
(pUpper==pEnd) {
break;
}
y2 = yd2(pLower, pUpper);
}
//
decrement lower, reset upper
--pLower;
if
(pLower<pBase) {
break;
}
pUpper = pMidPoint;
y2 = yd2(pLower, pUpper);
}
return
minDist;
}
inline float BoundaryCheckZ(float
minDist, PointC *pBase, PointC *pMidPoint, PointC *pEnd) {
PointC *pUpper(pMidPoint);
PointC *pLower(pMidPoint-1);
register
float
y2(zd2(pLower,pUpper));
while
(y2<minDist) { // outer loop bumps lower
while
(y2<minDist) { // inner loop bumps upper
register
float dist(y2+xydist2(pLower,pUpper));
if
(dist < minDist) {
minDist=dist;
}
++pUpper;
if
(pUpper==pEnd) {
break;
}
y2 = zd2(pLower, pUpper);
}
//
decrement lower, reset upper
--pLower;
if
(pLower<pBase) {
break;
}
pUpper = pMidPoint;
y2 = zd2(pLower, pUpper);
}
return
minDist;
}
/******************************************************************************/
/*! GetMinDist
\brief Returns the distance between the two
"Nearest Neighbors" in the Cloud.
\param pCloud
- CloudC pointer
\return - float.
I've had a copy of Preparata
& Shamos since I was 19, and that was a
really long time ago. (it was
brand new) When I read this challenge I pulled
out my copy (& de Berg et
al) to look for ideas.
The original 'optimal' solution
by Mike Shamos is a mess. Since it is a O(NlogN)
solution, O(N) operations are
treated as being free. With the sizes of point
clouds we are to expect for this
problem, anything that walks all the points
could easily make an
"optimal" solution take much longer than an efficient,
non-optimal solution. I think
that this is a really interesting choice of a
problem because of that
dichotomy between academic optimality and real world
hacker speed freak efficiency.
Like my other solutions, this
one takes the idea of sorting on a major axis to
limit the scope of the search.
This one also attempts to mimic the divide and
conquer search that is used in
Shamos et. al., but folding a recursive subdivision
into an iterative one. In the
first iteration of the loop each even odd pair
is checked for min distance.
Successive iterations merge the results into succesively
larger regions that have been
searched. When merging regions, it is only necessary
to check for a pair that spans
the boundary between two regions that is closer
than current minimum. Shamos
does a lot of fancy stuff to limit the scope of
the merge step. I just made a
pretty straightforward iterative search again limited
by the current working min
distance.
My test results show that it is
largely random which axis is best to sort on
for a given model data set. The
quicksort takes the same amount of time for all
the models I tested with. The
best results happen when a close pair is found
early so the remaining search
doesn't have to consider as many points. It also
is random whether this divide
and conquer search is faster than the sort-limited
brute force search used in my
other submissions.
My remaining strategy is the
shotgun approach: Each of my remaining submissions
sorts on a different axis and
uses a different search algorithm in the hope
that one of these is the best
match for the sample data you test us against. No
one combination of these was the
best over all my test input data sets.
This one sorts on the x axis and
uses the divide and conquer search described
above.
*/
/******************************************************************************/
float
GetMinDist(CloudC *pCloud)
{
PointC * pBuffer = pCloud->GetPointBuffer();
unsigned
long uiNumPoints = pCloud->GetNumPoints();
PointC * pBufferEnd = pBuffer +
uiNumPoints;
compareY compare;
quick_sort(pBuffer, pBufferEnd, compare);
float
minDist = FLT_MAX;
for (unsigned long
nStep=2; nStep<uiNumPoints; nStep <<=1) {
PointC * pStart = pBuffer;
PointC * pMid = pBuffer + (nStep>>1);
PointC * pEnd = pBuffer + nStep;
for
(; pEnd <= pBufferEnd; pStart+=nStep, pMid+=nStep, pEnd+=nStep) {
float
nMin = BoundaryCheckY(minDist, pStart, pMid, pEnd);
if
(nMin<minDist) {
minDist = nMin;
}
}
//
remainder - numPoints is not a multiple of current power of 2
if
(pBufferEnd>pMid) {
float
nMin = BoundaryCheckY(minDist, pStart, pMid, pBufferEnd);
if
(nMin<minDist) {
minDist = nMin;
}
}
}
// return the
distance between the two closest points
return
sqrtf(minDist); // not actually correct
}
} // end namespace