In a vector, how to remove neighbours too close from one another

55 views (last 30 days)
Hi everyone,
I searched for quite a while to avoid asking a question that was already answered. I found some resembling questions but not exactly what i want.
I will try to be as clear as possible. Here is what i want to do. I have a vector x. For each element in this vector, i want to analyze the distance with the next element, and remove the next one if too close (below a certain threshold). But i want to do that avoiding a for-loop if possible.
Let's have a quick example. Given the x vector below, i want to remove the points if the distance is below the threshold 10.
x = [1 6 12 17 23 25 34]
If i use the function diff, this is not entirely satisfying.
diff(x)<10
will return me a logical array, where it is true all the time. Indeed, the distance between 1 and 6 is below 10, same between 6 and 12, and so on. So the command
x(diff(x)<10)
will give me an empty vector. Whereas, what i would like to have is the output
x = [1 12 23 34]
because once i remove a point, the next one is different and the threshold might be satisfied.
I could do that with a for loop, however i would prefer to avoid it, as i know that it is not the best practice for efficient computation, and use built-in Matlab functions.
I tried to come up with some ideas combining some function such as diff, cumsum, movsum to do it, but nothing satisfying.
Maybe someone has some good ideas? I guess i am not the first one trying to solve that problem. But maybe the for loop is inevitable after all.
Thanks in advance for your inputs.
All the best.
Guillaume

Accepted Answer

Jan
Jan on 12 Feb 2021
Edited: Jan on 13 Feb 2021
I assume a C-mex function is the best solution:
// File: MovThresh.c
// Compile: mex -O MoveThresh.c
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
size_t n;
double *x, *xEnd, *y, *y0, Thresh, C;
if (nrhs != 2) {
mexErrMsgIdAndTxt("JS:MovThresh:BadNInput", "2 inputs needed.");
}
if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || !mxIsNumeric(prhs[1])) {
mexErrMsgIdAndTxt("JS:MovThresh:BadInput",
"1st input must be a real double, 2nd a numerical value.");
}
Thresh = mxGetScalar(prhs[1]);
n = mxGetNumberOfElements(prhs[0]);
if (n == 0) { // Early return for empty input:
plhs[0] = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxREAL);
return;
}
// Create output:
plhs[0] = mxCreateUninitNumericMatrix(1, n, mxDOUBLE_CLASS, mxREAL);
// Get pointers to data:
#if MX_HAS_INTERLEAVED_COMPLEX
x = mxGetDoubles(prhs[0]);
y = mxGetDoubles(plhs[0]);
#else
x = mxGetPr(prhs[0]);
y = mxGetPr(plhs[0]);
#endif
// Remember pointer to output and limit for input:
y0 = y;
xEnd = x + n;
// Loop over data for moving thresholding:
*y++ = *x;
C = *x + Thresh;
while (++x < xEnd) {
if (*x > C) {
*y++ = *x;
C = *x + Thresh;
}
}
// Crop unused elements of the output:
mxSetN(plhs[0], y - y0);
return;
}
Methods with loops in Matlab:
function Test
X = cumsum(randi([0, 20], 1, 1e7)); % Test data
tic; Y1 = MovThresh(X, 10); toc
tic; Y2 = MovThresh_M1(X, 10); toc
tic; Y3 = MovThresh_M2(X, 10); toc
isequal(Y1, Y2, Y3)
end
% ***************************************
function Y = MovThresh_M1(X, Thresh)
if isempty(X), Y = []; return; end
M = false(size(X));
M(1) = true;
C = X(1);
for iX = 2:numel(X)
if X(iX) - C > Thresh
M(iX) = true;
C = X(iX);
end
end
Y = X(M);
end
% ***************************************
function Y = MovThresh_M2(X, Thresh)
if isempty(X), Y = []; return; end
Y = zeros(size(X));
iY = 1;
Y(iY) = X(1);
for iX = 2:numel(X)
if X(iX) - Y(iY) > Thresh
iY = iY + 1;
Y(iY) = X(iX);
end
end
Y = Y(1:iY);
end
Timings:
Elapsed time is 0.042547 seconds. C-Mex
Elapsed time is 0.133265 seconds. Logical indexing
Elapsed time is 0.101082 seconds. Collect double and crop
Note: The logical indexing is not implemented efficiently in Matlab. Use FEX: CopyMask (MEX) for a speed up:
% Replace: Y = X(M); by
Y = CopyMask(X, M).';
Elapsed time is 0.099813 seconds. Logical indexing with CopyMask

More Answers (2)

Bruno Luong
Bruno Luong on 12 Feb 2021
Edited: Bruno Luong on 12 Feb 2021
Using stock function
>> uniquetol([1 6 12 17 23 25 34],10,'DataScale',1)
ans =
1 12 23 34
Jan's various codes are faster. Make one wonder what uniquetol does internally to waste that much speed.
  3 Comments
Bruno Luong
Bruno Luong on 13 Feb 2021
Yes, I run the benchmark also on my side thus my comment. Wonder why uniquetol is that inefficient.

Sign in to comment.


Sharareh J
Sharareh J on 12 Feb 2021
Edited: Sharareh J on 12 Feb 2021
You can try this one:
x = [1 6 12 17 23 25 34];
i = length(x);
j = 0;
while i ~= j
j = i;
x(find(diff(x)<10, 1) + 1) = [];
i = length(x);
end

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!