nrnr / nrnr.c

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



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



static void SortIndices(unsigned short* array, unsigned short numPoints, const float* points, size_t axis, size_t dimensions);
static struct range BisectIndices(const unsigned short* array, const float* points, unsigned short numPoints, size_t dimensions, 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(unsigned short) * numPoints * 6;
    struct NrNrHood* hood = malloc(allocSize);
    hood->numPoints = numPoints;
    hood->dimensions = dimensions;
    hood->points = points;
    hood->chunks = (numPoints + USHRT_MAX - 1) / USHRT_MAX;
    hood->indices = (unsigned short*)(hood + 1);
    hood->nextIndices = hood->indices + numPoints * 3;

    // Allocate temporary buffer for storing the inverse sorted array
    // Find inverse mappifng to sort index (point to index)
    // This goes into a temporary set of inverse indices
    size_t inverseSize = USHRT_MAX;
    if (inverseSize > numPoints) {
        inverseSize = numPoints;
    }
    unsigned short* inverse = (unsigned short*)malloc(sizeof(unsigned short) * inverseSize);

    // Break the data into small chunks
    unsigned short* chunkIndices = hood->indices;
    unsigned short* chunkNextIndices = hood->nextIndices;
    const float* chunkPoints = points;
    size_t indexStep = hood->dimensions * USHRT_MAX;
    for (size_t c = 0; c < hood->chunks; ++c,
                chunkIndices += indexStep, chunkNextIndices += indexStep, chunkPoints += indexStep) {
        size_t chunkSize = numPoints - (c * USHRT_MAX);
        if (chunkSize > USHRT_MAX) {
            chunkSize = USHRT_MAX;
        }

        for (unsigned short i = 0; i < chunkSize; ++i) {
            chunkIndices[i] = i;
        }

        // Copy ordered point numbers from dimension 0 into all others
        for (size_t d = 1; d < dimensions; ++d) {
            unsigned short* initSrc = chunkIndices + (d - 1) * chunkSize;
            unsigned short* initDst = chunkIndices + d * chunkSize;
            memcpy(initDst, initSrc, sizeof(unsigned short) * chunkSize);
        }

        // Sort indices for each dimension
        for (size_t d = 0; d < dimensions; ++d) {
            unsigned short* array = chunkIndices + d * chunkSize;
            SortIndices(array, chunkSize, chunkPoints, d, dimensions);
        }

        // 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
            unsigned short* ind = chunkIndices + chunkSize * d;

            for (unsigned short i = 0; i < chunkSize; ++i) {
                inverse[ind[i]] = i;
            }

            // Fill out the next array, move 'ind' to the array coming from
            unsigned short* next = chunkNextIndices + chunkSize * d;
            if (d) {
                ind = chunkIndices + chunkSize * (d - 1);
            } else {
                ind = chunkIndices + chunkSize * (dimensions - 1);
            }
            for (unsigned short i = 0; i < chunkSize; ++i) {
                next[i] = inverse[ind[i]];
            }
        }
    }

    free(inverse);
    return hood;
}


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



size_t NrNrSearchBlock(const struct NrNrHood* hood, size_t** results, const float* minimum, const float* maximum) {
    struct range* ranges = alloca(sizeof(ranges) * hood->dimensions);

    size_t* workAxis = alloca(sizeof(size_t) * hood->dimensions);;
    unsigned short** workNexts = alloca(sizeof(unsigned short*) * hood->dimensions);
    workAxis[0] = hood->dimensions; // unused (will crash if used)
    workNexts[0] = NULL; // unused (will crash if used)

    size_t* common = (size_t*)malloc(sizeof(size_t) * hood->numPoints);
    size_t* found = common;

    // Break the data into small chunks
    unsigned short* chunkIndices = hood->indices;
    unsigned short* chunkNextIndices = hood->nextIndices;
    const float* chunkPoints = hood->points;
    size_t chunkOffset = 0;
    size_t totalEmptyChunks = 0;
    size_t indexStep = hood->dimensions * USHRT_MAX;
    for (size_t c = 0; c < hood->chunks; ++c,
                chunkIndices += indexStep, chunkNextIndices += indexStep, chunkPoints += indexStep, chunkOffset += USHRT_MAX) {
        size_t chunkSize = hood->numPoints - (c * USHRT_MAX);
        size_t chunkFound = 0;
        if (chunkSize > USHRT_MAX) {
            chunkSize = USHRT_MAX;
        }

        // Find ranges of indices for each dimension
        int emptyBisect = 0;
        for (size_t d = 0; d < hood->dimensions; ++d) {
            unsigned short* array = chunkIndices + chunkSize * d;
            ranges[d] = BisectIndices(array, chunkPoints, chunkSize, hood->dimensions, minimum[d], maximum[d], d);
            if (ranges[d].start == ranges[d].end) {
                emptyBisect = 1;
                totalEmptyChunks++;
                break;
            }
        }
        if (emptyBisect) {
            continue;
        }

        // Find dimension with the least number of matches
        size_t minAxis = 0;
        unsigned short minSize = ranges[0].end - ranges[0].start;
        for (size_t d = 0; d < hood->dimensions; ++d) {
            unsigned short size = ranges[d].end - ranges[d].start;
            if (size < minSize) {
                minSize = size;
                minAxis = d;
            }
        }
        // To eliminate as much code as possible for the search loop,
        // precompute several needed values

        unsigned short* minIndices = chunkIndices + chunkSize * minAxis;
        for (size_t d = 1; d < hood->dimensions; ++d) {
            workAxis[d] = minAxis + d;
            if (workAxis[d] >= hood->dimensions) {
                workAxis[d] -= hood->dimensions;
            }
            workNexts[d] = chunkNextIndices + chunkSize * workAxis[d];
        }

        // The main search loop. All searching time is spent here.

        for (unsigned short minStep = ranges[minAxis].start; minStep < ranges[minAxis].end; ++minStep) {
            unsigned short workStep = minStep;
            int match = 1;
            for (size_t d = 1; d < hood->dimensions; ++d) {
                unsigned short pt = workNexts[d][workStep];
                if (!(match = (pt >= ranges[workAxis[d]].start && pt < ranges[workAxis[d]].end)))
                    break;
                workStep = pt;
            }
            if (match) {
                *found++ = minIndices[minStep] + chunkOffset;
                ++chunkFound;
            }
        }
    }

    // Handle return values
    size_t totalFound = found - common;
    if (!totalFound) {
        free(common);
        *results = NULL;
        return 0;
    }

    // Realloc common to the smaller, found size
    if (totalFound < hood->numPoints) {
        common = realloc(common, sizeof(size_t) * totalFound);
    }

    *results = common;
    return totalFound;
}



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



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

    // Insertion sort for small ranges
    if (numPoints < 16) {
        for (unsigned short i = 1; i < numPoints; ++i) {
            unsigned short point = array[i];
            float position = points[point * dimensions + axis];
            unsigned short 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];
    unsigned short *left = array;
    unsigned short *right = array + numPoints - 1;
    while (left <= right) {
        while (points[(*left) * dimensions + axis] < position)
            left++;
        while (points[(*right) * dimensions + axis] > position)
            right--;
        if (left <= right) {
            unsigned short 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 unsigned short* array, const float* points, unsigned short numPoints, size_t dimensions, float minVal, float maxVal, size_t axis) {
    struct range range;
    const unsigned short *fullStart, *fullEnd;
    const unsigned short *start;
    const unsigned short *end;

    fullStart = start = array;
    fullEnd = end = fullStart + numPoints;

    float startValue = points[array[0] * dimensions + axis];
    if (startValue > maxVal) {
        range.start = range.end = numPoints;
        return range;
    }

    float finalValue = points[array[numPoints - 1] * dimensions + axis];
    if (finalValue < minVal) {
        range.start = range.end = 0;
        return range;
    }

    while (start < end) {
        unsigned short dist = end - start;
        const unsigned short* half = start + (dist / 2);
        float halfValue = points[*half * dimensions + axis];
        if (minVal <= halfValue) {
            end = half;
            if (finalValue < halfValue) {
                fullEnd = half;
                finalValue = points[*fullEnd * dimensions + axis];
            }
        } else {
            start = half + 1;
        }
    }
    range.start = start - fullStart;

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

    return range;
}
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.