optimizing: fastest way to detect if ANY points are within distance N of points?

16 次查看（过去 30 天）
Carson Purnell
Carson Purnell 2023-2-12
评论： Huiyan Liang 2023-12-9
I'm creating a volume from a very large number of points, created by accumulating a large number of point cloud objects in a box without overlapping any existing points. Everything works so far. But testing if the new cloud is touching anything already accumulated is very slow - each cloud has >1e5 points and the whole thing can accumulate a vast number of points, creating tons of data to sift through.
I'm using a KDTreesearcher object to search the existing points with rangesearch(), after pre-pruning to points within the maximal distance that a new cloud could have points within. It's dramatically faster than pdist2 or looping through the points but still the vast majority of runtime. Making a new KDT when needed is the slower part, though faster than the exhaustive search method.
Is there something faster? Some hidden tool on the FEX that works like a short-circuit rangesearch rather than exhaustively searching each point? I can't find anything even close other than older implementations of KDT from before it was implemented in matlab.
  4 个评论
显示 2更早的评论
James Tursa
James Tursa 2023-5-16
编辑：James Tursa 2023-5-16
Seems like a fairly trivial mex routine for short-circuiting. How are the data represented? Double arrays, or ...? Is there any ordering (e.g., sorting) to the data?
Carson Purnell
Carson Purnell 2023-5-16
Data are nx3 double arrays, no ordering of any kind.
Here's my current testing subfunction. c is the large initial set of points, pts is the small set being tested.
function err = proxtest(c,pts,tol)
l = min(pts,[],1)-tol; h = max(pts,[],1)+tol; %low and high bounds per dimension
ix = c>l & c<h;
ix = find(sum(ix,2)>2); %prefiltering bottleneck
err=0;
if ~isempty(ix) %skip KDT when no overlap is possible
buck = round( size(c,1)/1650 );
modeltree = KDTreeSearcher(c(ix,:),'Bucketsize',buck);
[~,d] = rangesearch(modeltree,pts,tol,'SortIndices',0); %
d = [d{:}]; if any(d<tol), err=1; end %check output distances
end
end

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Kartik
Kartik 2023-5-15
Hi,
One potential approach to make the search faster is to use a bounding box to filter out points that are definitely outside the search radius before performing a more precise search with rangesearch. This can reduce the number of points that are included in the search, potentially improving run time. Another approach is to divide the existing points into smaller sub-sections and perform rangesearch for each sub-section. This can help reduce the search time while still being more precise than using pdist2.
MathWorks provides useful functions for point cloud processing and spatial indexing, which can be useful for performing efficient operations on large point sets. Specifically, the `pointCloud` and `pointCloudPartition` functions can be used for operations on point clouds. The `KDTreeSearcher` and `rangesearch` functions can be used for efficient searching of points. Additionally, the `BoundingBox` function can be used to define a bounding box region of interest.
Here are some documentation links that are relevant to your task:
  1. `KDTreeSearcher`: https://www.mathworks.com/help/stats/kdtreesearcher.html?s_tid=doc_ta
  2. `rangesearch`: https://www.mathworks.com/help/stats/exhaustivesearcher.rangesearch.html?searchHighlight=rangesearch&s_tid=srchtitle_rangesearch_2
  3. `BoundingBox`: https://www.mathworks.com/help/matlab/ref/polyshape.boundingbox.html?searchHighlight=BoundingBox&s_tid=srchtitle_BoundingBox_1
  4. `pointCloud`: https://www.mathworks.com/help/vision/ref/pointcloud.html
  2 个评论
Carson Purnell
Carson Purnell 2023-5-15
I'm already pruning the points outside the search area before using rangesearch. I need something better - short circuit logic, a non-KDT option that can mutate points quickly or search faster, or to MEX my own prefilter function.
So unless matlab has a cumulative octree that I missed it appears that I need some wacky FEX assistance.
Huiyan Liang
Huiyan Liang 2023-12-9
Hello! Could you please clarify the mean of filtering out points and how to realize it specifically? Does it mean to filter out points outside the search radius for each query point?

请先登录，再进行评论。

James Tursa
James Tursa 2023-5-25
编辑：James Tursa 2023-5-25
Here is a naive brute force mex routine you can try out. It makes no attempts at cache optimization or multi-threading, and is probably sensitive to the relative sizes of the matrices involved. But it does short circuit. Maybe for your particular use case it will help with timing. Or if your data were pre-sorted in some way then maybe we could take advantage of that with a different approach.
// proxtest_mex.c tests if any 3D points of two matrices are withing tol of each other
//
// FLAG = proxtest_mex(A,B,tol)
//
// A = Nx3 full real double matrix
// B = Mx3 full real double matrix
// tol = positive tolerance scalar
// FLAG = 0 (nothing close)
// 1 (at least one point in B close to a point in A)
/* Includes ----------------------------------------------------------- */
#include "mex.h"
/* Gateway ------------------------------------------------------------ */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize N, M, i, j;
double tol2;
double *Ax, *Ay, *Az, *Bx, *By, *Bz, *BX, *BY, *BZ;
// Argument checks
if( nlhs > 1 ) {
mexErrMsgTxt("Too many outputs.");
}
if( nrhs != 3 ) {
mexErrMsgTxt("Need exactly three inputs.");
}
if( !mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1]) ||
mxIsSparse(prhs[0]) || mxIsSparse(prhs[1]) ||
mxIsComplex(prhs[0]) || mxIsComplex(prhs[1]) ||
mxGetNumberOfDimensions(prhs[0]) != 2 ||
mxGetNumberOfDimensions(prhs[1]) != 2 ||
mxGetN(prhs[0]) != 3 || mxGetN(prhs[1]) != 3 ) {
mexErrMsgTxt("First two inputs must be Nx3 and Mx3 full real double 2D matrices.");
}
if( !mxIsNumeric(prhs[2]) || mxGetNumberOfElements(prhs[2]) != 1 ) {
mexErrMsgTxt("Third input must be numeric scalar.");
}
// Get sizes and data pointers
N = mxGetM(prhs[0]);
M = mxGetM(prhs[1]);
Ax = mxGetData(prhs[0]); Ay = Ax + N; Az = Ay + N;
BX = mxGetData(prhs[1]); BY = BX + M; BZ = BY + M;
tol2 = mxGetScalar(prhs[2]);
if( tol2 <= 0.0 ) {
mexErrMsgTxt("Third input tolerance must be positive.");
}
tol2 = tol2 * tol2;
// Brute force checking
for( i=0; i<N; i++ ) {
Bx = BX; By = BY; Bz = BZ;
for( j=0; j<M; j++ ) {
if( (*Ax - *Bx)*(*Ax - *Bx) + (*Ay - *By)*(*Ay - *By) + (*Az - *Bz)*(*Az - *Bz) < tol2 ) {
plhs[0] = mxCreateDoubleScalar(1.0);
return; // short-circuit return
}
Bx++; By++; Bz++;
}
Ax++; Ay++; Az++;
}
plhs[0] = mxCreateDoubleScalar(0.0);
}

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by