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

diff()differenceMATLABneighborneighboursthresholdvector

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

Best Answer

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
Related Question