Home > Net >  How to fast retrieve all spheres containing a certain spatial point
How to fast retrieve all spheres containing a certain spatial point

Time:10-13

Imagine that I have an array of balls in the N-dimensional space. Each ball has a certain radius and location. Just to give a simplified example of circles in a 2D plain as follows:

  1. Center at (0, 0), radius: 5.5
  2. Center at (3.5, 4), radius: 12.1
  3. Center at (10.5, -3.2) radius: 3.1 ...

Hope it explains. This list can be very large. Every ball (in 2D, it's a circle) has a certain radius and location. Balls can overlap (well, geometric objects, not physical ones).

My challenge is, given a point P in the space (say, (12, 4) in the 2D example), I want to quickly search for all balls in the given list (potentially very large list) which contain the point P (distance between P and the center of a ball is shorter than the radius of that ball).

Since the list of balls is long, I can't afford scanning the list for each search. I need some index on this list. However, the location of point P is arbitrary. I wonder whether there's any indexing methods which can support this search?

Thanks very much!

CodePudding user response:

If you can allow for some constraints on the max radius, you might be able to use a modified Vantage Point Tree.

Normally a VP Tree partitions nth-dimensional points and you can query for the closest point to your search point in Log N time.

You'd have to modify it to store the radius and then make the query recursively keep finding the nth nearest neighbor to your search point and then verify the search point is within the radius until your "max radius" distance threshold is exceeded.

Here's a standard VP Tree implementation in JS: https://codepen.io/louisricci/pen/JYaXKb?editors=001


There's other metric tree algorithms like KD Tree and Ball Tree that do a similar job. https://towardsdatascience.com/tree-algorithms-explained-ball-tree-algorithm-vs-kd-tree-vs-brute-force-9746debcd940


Source of the VP Tree code pen linked above, might be a good starting point for you.

"use strict";
const VPLearn = (function() {

const VP_LEAF_SIZE = 20;
const VP_CONSTANT = 1.7;
const Z_SCORE_95 = 1.96;
const CONFIDENCE_INTERVAL_10 = 0.10;
const SAMPLE_CONSTANT = ((Z_SCORE_95 * Z_SCORE_95) * 0.5 * (1 - 0.5)) / (CONFIDENCE_INTERVAL_10 * CONFIDENCE_INTERVAL_10);

function find(arr, d) {
    var l = 0, h = arr.length - 1, i; 
    while(l <= h) { 
        i = (l   h) >> 1; 
        if(arr[i].d < d) {
            l =   i; 
        } else if(arr[i].d > d) {
            h = --i; 
        } else {
            return i;
        } 
    } 
    return -(  l); 
}

class VPValue {
    /**
    @param {!number} index
    @param {!Array.<number>} value
    @param {number=} distance
     */
    constructor(index, value, distance) {
        this.i = index;
        this.v = value;
        this.d = distance||Infinity;
    }
}
class VPNode {
    /**
    @param {number=} index
    @param {number=} radius
    @param {VPNode=} near
    @param {VPNode=} far
    @param {Array.<VPValue>=} values
     */
    constructor(index, radius, near, far, values) {
        this.i = index||0;
        this.r = radius||0;
        this.n = near||null;
        this.f = far||null;
        this.v = values||null;
    }
}
class CustomHeap {
    /**
    @param {!number} size
     */
    constructor(size) {
        this.size = size;
        this.data = [];
    }
    peek() {
        return this.data[0];
    }
    /**
    @param {!VPValue} el
     */
    add(el) {
        if(this.size === this.data.length && this.peek().d <= el.d) {
            return;
        }
        let pos = this.data.length,
            parent;
        this.data.push(el);
        while(pos > 0) { // bubble up
            parent = (pos - 1) >> 1;
            if(this.data[parent].d < this.data[pos].d) {
                const tmp = this.data[parent];
                this.data[parent] = this.data[pos];
                this.data[pos] = tmp;
                pos = parent;
            } else {
                break;
            }
        }
        if(this.size < this.data.length) {
          this.remove();
        }
    }
    /**
    @return {?VPValue}
     */
    remove() {
        if(this.data.length === 0) {
            return null;
        }
        const result = this.data[0],
            head = this.data.pop();
        if(this.data.length > 0) {
            this.data[0] = head;
            let pos = 0, 
                last = this.data.length - 1,
                l, r, m;
            do { // bubble down
                l = (pos << 1)   1;
                r = l   1;
                m = pos;
                if(l <= last && this.data[m].d < this.data[l].d) {
                    m = l;              
                }
                if(r <= last && this.data[m].d < this.data[r].d) {
                    m = r;              
                }
                if(m !== pos) {
                    const tmp = this.data[m];
                    this.data[m] = this.data[pos];
                    this.data[pos] = tmp;
                    pos = m;
                } else {
                    break;
                }
            } while(true);
        }
        return result;
    }
}
class CustomHeap2 {
    /**
    @param {!number} size
     */
    constructor(size) {
        this.size = size;
        this.data = [];
    }
    peek() {
        return this.data[this.data.length - 1];
    }
    /**
    @param {!VPValue} el
     */
    add(el) {
        if(this.size === this.data.length && this.peek().d <= el.d) {
            return;
        }
        const index = find(this.data, el.d);
        this.data.splice((index < 0 ? ~index : index), 0, el);
        if(this.size < this.data.length) {
            this.remove();
        }
    }
    /**
    @return {?VPValue}
     */
    remove() {
        const result = this.data.pop();
        return result;
    }
}

/**
@param {!Array.<number>} v
@return {!number}
 */
function magnitude(v) {
    return Math.sqrt(v.reduce((p, c) => p   c * c, 0));
}
/**
@param {!Array.<number>} v
@return {!Arra.<number>}
 */
function scale(v, s) {
    return v.map(vi => s * vi);
}
/**
@param {!Array.<number>} v
@return {!Arra.<number>}
 */
function weight(v, w) {
    return v.map((vi, i) => w[i] * vi);
}
/**
@param {!Array.<number>} v
@return {!Arra.<number>}
 */
function normalize(v) {
    const dist = magnitude(v);
    return v.map(vi => vi / dist); 
}
/**
@param {!Array.<number>} v1
@param {!Array.<number>} v2
@return {!number}
 */
function distanceSquared(v1, v2) {
    let dist = 0;
    for(let i=0; i<v1.length; i  ) {
        const diff = v1[i] - v2[i];
        dist  = diff * diff;
    }
    return dist;
}
/**
@param {!Array.<number>} v1
@param {!Array.<number>} v2
@return {!number}
 */
function distance(v1, v2) {
    return Math.sqrt(distanceSquared(v1, v2));
}
/**
@param {!Array.<number>} v
@return {!number}
 */
function mean(v) {
    return v.reduce((p, c) => p   c) / v.length;
}
/**
@param {!Array.<number>} v
@return {!number}
 */
function sigmaSquared(v) {
    const mu = mean(v);
    return v.reduce((p, c) => p   Math.pow(c - mu, 2), 0);
}
/**
@param {!Array.<number>} v
@return {!number}
 */
function standardDeviation(v) {
    const mu = mean(v);
    return Math.sqrt(v.reduce((p, c) => p   Math.pow(c - mu, 2), 0) / v.length);
}
/**
@param {!number} populationCount
@return {!number}
 */
function calcSampleSize(populationCount) {
    return (SAMPLE_CONSTANT / (1   ((SAMPLE_CONSTANT - 1) / populationCount)))|0;
}
/**
@param {!number} populationCount
@param {!number} size
@return {!Array.<number>}
 */
function sample(populationCount, size) {
    const set = new Set();
    for(let i=0, pc=populationCount; i<size; i  ) {
        set.add(Math.floor(Math.random() * pc)|0);
    }
    return Array.from(set);
}
/**
@param {!Array.<Array.<number>>} vList
@return {!Array.<number>}
 */
function centeroid(vList) {
    const sum = Array(vList[0].length).fill(0);
    vList.forEach(vi => vi.forEach((vij, j) => sum[j]  = vij));
    return scale(sum, 1 / vList.length);    
}

/**
@param {!Array.<Array.<number>>} vList
@return {!Array.<VPValue>}
 */
function pickVantagePoints(vList) {
    const sampleSize = Math.min(vList.length, calcSampleSize(vList.length)),
        count = Math.max(1, Math.log(vList.length) / Math.log(VP_CONSTANT))|0,
        vpSet = new Set(),      
        vpList = [];
    let bestIndex, bestValue;
    // pick root vantage point, candidate with largest std dev from samples
    const candidates = sample(vList.length, sampleSize);
    bestIndex = -1;
    bestValue = 0;
    for(let i=0; i<candidates.length; i  ) {
        const index = candidates[i],
            candidate = vList[index],
            stdDev = sigmaSquared(sample(vList.length, sampleSize).map(vi => distance(vList[vi], candidate)));
        if(stdDev >= bestValue) {
            bestIndex = index;
            bestValue = stdDev;
        }
    }
    vpSet.add(bestIndex);
    vpList.push(new VPValue(bestIndex, vList[bestIndex]));
    // pick 2nd vantage point, furthest sample from 1st vantage point
    if(count > 1) {
        bestIndex = -1;
        bestValue = 0;
        sample(vList.length, sampleSize).forEach(vi => {
            const dist = distanceSquared(vpList[0].v, vList[vi]);
            if(dist >= bestValue) {
                bestIndex = vi;
                bestValue = dist;
            }       
        });
        vpSet.add(bestIndex);
        vpList.push(new VPValue(bestIndex, vList[bestIndex]));
    }
    // pick remaining vantage points, furthest from other vantage points
    while(vpSet.size < count) {
        bestIndex = -1;
        bestValue = 0;
        sample(vList.length, sampleSize).forEach(vi => {
            const meanDistSq = mean(vpList.map(vpi => distanceSquared(vpi.v, vList[vi])));
            if(meanDistSq >= bestValue) {
                bestIndex = vi;
                bestValue = meanDistSq;
            }       
        });
        if(!vpSet.has(bestIndex)) {
            vpSet.add(bestIndex);
            vpList.push(new VPValue(bestIndex, vList[bestIndex]));
        }       
    }
    //
    return vpList;
}

/**
@param {!Array.<VPValue>} vpList
@param {!Array.<Array.<number>>} vList
@return {!VPNode}
 */
function buildVantagePointTree(vpList, vList) {
    const vnList = vList.map((vi, i) => new VPValue(i, vi)),
        root = new VPNode(),
        queue = [root, vnList];
    let qIndex = 0, node, list, bestIndex, bestValue, pivot, radius, near, far;
    // partition data
    while(qIndex < queue.length) { 
        node = queue[qIndex  ];
        list = queue[qIndex  ];
        // leaf node
        if(list.length <= VP_LEAF_SIZE) {
            node.v = list;
            continue;
        }
        // root node
        if(node === root) {
            bestIndex = 0;
        } else { // choose best vp
            bestIndex = -1;
            bestValue = -1;
            vpList.forEach((vpi, i) => {
                const vpiValue = vpi.v,
                    vpIndex = i,
                    samples = sample(list.length, calcSampleSize(list.length)),
                    tmp = samples.map(idx => distance(list[idx].v, vpiValue)),
                    stdDev = sigmaSquared(tmp);
                if(stdDev > bestValue) {
                    bestIndex = vpIndex;
                    bestValue = stdDev;
                }
            });
        }
        // calc distances from vp to all elements
        list.forEach(vni => vni.d = distance(vni.v, vpList[bestIndex].v));
        // sort highest to lowest distance      
        list.sort((a, b) => b.d - a.d);
        // pivot
        pivot = Math.max(0, (list.length >> 1) - 1);
        radius = list[pivot].d;
        while(pivot   1 < list.length && radius === list[pivot   1].d) {
            pivot  ;
        }
        // near and far partitions
        // far are all GREATER THAN OR EQUAL TO RADIUS
        far = list.slice(0, pivot   1);     
        // near are all LESS THAN RADIUS
        near = list.slice(pivot   1);
        // node
        node.i = bestIndex;
        node.r = radius;
        if(far.length > 0) {
            node.f = new VPNode();
            queue.push(node.f, far);
        }
        if(near.length > 0) {
            node.n = new VPNode();
            queue.push(node.n, near);
        }
    }
    return root;
}

class VPLearn {
    /**
    @param {boolean=} normalizeInputs Normalize all input vectors
    @param {Array.<number>=} customWeights Multiply components of input vectors by corresponding weight 
     */
    constructor(normalizeInputs, customWeights) {
        this._inputs = [];
        this._targets = [];
        this._vantagePoints = [];
        this._tree = null;
        this._weights = customWeights || null;
        this._normalize = normalizeInputs || false;
    }
    /**
    Add an input-target pair.
    @param {Array.<number>} input Training vector
    @param {*} target Training target
     */
    train(input, target) {      
        if(this._weights !== null) {
            input = weight(input, this._weights);
        }
        if(this._normalize === true) { 
            input = normalize(input);
        }
        this._inputs.push(input);
        this._targets.push(target);
        this._tree = null;
    }
    /**
    Add multiple training objects.
    @param {{input:Array.<number>, target:*}} trainingData
     */
    trainRange(trainingData) {
        trainingData.forEach(data => this.train(data.input, data.target));
    }
    /**
    Compile training data into a vantage point tree.
     */
    compile() {
        this._vantagePoints = pickVantagePoints(this._inputs);
        this._tree = buildVantagePointTree(this._vantagePoints, this._inputs);
    }
    /**
    Find the nearest input values and return the corresponding targets.
    @param {Array.<number> input Search vector
    @param {number=} nearest Nth nearest results
    @param {Array.<{i:number, v:Array.<number>, d:number}>} inputResults OUT array of nth nearest tree leafs: (i)ndex, (v)alue, and (d)istance from input
    @returm {Array.<*>} Nearest target values
     */
    query(input, nearest = 1, inputResults = null) {
        if(this._weights !== null) {
            input = weight(input, this._weights);
        }
        if(this._normalize === true) { 
            input = normalize(input);
        }
        if(this._tree === null) {
            this.compile(); 
        }       
        const vpList = this._vantagePoints.map(vpi => new VPValue(vpi.i, vpi.v)),
            stack = [this._tree],
            nth = new CustomHeap(nearest),
            result = [];
        let tau = Number.MAX_VALUE,
            node, vp, dist, tmp;
        while(stack.length > 0) {
            node = stack.pop();
            // value node
            if(node.v !== null) {
                node.v.forEach(vi => {
                    vi.d = distance(input, vi.v);
                    nth.add(vi);
                    if(nth.data.length === nearest) {
                        tau = nth.peek().d;
                    }
                });
                continue;
            }
            // vantage point
            vp = vpList[node.i];
            if(vp.d !== Infinity) {
                dist = vp.d;
            } else {
                vp.d = dist = distance(input, vp.v);
            }
            if(dist < node.r) {
                if(dist < node.r   tau && node.n !== null) {
                    stack.push(node.n);
                }
                if(dist >= node.r - tau && node.f !== null) {
                    stack.push(node.f);
                }
            } else {
                if(dist >= node.r - tau && node.f !== null) {
                    stack.push(node.f);
                }
                if(dist < node.r   tau && node.n !== null) {
                    stack.push(node.n);
                }
            }
        }
        // return the target values in closest to furtherest order of the corresponding input       
        if(Array.isArray(inputResults)) {
            while((tmp = nth.remove()) !== null) {
                result.push(JSON.parse(JSON.stringify(this._targets[tmp.i])));
                inputResults.push(JSON.parse(JSON.stringify(tmp)));
            }
            inputResults.reverse();
            return result.reverse();
        } else {
            while((tmp = nth.remove()) !== null) {
                result.push(JSON.parse(JSON.stringify(this._targets[tmp.i])));
            }
            return result.reverse();
        }
    }
}

return VPLearn;
})();

document.body.innerHTML  = '<pre style="overflow:scroll; height:10em;"></pre>';
const outp = document.querySelector("pre");

const compass = new VPLearn();
compass.trainRange([
    {input: [0], target:"N"},
    {input: [22.5], target:"NNE"},
    {input: [45], target:"NE"},
    {input: [77.5], target:"ENE"},
    {input: [90], target:"E"},
    {input: [112.5], target:"ESE"},
    {input: [135], target:"SE"},
    {input: [157.5], target:"SSE"},
    {input: [180], target:"S"},
    {input: [202.5], target:"SSW"},
    {input: [225], target:"SW"},
    {input: [247.5], target:"WSW"},
    {input: [270], target:"W"},
    {input: [292.5], target:"WNW"},
    {input: [315], target:"NW"},
    {input: [337.5], target:"WNW"},
    {input: [360], target:"N"}
    ]);
compass.compile();
outp.innerHTML  = "Compass (Degrees -> Direction)\n";
for(let i=0; i<=360; i  ) {
    outp.innerHTML  = `${i} -> ${compass.query([i])[0]}\n`;
}

const colors = new VPLearn();
colors.trainRange([
{input: [240,248,255], target: "aliceblue"},
{input: [250,235,215], target: "antiquewhite"},
{input: [0,255,255], target: "aqua"},
{input: [127,255,212], target: "aquamarine"},
{input: [240,255,255], target: "azure"},
{input: [245,245,220], target: "beige"},
{input: [255,228,196], target: "bisque"},
{input: [0,0,0], target: "black"},
{input: [255,235,205], target: "blanchedalmond"},
{input: [0,0,255], target: "blue"},
{input: [138,43,226], target: "blueviolet"},
{input: [165,42,42], target: "brown"},
{input: [222,184,135], target: "burlywood"},
{input: [95,158,160], target: "cadetblue"},
{input: [127,255,0], target: "chartreuse"},
{input: [210,105,30], target: "chocolate"},
{input: [255,127,80], target: "coral"},
{input: [100,149,237], target: "cornflowerblue"},
{input: [255,248,220], target: "cornsilk"},
{input: [220,20,60], target: "crimson"},
{input: [0,255,255], target: "cyan"},
{input: [0,0,139], target: "darkblue"},
{input: [0,139,139], target: "darkcyan"},
{input: [184,134,11], target: "darkgoldenrod"},
{input: [169,169,169], target: "darkgray"},
{input: [0,100,0], target: "darkgreen"},
{input: [189,183,107], target: "darkkhaki"},
{input: [139,0,139], target: "darkmagenta"},
{input: [85,107,47], target: "darkolivegreen"},
{input: [255,140,0], target: "darkorange"},
{input: [153,50,204], target: "darkorchid"},
{input: [139,0,0], target: "darkred"},
{input: [233,150,122], target: "darksalmon"},
{input: [143,188,143], target: "darkseagreen"},
{input: [72,61,139], target: "darkslateblue"},
{input: [47,79,79], target: "darkslategray"},
{input: [0,206,209], target: "darkturquoise"},
{input: [148,0,211], target: "darkviolet"},
{input: [255,20,147], target: "deeppink"},
{input: [0,191,255], target: "deepskyblue"},
{input: [105,105,105], target: "dimgray"},
{input: [30,144,255], target: "dodgerblue"},
{input: [178,34,34], target: "firebrick"},
{input: [255,250,240], target: "floralwhite"},
{input: [34,139,34], target: "forestgreen"},
{input: [255,0,255], target: "fuchsia"},
{input: [220,220,220], target: "gainsboro"},
{input: [248,248,255], target: "ghostwhite"},
{input: [255,215,0], target: "gold"},
{input: [218,165,32], target: "goldenrod"},
{input: [128,128,128], target: "gray"},
{input: [0,128,0], target: "green"},
{input: [173,255,47], target: "greenyellow"},
{input: [240,255,240], target: "honeydew"},
{input: [255,105,180], target: "hotpink"},
{input: [205,92,92], target: "indianred"},
{input: [75,0,130], target: "indigo"},
{input: [255,255,240], target: "ivory"},
{input: [240,230,140], target: "khaki"},
{input: [230,230,250], target: "lavender"},
{input: [255,240,245], target: "lavenderblush"},
{input: [124,252,0], target: "lawngreen"},
{input: [255,250,205], target: "lemonchiffon"},
{input: [173,216,230], target: "lightblue"},
{input: [240,128,128], target: "lightcoral"},
{input: [224,255,255], target: "lightcyan"},
{input: [250,250,210], target: "lightgoldenrodyellow"},
{input: [144,238,144], target: "lightgreen"},
{input: [211,211,211], target: "lightgrey"},
{input: [255,182,193], target: "lightpink"},
{input: [255,160,122], target: "lightsalmon"},
{input: [32,178,170], target: "lightseagreen"},
{input: [135,206,250], target: "lightskyblue"},
{input: [119,136,153], target: "lightslategray"},
{input: [176,196,222], target: "lightsteelblue"},
{input: [255,255,224], target: "lightyellow"},
{input: [0,255,0], target: "lime"},
{input: [50,205,50], target: "limegreen"},
{input: [250,240,230], target: "linen"},
{input: [255,0,255], target: "magenta"},
{input: [128,0,0], target: "maroon"},
{input: [102,205,170], target: "mediumaquamarine"},
{input: [0,0,205], target: "mediumblue"},
{input: [186,85,211], target: "mediumorchid"},
{input: [147,112,219], target: "mediumpurple"},
{input: [60,179,113], target: "mediumseagreen"},
{input: [123,104,238], target: "mediumslateblue"},
{input: [0,250,154], target: "mediumspringgreen"},
{input: [72,209,204], target: "mediumturquoise"},
{input: [199,21,133], target: "mediumvioletred"},
{input: [25,25,112], target: "midnightblue"},
{input: [245,255,250], target: "mintcream"},
{input: [255,228,225], target: "mistyrose"},
{input: [255,228,181], target: "moccasin"},
{input: [255,222,173], target: "navajowhite"},
{input: [0,0,128], target: "navy"},
{input: [253,245,230], target: "oldlace"},
{input: [128,128,0], target: "olive"},
{input: [107,142,35], target: "olivedrab"},
{input: [255,165,0], target: "orange"},
{input: [255,69,0], target: "orangered"},
{input: [218,112,214], target: "orchid"},
{input: [238,232,170], target: "palegoldenrod"},
{input: [152,251,152], target: "palegreen"},
{input: [175,238,238], target: "paleturquoise"},
{input: [219,112,147], target: "palevioletred"},
{input: [255,239,213], target: "papayawhip"},
{input: [255,218,185], target: "peachpuff"},
{input: [205,133,63], target: "peru"},
{input: [255,192,203], target: "pink"},
{input: [221,160,221], target: "plum"},
{input: [176,224,230], target: "powderblue"},
{input: [128,0,128], target: "purple"},
{input: [255,0,0], target: "red"},
{input: [188,143,143], target: "rosybrown"},
{input: [65,105,225], target: "royalblue"},
{input: [139,69,19], target: "saddlebrown"},
{input: [250,128,114], target: "salmon"},
{input: [244,164,96], target: "sandybrown"},
{input: [46,139,87], target: "seagreen"},
{input: [255,245,238], target: "seashell"},
{input: [160,82,45], target: "sienna"},
{input: [192,192,192], target: "silver"},
{input: [135,206,235], target: "skyblue"},
{input: [106,90,205], target: "slateblue"},
{input: [112,128,144], target: "slategray"},
{input: [255,250,250], target: "snow"},
{input: [0,255,127], target: "springgreen"},
{input: [70,130,180], target: "steelblue"},
{input: [210,180,140], target: "tan"},
{input: [0,128,128], target: "teal"},
{input: [216,191,216], target: "thistle"},
{input: [255,99,71], target: "tomato"},
{input: [64,224,208], target: "turquoise"},
{input: [238,130,238], target: "violet"},
{input: [245,222,179], target: "wheat"},
{input: [255,255,255], target: "white"},
{input: [245,245,245], target: "whitesmoke"},
{input: [255,255,0], target: "yellow"},
{input: [154,205,50], target: "yellowgreen"} 
]);
colors.compile();
console.log(colors._tree);
console.log(colors._vantagePoints);
document.body.innerHTML  = '<input id="cvR" type="number" value="255"/><input id="cvG" type="number" value="255"/><input id="cvB" type="number" value="255"/><input type="button" id="btnCV" value="Top 2 Colors"/><div id="CV"></div>';
const MAX_COLOR_DIST = Math.sqrt(255*255   255*255   255*255);
document.querySelector("#btnCV").addEventListener("click", e => {
    const rgb = [
        parseInt(document.querySelector("#cvR").value),
        parseInt(document.querySelector("#cvG").value),
        parseInt(document.querySelector("#cvB").value)
    ];
    const meta = [];
    const result = colors.query(rgb, 2, meta);
    document.querySelector("#CV").innerHTML = `<div style="border:4px solid rgb(${rgb});" >rgb(${rgb}) : "${result.toString()} ${(MAX_COLOR_DIST - meta[0].d) * 100 / MAX_COLOR_DIST}%, ${(MAX_COLOR_DIST - meta[1].d) * 100 / MAX_COLOR_DIST}%"</div>`;
}, false);
  • Related