Application Of Mask in piv_FFTmulti

[ABHIMANYU DUBEY]

Dear William, 

I've noticed that in the function pic_FFTmulti(), the mask is being applied at the end of the computation of the correlation matrix. Once the correlation matrix has been evaluated, the following snippet of code removes the vectors at the points which lie inside the mask.

%apply mask

       masked_xs = (minix:step:maxix) + round(interrogationarea/2);

       masked_ys = (miniy:step:maxiy) + round(interrogationarea/2);

       typevector(mask(masked_ys, masked_xs)) = 0;

       result_conv(:, :, mask(masked_ys, masked_xs)) = 0;

       if multipass == passes

           correlation_map(mask(masked_ys, masked_xs)) = 0;

       end

What could happen because of this is the vectors which lie closer to the mask will not be accurate. Kindly correct me if I'm wrong here.

We could improve this slightly by eliminating the pixels inside the mask each time before we compute the correlation matrix. Maybe something like, 

    N = size(img1, 3);

   n = size(img1, 1) * size(img1, 2);

   a = reshape(img1, [n N]);

   b = reshape(img2, [n N]);

   a( mask(:) ) = [] ; 

   b( mask(:) ) = [] ; 

   mean_a = sum(a) / n;

   mean_b = sum(b) / n;

   corr_map = zeros(N, 1);

   for i=1:N

       % All this is a long, but fast way of calculating

       % corr_map(i) = corr2(img1(:,:,i), img2(:,:,i))

       a_ = a(:,i) - mean_a(i);

       b_ = b(:,i) - mean_b(i);

       corr_map(i) = sum(a_.*b_) / sqrt(sum(a_.*a_) * sum(b_.*b_));

   end

The above snippet of code is to get my idea across. Of course, the exact implementation would be slightly different,

If this interests you, I would be happy to discuss the implementation details and create a pull request on Git Hub after incorporating those changes in the code. 

Thanks,

Regards,

Abhimanyu

[thielickeoptolution]

Hi Abhimanyu,

I tried something similar before, and that part of the code is still available in piv_fft_multi.m (as a comment), somewhere near line 141:

%{
        %Improve masking?
        max_img_value=(max(image1_roi(:))+max(image2_roi(:)))/2;
        noise_mask1=rand(size(image1_roi))*max_img_value*0;
        noise_mask2=rand(size(image1_roi))*max_img_value*0;
        image1_roi(mask==1)=0;
        image2_roi(mask==1)=0;
        noise_mask1(mask==0)=0;
        noise_mask2(mask==0)=0;
        image1_roi=image1_roi+noise_mask1;
        image2_roi=image2_roi+noise_mask2;
        %keyboard
        disp('XXX')
        %}

Here, I replaced the parts that are masked with random noise to get zero correlation. I think I read somewhere that this would improve the output..

In this comment, I multiplied the generated noise with zero, just to check whether noise on or off makes any difference. I don't remember what I did exactly, and I think I didn't have a test case where this made a difference. It would be nice if you test and improve this part!

[ABHIMANYU DUBEY]

Dear William, 

Thanks for your input. 

Yeah, I had observed the part of the code you mentioned here but couldn't determine why it was commented.  Replacing the parts that are masked with random noise to get zero correlation would be an excellent way to approach this. I'll work on this and report any changes. 

Also, is there any specific reason you have introduced masks as a cell array containing the boundary pixels of the masked area? I'm trying to get PIVlab to work with dynamic masks. So passing only the boundary pixels and then reconstructing the mask like this by looping over each individual connected component is too much work. Of course, if we have a constant mask for all the image pairs, then this would be an excellent implementation that greatly reduces the amount of memory used. I assume you had this in mind while developing the code. :-)

            A=im2bw(A,0.5); %#ok<*IM2BW>

           A1=zeros(size(A));

           A2=A1;A3=A1;A4=A1;

           %cut mask in 4 pieces to minimize parent / child / hole problems in masks

           rowshalf=round(size(A,1)/2);

           colshalf=round(size(A,2)/2);

           %A1(1:rowshalf,1:colshalf) = A(1:rowshalf,1:colshalf);%top left part

           %A2(1:rowshalf,colshalf+1:end) = A(1:rowshalf,colshalf+1:end);%top right half

           %A3(rowshalf+1:end,1:colshalf) = A(rowshalf+1:end,1:colshalf); % lower left part

           %A4(rowshalf+1:end,colshalf+1:end) = A(rowshalf+1:end,colshalf+1:end); %lower right part

           A1(1:rowshalf,1:colshalf) = A(1:rowshalf,1:colshalf);%top left part

           A2(1:rowshalf,colshalf:end) = A(1:rowshalf,colshalf:end);%top right half

           A3(rowshalf:end,1:colshalf) = A(rowshalf:end,1:colshalf); % lower left part

           A4(rowshalf:end,colshalf:end) = A(rowshalf:end,colshalf:end); %lower right part

           %A(:,round(size(A,2)/2))=0;

           %A(round(size(A,1)/2),:)=0;

           %B=A;

           %B=im2bw(abs(A-B));

           bwbound=[bwboundaries(A1); ...

bwboundaries(A2) ; ...

bwboundaries(A3) ; ...

bwboundaries(A4)];

           %bwbound=bwboundaries(A);

           importmaskx=cell(size(bwbound,1),1);

           importmasky=importmaskx;

           for i=1:size(bwbound,1)

               temp=bwbound{i,1};

               importmasky{i,1}=temp(1:3:end,1); % ABHI: changed 1:10:end to 1:3:end

               temp=bwbound{i,1};

               importmaskx{i,1}=temp(1:3:end,2); % ABHI: changed 1:10:end to 1:3:end

           end

           maskiererx=retr('maskiererx');

           maskierery=retr('maskierery');

I'll modify the code to read/input mask in piv_analysis() as a logical array and then send it to pin_FFTmulti() as an additional input. I'll also modify piv_FFTmulti() to work with this. Kindly let me know if you have any better ideas to approach this.

So I'll first work on implementing this with the command line. Once I'm done with this, I may need your help to incorporate this with "PIVlab_gui.m". This m-file is quite a zig-saw puzzle :-). Since this is one of many projects I'm currently working on, it may take me some time to finish. I'll get back to you once I've finished the command-line part, and we'll take it up from there. 

 

Thanks once again for your support. 

With warm regards,

Abhimanyu 

[Zuhaib Nissar]

Hi William,

Thank you for developing PIVlab and supporting the research community. Your contributions are much appreciated.

I am currently using PIVlab for time-resolved PIV analysis of a batch of 1000 images. I would like to apply a logical mask to the entire batch and have found some very useful conversations on this topic. I have modified PIVlab_commandline.m and piv_FFTmulti() to pass a logical image matrix (Binary_mask) as an argument in mask settings instead of concatenating xy indices (see below).

s{4,1}= 'Mask'; s{4,2}=[Binary_mask]; % If needed, generate via: imagesc(image); [temp,Mask{1,1},Mask{1,2}]=roipoly;

I am not sure if I have followed it correctly. I would expect 0 values in typevector_filt, which is not the case. Could you please let me know where I went wrong.

Best wishes,

Zuhaib

[Zuhaib Nissar]

Hello Both,

I have found solution to the problem asked earlier. I have revised the function piv_analysis() to apply mask on the image pair before piv_FFTmulti() is called.

Many thanks.

Zuhaib

[ABHIMANYU DUBEY]

Hello Zuhaib, 

I'm glad you were able to resolve it. I didn't respond to your query because I didn't understand what you were trying to ask in the first place. 

Regarding my case, I was trying to adapt the PIVlab code to incorporate dynamic masking in the images. I deal with turbulent particle-gas suspensions in a channel, and hence, to the best of my knowledge, I can only use PIVlab indirectly as it stands currently. Simply because, in my case, every pair of images will have a unique mask since the position of particles changes from image to image

Also, I observed that in the function pic_FFTmulti(), the mask is being applied at the end of the computation of the correlation matrix. Once the correlation matrix has been evaluated, the code removes the vectors at the points inside the mask. To me, that didn't felt right. I think these pixels should have been removed from the matrix before they were used to compute the correlation matrix. 

I've incorporated those changes locally in my code. But I still need to share it here. After incorporating these changes, when I superimposed my PIV profiles on the DNS profiles obtained for the same Re, St, and volume fraction, they were not overlapping within the error bar. I wanted to know if my code was working correctly, whether I needed to improve it further, or if there was some issue with the DNS code itself. So my colleague (whose DNS profiles I had used) and I are trying to resolve this. I'm quite confident with the DNS code because it's been developed in the lab for 20 years by multiple PhD students, and they have had multiple publications with it in JFM, with the most recent one in 2020. But regardless, we'll need to resolve this issue.  

That's why I have yet to reply to this thread. 

I hope that helps, 

With warm regards,

Abhimanyu

[William Thielicke (Shrediquette)]

Dear  Abhimanyu,

in my test cases, the masked areas mostly were shadows that do not have any correlation anyway. This is why I couldn't see any significant difference when masking the pixel image before tha analysis. If you have test cases where this makes a difference, then this would be interesting to see and to implement. Maybe open an issue on github?

[William Thielicke (Shrediquette)]

Just to summarize for all people that measure close to boundaries:

For everyone willing to contribute to the best method to measure close to masked boundaries, please test this small modification code and report your opinion:

in PIV_FFT_multi.m, there is a commented out section that potentially improves the measurements close to masked regions. Search for the string %Improve masking? in that file and replace the commented out code with the following (latest) version:

%Improve masking?

max_img_value=(max(image1_roi(:))+max(image2_roi(:)))/2;
random_number_state=rng;rng('default');rng(1);
noise_mask1=rand(size(image1_roi))*max_img_value*1;
noise_mask2=rand(size(image1_roi))*max_img_value*1;
rng(random_number_state);

image1_roi(mask==1)=0;
image2_roi(mask==1)=0;
noise_mask1(mask==0)=0;
noise_mask2(mask==0)=0;
image1_roi=image1_roi+noise_mask1;
image2_roi=image2_roi+noise_mask2;