In a vector, how to remove neighbours too close from one another
13 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Guillaume FOURNIER
el 12 de Feb. de 2021
Comentada: Guillaume FOURNIER
el 13 de Feb. de 2021
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
0 comentarios
Respuesta aceptada
Jan
el 12 de Feb. de 2021
Editada: Jan
el 13 de Feb. de 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
Más respuestas (2)
Bruno Luong
el 12 de Feb. de 2021
Editada: Bruno Luong
el 12 de Feb. de 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 comentarios
Bruno Luong
el 13 de Feb. de 2021
Yes, I run the benchmark also on my side thus my comment. Wonder why uniquetol is that inefficient.
Sharareh J
el 12 de Feb. de 2021
Editada: Sharareh J
el 12 de Feb. de 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
Ver también
Categorías
Más información sobre Matrix Indexing en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!