How to work around memory errors?

Hi,

I am trying to create an object represented by a matrix of coordinate points. This object will be a subset of a larger surface, say blue. What I have is two matrices of coordinate points for two surfaces, blue and red, one inside the other (red is inside blue), and I am trying to create a third surface that is points of blue within a radius of points of the boundary points of red. The code I have so far is:

idx = true(1,size(blue,2));
for pointb=blue
    pb = repmat(pointb,1,length(boundred));
    idx = idx & sum((pb-boundred).^2)<r^2;
end
A = blue(idx);

where A is the matrix of coordinate points from blue.

This gives me the error:

Error using repmat
Out of memory. Type HELP MEMORY for your options.

which I don't know how to get around. My matrices are very large matrices, one ~6000x3 and the other ~850,000x3. Is there a way to smooth the matrices without losing too much information? Or another method to deal with large data and memory restrictions? Thanks in advance for any help!

Answer 1

[James Tursa]

I guess the obvious suggestion is don't use repmat, but rewrite the sum((pb-boundred).^2) calculation as a loop instead. What are the dimensions of blue and boundred?

[Michael J]

My red surface is ~6000x3 and my blue matrix is ~850000x3. I initially wrote a loop to do this but it took too long to compute. The code for the loop was

A = [];
for pointr = boundred'
    for pointb = blue'
        if sqrt((pointb(1)-pointr(1))^2 + (pointb(2)-pointr(2))^2 + (pointb(3)-pointr(3))^2) < r^2
            A = [A pointb];
        end
    end
end

[James Tursa]

Some issues with this loop. First, you are dynamically increasing the size of A within the loop, which can slow the whole thing down considerably if there are a lot of points to put in A. Second, points in blue could appear multiple times in A if they are near multiple points in boundred, which I don't think is your intent. Third, some calculations that are invariant to the loop index (e.g., transpose blue' and r^2) are done multiple times in a loop when it could have been done once outside the loop. Fourth, you have sqrt(sum-or-squares) < r^2 ... I think you intended to drop the sqrt because you are comparing the calculation to r^2 and not r. In any event, how fast does this run:

r2 = r^2;          % Square of radius
m = size(blue,1);  % Number of blue points
s = false(m,1);    % Pre-allocate logical flag array
for k=1:m          % For every point in blue
    d = bsxfun(@minus,boundred,blue(k,:)); % Difference between blue point and all red points
    s(k) = any(sum(d.*d,2) < r2); % Set flag if any of the difference magnitudes is < radius
end
A = blue(s,:); % Construct answer with logical indexing

If this still isn't fast enough, it might be sped up by converting it to a mex routine where the calculations can easily be short-circuited on the first match. Pre-sorting the points in some fashion might also help this approach.

[Michael J]

Thank you for explaining how my code can be improved on and helping me with a solution! It worked in about 3 minutes! Much faster than the several hours my original code was taking.

[Stephen23]

@Michael J: the improvements that James Tursa made are well documented, you should read and learn these for your future code:

http://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

http://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html