Search code examples
matlabindexingmatchbsxfun

bsxfun vs repmat using different sized matrix inputs to match indices


I know there are already a lot of questions regarding how bsxfun is supposed to be faster than repmat, so I hope this question isn't too redundant.

I am using a short code with repmat. The purpose of this code is to identify indices between A and B where the row vectors are the same so I can later remove them, but using repmat takes far too long using Matlab R2016a.

[~,i] = min(abs(dot(A',repmat(B,length(A),1)')./sqrt(sum(A'.^2))./norm(B)-1));

Here A and B are not necessarily the same size, so using repmat doesn't seem to work in general anyway. Currently the size of A and B are 8020x3 and 21615x3. I have found that this method of making the sizes of A and B the same works but I am open to other methods as well.

if isequal(size(A),size(B))==1

else
    a = size(A,1);
    b = size(B,1);

    if a<b
        A = [A;nan(abs(B-A),3)];
    else
        B = [B;nan(abs(B-A),3)];
    end
end

I was looking at the bsxfun documentation and on this Mathworks site https://www.mathworks.com/matlabcentral/answers/297088-speed-up-indexing-repmat-operation and it seems like I should be able to replace repmat(B,length(A),1)' with bsxfun(@times, B, length(A)) but they do not output the same final array, so I must be doing something incorrect.

If possible I would like to modify my code above so that it takes less time to calculate, and also have A and B be different sizes such as given above if possible. If necessary, I can still work with A and B being the same size but they will be both be larger, in this case 571000x3 so the calculations will take longer for this reason as well.

Addition: In the comments, it was discussed whether intersect would work. Theoretically, it should since A (black) and B are defined as being inside or outside the isosurface (red) within a larger volume T using the function inpolyhedron https://www.mathworks.com/matlabcentral/fileexchange/37856-inpolyhedron-are-points-inside-a-triangulated-volume-#comments. For some reason, the output of inpolyhedron is a set of the points within the isosurface with some scattered points outside as well. When using inpolyhedron to find the values outside the isosurface, these "extra" scattered points are not included in the set B, so intersect finds an empty set when inputting A and B. Points inside the isosurface


Solution

  • So it turns out that the creator of inpolyhedron has a "patch" that mostly corrects this issue.

    Essentially you need to move the points by a small amount.

    Let's say all of my points are in an array called TestPoints.

    TestPoints = combvec(x',y',z')';  % x, y, and z are column vectors
    NudgedPoints = TestPoints + repmat([1e-10 1e-10 0],length(TestPoints),1);
    

    Create a surface with faces and vertices from the data

    [F,V] = isosurface(X,Y,Z,A,-3);   % you can use any value, and X, Y, and Z are positions from a meshgrid
    in = inpolyhedron(F,V,NudgedPoints,'FlipNormals',false);
    TestPoints_in = TestPoints(in,:);
    

    All credit to this solution goes to Sven the creator of inpolyhedron. You can get in touch with him here.