1. Peter Shinners
  2. nrnr

Source

nrnr / nrnr.c

#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;
};



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



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 * dimensions);

    for (size_t d = 0; d < dimensions; ++d) {
        size_t* ind = hood->indices + numPoints * d;
        size_t* inv = inverse + numPoints * d;
        for (size_t i = 0; i < numPoints; ++i) {
            inv[ind[i]] = i;
        }
    }

    // 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) {
        size_t* next = hood->nextIndices + numPoints * d;
        size_t* inv = inverse + numPoints * d;
        size_t* ind;
        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] = inv[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);
    size_t* found = common;

    for (size_t minStep = ranges[minAxis].start; minStep < ranges[minAxis].end; ++minStep) {
        size_t d;
        size_t workStep = minStep;
        for (d = 1; d < hood->dimensions; ++d) {
            size_t workAxis = minAxis + d;
            if (workAxis >= hood->dimensions) {
                workAxis -= hood->dimensions;
            }
            size_t pt = hood->nextIndices[hood->numPoints * workAxis + workStep];
            if (pt < ranges[workAxis].start || pt > ranges[workAxis].end) {
                break;
            }
            workStep = pt;
        }
        if (d == hood->dimensions) {
            *found = hood->indices[hood->numPoints * minAxis + minStep];
            ++found;
        }
    }

    // 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;

    // 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;
        } 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;
}