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

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

Answer 1

[Jan]

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

[Guillaume FOURNIER]

Thanks a lot for your detailed answer !

Answer 2

[Bruno Luong]

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.

[Jan]

For the test data of length 1e7 uiniquetol needs 1.28 seconds, while the Mex version needs 0.048 s and the Matlab codes 0.10 s.

[Bruno Luong]

Yes, I run the benchmark also on my side thus my comment. Wonder why uniquetol is that inefficient.

[Guillaume FOURNIER]

Thanks a lot for your answer ! I did not know that function.

Answer 3

[Sharareh J]

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

[Guillaume FOURNIER]

Thanks for your answer