How to simplify this piece of code ( maybe vectorizing?)

Question

AliHg 2020-9-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/593680-how-to-simplify-this-piece-of-code-maybe-vectorizing

评论： Image Analyst 2020-9-16

在 MATLAB Online 中打开

Hi.

I have a binary matrix entitled "Anything" and implemented the following steps on it:

Using bwlabel(Anything,8) function, I evaluated different clusters and put them in "X"
Using repmat function, I duplicated "X", 2*2 times and put it in "Repsol"
Using repmat function, I duplicated "Anything", 2*2 times and rewrited "Anything"
Next, I Preallocated a matrix named "p"
Evaluating the probability of positioning random vectors within "Anything" in a way that both head and tail of each vector be placed in specific clusters, is the ultimate goal of this code

I wrote the following piece of code for these purposes. I have to mention that, the part with the nested loops were my initial code and someone in this forum simplified it a little bit a year ago. I commented this simplification. Although, this simplification made my code faster, I need it to become many times faster. So, my question is regarding simplifying this code.

Any suggestion would be appriciated.

% Anything = a binary 2D matrix
X = bwlabel(Anything,8);
s = size(X);
RepSol = repmat(X,2,2);
Anything = repmat(Anything,2,2);
% Preallocation
p = zeros(s);
% for m = 1:N
%     for n = 1:N
%         cond = RepSol(1:s(1),1:s(2)) == RepSol((1:s(1))+m-1,(1:s(2))+n-1);
%         P = Anything(1:s(1),1:s(2)) .* Anything((1:s(1))+m-1,(1:s(2))+n-1);
%         res = cond .* P ;
%         p(m,n) = sum(res(:));
%     end
% end
for m = 1:s(1)
    for n = 1:s(1)
        for i = 1:s(1)
            for j = 1:s(1)
                if RepSol(i,j) == RepSol(m+i-1,n+j-1)
                    p(m,n) = p(m,n) + Anything(i,j) * Anything(m+i-1,n+j-1);
                end
            end
        end
    end
end
Result = p / numel(X);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Image Analyst 2020-9-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/593680-how-to-simplify-this-piece-of-code-maybe-vectorizing#answer_494548

To determine the likelihood of a head or tail being in a particular blob, you can simply use the area fraction. So to have them both in a blob (head and tail in the same blob) is the area fraction squared.

Now you replicated the binary image and the labeled image so each label now is in 4 blobs instead of just one unique blob. One watchout is what happens at the edge where a blob with one label (on the left edge) might touch a blob with a different label from the right side of the image. Do any blobs touch the side of the image? Anyway...

To learn the probability of the head being in one of the 4 blobs of a particular label, and the tail being in one of the other blobs with the same label would require you to relabel the replicated binary image so you'll know if the head is in blob #1 of (say) blob with label 3, and the tail is in blob #2 of the blob #3. So you'll use both RepSol and the other labeled image of the replicated logical image. Then you can either use the area fraction method I mentioned, or use a Monte Carlo approach like you're sort of suggesting, but didn't actually do. Either should work. I just don't know why you're doing it in the first place. Why? What is the use case?

4 个评论
显示 2更早的评论隐藏 2更早的评论

Image Analyst 2020-9-15

在 MATLAB Online 中打开

Yes, it would. I'm just not sure why the image needs to be replicated 2-by-2. Why?

To do it non-Monte Carlo, imagine a simple case where you had only two blobs. Blob #1 has 10% of the area of the whole image, and blob #2 has 30% of the area. So if you were to randomly pick one endpoint of the line, what's the probability it would be in blob #1? 10%, right? And what's the probability the second point would be in blob #2? 30%, right? So the probability that one endpoint is in blob#1 and one endpoint is in blob #2 is 0.1 * 0.3 = 0.03, right? No need to actually do the random endpoint picking unless you want to.

If you want to do a Monte Carlo method, just have a loop that iterates millions of times getting each endpoint. Then at each endpoint check the label in your labeled image. The label is either the row or column of your 2-D histogram that you'll be building up. Something like (untested):

[labeledImage, numBlobs] = bwlabel(binaryImage);
[rows, columns, numColorChannels] = size(binaryImage);
numExperiments = 1000000;
hist2d = zeros(numBlobs, numBlobs);
for k = 1 : numExperiments
    % Get two random endpoints.
    r = randi(rows, 2, 1);
    c = randi(columns, 2, 1);
    % Find out what blob number they are in.
    label1 = labeledImage(r(1), c(1));
    label2 = labeledImage(r(2), c(2));
    % Increment the histogram for that pairing.
    if label1 > 0 && label2 > 0
        % Both endpoints are inside a blob.
        hist2d(label1, label2) = hist2d(label1, label2) + 1;
    end
end
imshow(hist2d, []); % Display the histogram.  Could use surf() if you want.

You might want to skip the case where label1 == label2 where they are both in the same blob.

if label1 == label2
    continue;
end

AliHg 2020-9-16

编辑：AliHg 2020-9-16

Thank you for your comprehensive suggestion.

I would definitely try to implement the aforementioned code. The reason for implementing "repmat" is to use periodic boundary condition, which is an inherent stipulation while dealing with correlation functions.

I have another question regarding the code I posted above. the nested for loops try to place different-sized vectors into the image and evaluate the probability of positioning them in a specific phase (say phase 1, which represents black areas), provided they are in the same cluster. These vectors will be randomly placed inside the binary image in all directions (x, y, diagonal, etc.). What if I want these vectors to be placed only in orthogonal directions, say "x" and "y" directions only? Which part of my code should be changed? One of the reasons for this alternation is to speed up the algorithm.

Thank you in advance

Image Analyst 2020-9-16

在 MATLAB Online 中打开

If you want blobs to "wrap around" and be considered the same blob, which normally they won't be, then you need to do the replication of the binary image ONLY. If you do the labeled image, you'll have a blob with two different labels in contact if the blob touched the edge, so don't bother to replicate the labeled image - we're going to relabel it anyway.
Next (after replication) you need to call imclearborder() to get rid of blobs touching the border. You don't need them anymore since they're now in the middle of the image because of the replication.
Next relabel the new, larger, border-cleared binary image. That way the blobs on the edges (which are now in the middle line and column of the image) will now have a single label number.

If you want to get random coordinates where you get only coordinates in the same line or column, you can do this:

    r = randi(rows, 2, 1);
    c = randi(columns, 2, 1);
    % Half the time make the row be the same and have the time make the column be the same.
    if rand(1) < 0.5
        % Make both row coordinates the same.
        % Line will be a horizontal line in one row.
        r(2) = r(1);
    else
        % Make both column coordinates the same.
        % Line will be a vertical line in one column.
        c(2) = c(1);
    end