Source

nrnr / nrnr.c

Peter Shinners 4a6c0a4 







Peter Shinners eccbc3b 

Peter Shinners 4a6c0a4 



Peter Shinners 2cfda6f 
















Peter Shinners 4a6c0a4 
































Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 

Peter Shinners 2cfda6f 



Peter Shinners 4a6c0a4 

Peter Shinners 2cfda6f 
Peter Shinners 4a6c0a4 

Peter Shinners 2cfda6f 
Peter Shinners 4a6c0a4 

Peter Shinners 2cfda6f 
Peter Shinners 4a6c0a4 

Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 
Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 

Peter Shinners 2cfda6f 
Peter Shinners 4a6c0a4 
























Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 

Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 



















Peter Shinners 1a75ba3 
Peter Shinners 4a6c0a4 

Peter Shinners 2cfda6f 


Peter Shinners 1a75ba3 

Peter Shinners 2cfda6f 





Peter Shinners 1a75ba3 
Peter Shinners 2cfda6f 


Peter Shinners 1a75ba3 
Peter Shinners 2cfda6f 

Peter Shinners eccbc3b 

Peter Shinners 2cfda6f 
Peter Shinners 1a75ba3 



Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 
Peter Shinners 2cfda6f 
Peter Shinners 1a75ba3 
Peter Shinners 4a6c0a4 

Peter Shinners 2cfda6f 
Peter Shinners 4a6c0a4 





































































Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 
Peter Shinners eccbc3b 


Peter Shinners 4a6c0a4 
Peter Shinners eccbc3b 


Peter Shinners 2cfda6f 

Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 

Peter Shinners eccbc3b 


Peter Shinners 4a6c0a4 
Peter Shinners 2cfda6f 





Peter Shinners 4a6c0a4 



Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 
Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 


Peter Shinners eccbc3b 


Peter Shinners 4a6c0a4 




Peter Shinners eccbc3b 
Peter Shinners 4a6c0a4 

#include <string.h>
#include <alloca.h>
#include <stdio.h> // temp
#include "nrnr.h"


// start is inclusive, end is exclusive
struct range {
    size_t start;
    size_t end;
};



static void SortIndices(size_t* array, int numPoints, const float* points, size_t axis, size_t dimensions);
static struct range BisectIndices(const struct NrNrHood* hood, float minVal, float maxVal, size_t axis);


#include <sys/time.h>
#include <time.h>
static float GetTime(void)
{
	static struct timeval timeval, first_timeval;
	gettimeofday(&timeval, 0);
	if(first_timeval.tv_sec == 0) {
		first_timeval = timeval;
		return 0;
	}
	unsigned int s = (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000.0;
    return s / 1000.0f;
}



struct NrNrHood* NrNrCreateNeighborhood(size_t numPoints, size_t dimensions, const float* points) {
    // All data is allocated in a single block. Allocate and then
    // aim some pointers into the space.

    size_t allocSize = sizeof(struct NrNrHood) + sizeof(size_t) * numPoints * 6;
    struct NrNrHood* hood = malloc(allocSize);
    hood->numPoints = numPoints;
    hood->dimensions = dimensions;
    hood->points = points;
    hood->indices = (size_t*)(hood + 1);
    hood->nextIndices = hood->indices + numPoints * 3;

    // Initialize indices with list of all points

    for (size_t i = 0; i < numPoints; ++i) {
        hood->indices[i] = i;
    }
    for (size_t d = 1; d < dimensions; ++d) {
        size_t* initSrc = hood-> indices + (d - 1) * numPoints;
        size_t* initDst = hood-> indices + d * numPoints;
        memcpy(initDst, initSrc, sizeof(size_t) * numPoints);
    }

    // Sort indices for each dimension

    for (size_t d = 0; d < dimensions; ++d) {
        size_t* array = hood->indices + d * numPoints;
        SortIndices(array, numPoints, points, d, dimensions);
    }

    // Find inverse mappifng to sort index (point to index)
    // This goes into a temporary set of inverse indices

    size_t* inverse = malloc(sizeof(size_t) * numPoints);

    // Build next arrays that map a point from one dimension to the next
    // Allows jumping from dim0 index N to dim1, and then to dim2, ...

    for (size_t d = 0; d < dimensions; ++d) {
        // Build an index of the 'inverse' mapping from point number to sorted index
        size_t* ind = hood->indices + numPoints * d;
        for (size_t i = 0; i < numPoints; ++i) {
            inverse[ind[i]] = i;
        }

        // Fill out the next array, move 'ind' to the array coming from
        size_t* next = hood->nextIndices + numPoints * d;
        if (d) {
            ind = hood->indices + numPoints * (d - 1);
        } else {
            ind = hood->indices + numPoints * (dimensions - 1);
        }
        for (size_t i = 0; i < numPoints; ++i) {
            next[i] = inverse[ind[i]];
        }
    }

    free(inverse);

    return hood;
}


struct NrNrHood* NrNrFreeNeighborhood(struct NrNrHood* hood) {
    if (!hood) {
        free(hood);
    }
    return NULL;
}



//size_t NrNrCountBlock(struct NrNrHood* hood, size_t most, float* minimum, float* maximum) {
//}


size_t NrNrSearchBlock(const struct NrNrHood* hood, size_t** results, const float* minimum, const float* maximum) {
    // Find ranges of indices for each dimension
    struct range* ranges = alloca(sizeof(ranges) * hood->dimensions);

    for (size_t d = 0; d < hood->dimensions; ++d) {
        size_t* array = hood->indices + hood->numPoints * d;
        ranges[d] = BisectIndices(hood, minimum[d], maximum[d], d);
    }

    // Find dimension with the least number of matches

    size_t minAxis = 0;
    size_t minSize = ranges[0].end - ranges[0].start;
    for (size_t d = 0; d < hood->dimensions; ++d) {
        size_t size = ranges[d].end - ranges[d].start;
        if (size < minSize) {
            minSize = size;
            minAxis = d;
        }
    }
    if (minSize == 0) {
        *results = NULL;
        return 0;
    }

    // Find values common to each dimension

    size_t* common = (size_t*)malloc(sizeof(size_t) * (minSize + hood->dimensions) + sizeof(size_t*) * hood->dimensions);
    size_t* found = common;

    // To eliminate as much code as possible fro the search loop, precompute several needed values

    size_t* workAxis = common + minSize;
    size_t** workNexts = (size_t**)workAxis + hood->dimensions;
    workNexts[0] = 0; // unused
    workAxis[0] = minAxis; // unused
    for (size_t d = 1; d < hood->dimensions; ++d) {
        workAxis[d] = minAxis + d;
        if (workAxis[d] >= hood->dimensions) {
            workAxis[d] -= hood->dimensions;
        }
        workNexts[d] = hood->nextIndices + hood->numPoints * workAxis[d];
    }

    // The main search loop. All searching time is spent here.
    size_t* minIndices = hood->indices + hood->numPoints * minAxis;

    float start = GetTime();
    for (size_t minStep = ranges[minAxis].start; minStep < ranges[minAxis].end; ++minStep) {
        size_t workStep = minStep;
        int match = 1;
        for (size_t d = 1; d < hood->dimensions; ++d) {
            size_t pt = workNexts[d][workStep];
            if (!(match = (pt >= ranges[workAxis[d]].start && pt <= ranges[workAxis[d]].end)))
                break;
            workStep = pt;
        }
        if (match) {
            *found++ = minIndices[minStep];
        }
    }
    printf("  SCAN: %.3f sec\n", GetTime() - start);

    // Handle return values

    size_t totalFound = found - common;
    if (!totalFound) {
        free(common);
        *results = NULL;
        return 0;
    }
    if (totalFound < minSize) {
        common = realloc(common, sizeof(size_t) * totalFound);
    }

    // realloc common to the smaller, found size
    *results = common;
    return totalFound;
}



size_t* NrNrFreeResults(size_t* results) {
    if (!results) {
        free(results);
    }
    return NULL;
}



void SortIndices(size_t* array, int numPoints, const float* points, size_t axis, size_t dimensions) {
    if (numPoints < 2) {
        return;
    }

    // Insertion sort for small ranges
    if (numPoints < 16) {
        for (size_t i = 1; i < numPoints; ++i) {
            size_t point = array[i];
            float position = points[point * dimensions + axis];
            size_t j;
            for (j = i; j > 0 && points[array[j - 1] * dimensions + axis] > position; --j) {
                array[j] = array[j - 1];
            }
            array[j] = point;
        }
        return;
    }

    // Quicksort (recursive) for most spans
    float position = points[array[numPoints / 2] * dimensions + axis];
    size_t *left = array;
    size_t *right = array + numPoints - 1;
    while (left <= right) {
        while (points[(*left) * dimensions + axis] < position)
            left++;
        while (points[(*right) * dimensions + axis] > position)
            right--;
        if (left <= right) {
            size_t t = *left;
            *left++ = *right;
            *right-- = t;
        }
    }

    SortIndices(array, right - array + 1, points, axis, dimensions);
    SortIndices(left, array + numPoints - left, points, axis, dimensions);
}



struct range BisectIndices(const struct NrNrHood* hood, float minVal, float maxVal, size_t axis) {
    struct range range;
    size_t *fullStart, *fullEnd;
    size_t *start;
    size_t *end;

    fullStart = start = hood->indices + hood->numPoints * axis;
    fullEnd = end = fullStart + hood->numPoints;

    float fullEndValue = hood->points[*fullEnd * hood->dimensions + axis];

    // TODO, early out when maxVal < start point or minVal > end point
    while (start < end) {
        size_t dist = end - start;
        size_t* half = start + (dist / 2);
        float halfValue = hood->points[*half * hood->dimensions + axis];
        if (minVal <= halfValue) {
            end = half;
#if 0
            if (fullEndValue < halfValue) {
                fullEnd = half;
                fullEndValue = hood->points[*fullEnd * hood->dimensions + axis];
            }
#endif
        } else {
            start = half + 1;
        }
    }
    range.start = start - fullStart;

    start = end;
    end = fullEnd;
    while (start < end) {
        size_t dist = end - start;
        size_t* half = start + (dist / 2);
        float halfValue = hood->points[*half * hood->dimensions + axis];
        if (maxVal >= halfValue) {
            start = half + 1;
        } else {
            end = half;
        }
    }
    range.end = start - fullStart;

    return range;
}