Openpiv Matlab: reshape question

[francesc...@edu.unito.it]

Hello,
I'm working with openpiv with Matlab 2013a: during the analysis of these specific attached images i found a strange behaviour concerning the interpolation of some filtered array in the middle flow.

I thought that the anomaly was in the script "inpaint_nans.m" but I understood, in a second time, that this is probably due to a wrong reshape of the filt_res variables in openpivgui.m before calling the interpolation in fill_holes (lines 565-566).

In details: since the reshape start to fill the columns of my new matrices, I guess that the two lines

u = reshape(res(:,3), numrows,numcols);
v = reshape(res(:,4), numrows,numcols);

create a grid in which each point has not always the right neighbors. This bring to a wrong interpolation because inpaint_nans will not take the right neighbors (I saw this problem because, only in this particular image, there is a strong difference between middle flow and boundary flow: for this reason some middle filtered velocity vector is replaced with an interpolated vector having a lower intensity in comparison with its real neighbors).

So I reshaped the matrices with:

u = reshape(res(:,3), numcols,numrows);
v = reshape(res(:,4), numcols,numrows);

Inverting numcols and numrows I have created the grid with all the right neighbors; now the results of the interpolation are better than the previous ones.

Is this reasonable or maybe there is something that I've not considered?

Thanks in advance,

Francesco Toselli

[Alex Liberzon]

Thanks a lot Francesco. It is indeed a serious bug. Could you fork the master and send us back the pull request with the fixed bug (better with a small test, e.g. an expected result of your two images) or at least raise an issue on openpiv-matlab repository and I'll fix it? 

Alex

------------------------

Indirizzo istituzionale di posta elettronica degli studenti e dei laureati dell'Università degli Studi di Torino

Official University of Turin email address for students and graduates 

--
You received this message because you are subscribed to the Google Groups "openpiv-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to openpiv-user...@googlegroups.com.
Visit this group at https://groups.google.com/group/openpiv-users.
For more options, visit https://groups.google.com/d/optout.

--

Alex Liberzon, Assoc. Prof.

Head of the Undergraduate Program, School of Mechanical Engineering, Tel Aviv University, Tel Aviv 69978, Israel

E-mail: ale...@tau.ac.il, Office: +972-3-640-8928, Lab: +972-3-640-6860 (telefax), www.eng.tau.ac.il/~alexlib

[Alex Liberzon]

Dear Francesco

You're absolutely correct - it's a bug. I took an opportunity to fix it and also improve a bit the visualization of the results such that the replaced arrows are marked red and the result is in green. Please test the present version 

https://github.com/OpenPIV/openpiv-matlab/archive/master.zip

To all users, it's strongly recommended to update your OPENPIV-MATLAB version to the latest trunk

https://github.com/OpenPIV/openpiv-matlab/archive/master.zip

 

[Alex Liberzon]

Hi Francesco, all the fixes and discussions are better to do on the forum - others might have a quicker and better answer than I do. Plus it's then preserved and will be helpful to others as well. So I copy the answer to the forum too (removing your private e-mail from it) 

On Tue, Mar 15, 2016 at 11:51 AM, francesco toselli wrote:

Dear Alex,

I'm sorry for my late reply.

Yes, you were right..I had a spatialbox version different from yours. Now I downloaded the most recent version from your github.

I didn't understand at point 3): I had to announce my findings on the google group? The ones you have already changed or next things I will find?

Yes, it's useful for the record of changes and reasons of changes. Github provides something called "Issues" - where each issue can get it's resolution and therefore all is clear. 

 

I'm writing you a single e-mail only because now I have some technical questions to do:

In matlab openpiv:
1) is 8 pixel the maximal resolution we can reach?

this is what we added to the list of options. I do not believe PIV can do much better, but if you're interested, it's easy to change. 

 

About python,version I'm trying to do the same analysis obtained with matlab to understand if the build of the modules, filters and inpaint_nans work as the matlab version.

Honestly I don't know very much of python environment and for this reason I'm stopping at every little different (I don't want to understand the scripts inside the modules but I would like to understand a few things):

1)What are the principal differences between matlab and python version? I mean: the steps sequence in processing images is the same?
Furthermore, I think I understood that there isn't a GUI or an analysis spatialbox for python yet isn't it?

the basic steps are the same between Matlab and Python versions if the settings are the same. Python version though has also two extensions that are not implemented in Matlab (just because we think it's not useful to develop something for the closed software anymore): Cython which is C written code with much faster execution time, linked to Python. (all the *.pyx files) 

 

2)I'm only at the first tutorial to understand how to process a pair of images but I saw that filter-module is different than matlab local filtering with the edge-kernel.

"u, v, mask = openpiv.validation.sig2noise_val( u, v, sig2noise, threshold = 1.3 )" 

This mask is referring for global vectors and is not a good filter for my pair of images (middle flow is quite bigger compared to the boundary one: with this threshold I will filter all middle vectors and I don't know how to filter locally with a edge-kernel convolution). Is this the only filter typology in python?

obviously the lack of documentation is what's "killing" us, but we have not many hands on that. There's just a start and the tutorials which are useful so far  http://openpiv.readthedocs.org/en/latest/index.html

 http://openpiv.readthedocs.org/en/latest/src/tutorial.html

signal to noise ratio is not a global filter, it's a value that is estimated during the cross-correlation and estimate of the highest peak. In some cases it takes the peak-to-the-mean value and in others peak-to-the-second-peak value. Above some threshold (or below it I think) we think that the correlation is too noisy and therefore the signal is not sufficient. 

 

Because I looked for another filter but I didn't find anything in the pdf guide with all modules...filtering the same pair of images(in attachment) with matlab and python I obtained different results (the python ones are not well filtered)

Yes one needs to select the same kind of parameters. The basic versions give the same results - we tested it in the past. It would be great if you could play with it and find the differences - we'll try to get to the bottom of it. 

 

3)I have still not found a way to process a sequence of images? I found in the tutorial-part2 an example to process in batch a list of image pairs but not a sequence of images like in matlab.

The python version can process the sequence just by using the loop through the names. 

The attached snipet of code can do much better: run the sequence on a multi-processing machine, using parallelism (or job distribution, basically). 

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>
# <codecell>
# Follow the example in tutorial 2: http://openpiv.readthedocs.org/en/latest/src/tutorial.html
# <codecell>
import sys, os, glob
sys.path.append(os.path.abspath(r'D:\OpenPIV\OpenPIV-0.11'))
import openpiv.tools
import openpiv.process
import openpiv.scaling
import openpiv.preprocess
import numpy as np
# import matplotlib.pyplot as plt
from skimage import io, img_as_float, exposure, data, img_as_uint
from skimage.filter import sobel
from skimage.util.dtype import dtype_range
from scipy.ndimage import median_filter, gaussian_filter, binary_fill_holes
from skimage.filter import rank, threshold_otsu
from skimage.morphology import disk
# <codecell>
def dynamic_masking(image,filter_size=7,threshold=0.005):
    """ Dynamically masks out the objects in the PIV images
    """
    imcopy = np.copy(image)
    image = exposure.rescale_intensity(img_as_float(image), in_range=(0, 1))
    blurback = gaussian_filter(image,filter_size)
    edges = sobel(blurback)
    blur_edges = gaussian_filter(edges,21)
    # create the boolean mask 
    bw = (blur_edges > threshold)
    bw = binary_fill_holes(bw)
    imcopy -= blurback
    imcopy[bw] = 0.0
    return imcopy #image
# <codecell>
def func( args ):
    """A function to process each image pair."""
    # this line is REQUIRED for multiprocessing to work
    # always use it in your custom function
    file_a, file_b, counter = args
    
    print file_a, file_b, counter
    filepath, filename = os.path.split(file_a)
    filedrive, filepath = os.path.splitdrive(filepath)
    filename = os.path.splitext(filename)
    filepath = os.path.join('d:',filepath)
    if not os.path.exists(filepath):
        os.makedirs(filepath)
    filename = os.path.join(filepath,filename[0])
    
    print filepath
    # typical parameters:
    window_size = 32 #pixels
    overlap = 16 # pixels
    search_area_size = 64 # pixels 
    frame_rate = 40 # fps
    # read images into numpy arrays
    frame_a  = openpiv.tools.imread( file_a )
    frame_b  = openpiv.tools.imread( file_b )
    
    # %%timeit
    im_a  = io.imread( file_a )
    im_b  = io.imread( file_b )
    
    # let's crop the region of interest
    frame_a =  im_a[600:1600,800:1400]
    frame_b =  im_b[600:1600,800:1400]
    # mask 
    frame_a = dynamic_masking(frame_a)
    frame_b = dynamic_masking(frame_b)
    # process again with the masked images, for comparison# process once with the original images
    u, v, sig2noise = openpiv.process.extended_search_area_piv(
                                                           frame_a.astype(np.int32) , frame_b.astype(np.int32), 
                                                           window_size = window_size,
                                                           overlap = overlap, 
                                                           dt=1./frame_rate, 
                                                           search_area_size = search_area_size, 
                                                           sig2noise_method = 'peak2peak')
    x, y = openpiv.process.get_coordinates( image_size = frame_a.shape, window_size = window_size, overlap = overlap )
    u, v, mask = openpiv.validation.global_val( u, v, (-100.,100.),(-100.,100.))
    u, v, mask = openpiv.validation.sig2noise_val( u, v, sig2noise, threshold = 1.6 )
    # u, v = openpiv.filters.replace_outliers( u, v, method='localmean', max_iter=1, kernel_size=2)
    # x, y, u, v = openpiv.scaling.uniform(x, y, u, v, scaling_factor = 96.52 )
    # save to a file
    openpiv.tools.save(x, y, u, v, mask, filename+'_%03d.txt' % counter, fmt='%8.7f', delimiter='\t')
# <codecell>
# '/Users/alex/Desktop/mark_23_06/5ppm_6hz/Camera1-0100.tif'
listdir = glob.glob(r'Z:\2.5PPM\6_86Hz\*')
for data_dir in listdir:
    task = openpiv.tools.Multiprocesser( data_dir = data_dir, pattern_a=r'Camera1-0*.tif', pattern_b = None)
    task.run( func = func, n_cpus = 1 )

or see here: 

https://bitbucket.org/snippets/alexlib/dK9bk

Have I mentioned that OpenPIV now runs on Google Cloud platform, Amazon Web Services or any other cloud (IBM, Microsoft Azure, etc.) - you can get infinite power for a very short time to finish 10,000 image pairs in 1 hr. 

 

I don't know if I was quite clear

Anyway this is not urgent so answer me only when you can. I don't want to steal time.

Thank you for your help and availability!

Best regards

Francesco T.

2016-03-07 18:56 GMT+01:00 Alex Liberzon <ale...@eng.tau.ac.il>:

see below

On Mon, Mar 7, 2016 at 5:31 PM, francesco toselli <francesc...@edu.unito.it> wrote:

Dear Alex,

About SPATIALTOOLBOX-OPENPIV-MATLAB, I attached 4 files concernig the same processed image: 2 files ( .txt and .vec) with units pix/dt; 2 files ( .txt and .vec) with units m/s.

All files have the first line of description (i.e., VARIABLES= "X pix", "Y pix", "U pix/dt", "V pix/dt", "CHC", ZONE I=33, J=63)

Can you try to load these 4 files with spatialbox gui?

I understood that:

1) Files.txt are not read because they have the first line of description (but I don't know where the script need to be modified for this bug)

TXT file should be depreciated if the new openpiv-matlab version is used. however we have people using old versions  of openpiv-matlab and openpiv-python so they use the TXT format. Maybe we shall add the note somewhere that the format has changed - in my case a dialog pops up and says - try VEC loader

 

2) If I load with spatialbox the file_pixdt.txt and the file_ms.txt (without first line of description) I have 2 different orientation of y axes (beware that only scale and dt factors are changed between the 2 openpiv run). Furthermore, analizing the file.txt in units m/s I can't use any option of the spatialbox GUI. Does it happen to you?

the loader in the spatialbox does some intelligent reading and if the position values are not integers, then it's millimeters (quite stupid I'd say as it could be meters or any other thing, but for some reason we left it as is).otherwise, I created the _old.txt files (attached) and when I read them one is from 100 to 1000 pixels roughly (from top to bottom, bottom is like small streamwise velocity, top is high velocity, e.g. boundary layer proflle) and the second one is from 0.01 mm to 0.08 mm or so, but oriented same way from top to bottom (which is probably incorrect - as smaller "y" or "z" values should be at the bottom) 

 is it what you're refering to? 

 

3) About file.vec, as I said in a previous mail, I can use flood only if I change the line in the script vecread.m ( inverting variables i and j at line 138: data = reshape(data,[j,i,columns]); )

see the new version of master on MY (alexlib) repository. 

https://github.com/alexlib/openpiv-spatial-analysis-toolbox/tree/master

If you could please announce your findings on the mailing list, you can tell them that it's fixed in the trunk on github and everyone is recommended to upgrade their versions of openpiv-matlab and openpiv-spatial-toolbox that now go synchronously. 

Thanks a lot for your careful tests. If you could also write things up - we'll make an official test for openpiv. I also need to do a lot of uploads to the official OpenPIV repository, just my hands are too short. Please, check that if you used the data from the https://github.com/OpenPIV/openpiv-spatial-analysis-toolbox, it's probably different from what I have in mine. Some bugs are laying there for quite some time.  

 

Can you test these 3 points?

Thanks in advance for all your help

Francesco T

2016-03-07 9:32 GMT+01:00 francesco toselli <francesc...@edu.unito.it>:

Ok..I am now able to run the tutorial script.

Since I don't know very well how python works (how it calls scripts etc..) I would like to show only the results image of tutorial: it seems to be slightly different to the image in the pdf tutorial file.

One vector is not filtered isn'it? (could be a problem of my python version or a difference in the pdf tutorial-part1?)

Thank you for your help

F.Toselli

2016-03-02 19:26 GMT+01:00 Alex Liberzon <ale...@eng.tau.ac.il>:

the documentation is also outdated. i opened a pull request to update the INSTALL. thanks for noticing. 

On Wed, Mar 2, 2016 at 3:57 PM, francesco toselli <francesc...@edu.unito.it> wrote:

I uploaded my different file.m in my 2 forks: https://github.com/francescotoselli/openpiv-matlab and https://github.com/francescotoselli/openpiv-spatial-analysis-toolbox

I don't know if you are already able to see them. (the changes in correlation.m is only one try to see results of normalized correlation but I don't know if that is right or not...I didn't pull the request because I don't want to share anything of wrong).

About python: I've read on INSTALL file that for mac was recommended the binary (32 or 64 bit) Enthought Python Distribution (EPD) so I downloaded it for that reason...anyway I can also download Anaconda (version 2.7 or 3.5 of python?)

I will also install Anaconda on my Window 10 if needed.

2016-03-01 22:30 GMT+01:00 Alex Liberzon <ale...@eng.tau.ac.il>:

you go to your github page https://github.com/francescotoselli/openpiv-matlab

you shall see a button "Compare" next to the last commit 

and then you'll get also the option to send Pull request to openpiv/openpiv-matlab repository

here is the tutorial: https://help.github.com/articles/creating-a-pull-request/

at the moment I did it for you, but we can also close the pull request and you do a new one. 

on all machines I strongly recommend to install Anaconda instead of Python(x,y) - it's much more uniform across platforms, much more convenient, etc. https://www.continuum.io/downloads 

it will help us also to identify whether the problems are the installation of the software. 

 

On Tue, Mar 1, 2016 at 5:51 PM, francesco toselli <francesc...@edu.unito.it> wrote:

I forked and created my openpiv-matlab-master but I never used git-hub before so I have to understand how to pull a request or make issue. Can you help me with some instructions?

Meanwhile I'm trying to use also the python version to see the differences between python and Matlab version (even if I never used python before)
On my laptop Windows I downloaded Python(x,y) and some modules: now I can run scripts..for example tutorials seems to run well

On the university machines (Python 2.7.10 on my Mac in local and Python 2.7.3 on the virtual machines Debian) I have also downloaded the python-folder from git-hub and added some missing modules with pip, but I can't run the tutorial scripts.
Especially for the Mac problem it tells me that module builtins.open is missing and I can't install it with pip because it tells me that can't find a version that satisfies the requirement builtins.

It could be a problem concerning python version?

Thank you for your help!

2016-02-29 21:09 GMT+01:00 Alex Liberzon <ale...@eng.tau.ac.il>:

sorry for not replying earlier - quite busy these days. I'll take a look into these issues sooner or later. Thanks a lot. I'd be very happy if you could at least fill the issue on github - it will keep us more focused on fixing bugs like in the spatialbox, for instance. 

regarind inpaint_nan - I admit that I never tried to go deep enough into the different methods. My tests were agains synthetic data and some real experiments and we were happy with the default result. I do not think it makes a lot of difference when we do interpolation of small "holes". When we interpolate large ones, everything is kind of "wrong". 

I think we can fix the filter, just need to dig into the kernels that have left-top-right-bottom settings. it's a modification of a default contrast enhancement one. 

best regards

Alex

ps Nice publication

On Thu, Feb 25, 2016 at 6:30 PM, francesco toselli <francesc...@edu.unito.it> wrote:

The images refer to a second phase to a microburst study.

The first phase of the study (with a different implementation: vertical current induced by "iceberg" fusion water) was published as

Physical simulation of atmospheric microbursts,, jgr,  Volume 119, Issue 11, 16 June 2014 , Pages 6292–6305
http://onlinelibrary.wiley.com/doi/10.1002/2013JD021243/full

This second "experimental package" (performed with vertical currents induced by denser water) has been enough analyzed but the article is still in progress. So there is not a direct reference.

About the new openpiv_matlab_master: now all seems good.

I have only some more questions concerning two things that I've not completely understood:

1) In the filter, using kernel, the boundary convolution is calculated with less point because out of boundaries the kernel takes null values. For this reason I guess that a lot of boundary vectors are filtered (even if they are good).

Do you think there is a way to fix it?

2) About interpolation I can't completely understand how inpaint_nans interpolator works (using method 0 or method 3 I always obtain the same values; moreover I can't understand if all methods are linear or not). I tried to understand the script but even on the web I didn't find a good explanation of this interpolation method.

Can you explain the basic theory behind it or send me some reference?

P.S. Talking about the reshape problem I think that in the Spatial_box package also the vecread.m line 138: data = reshape(data,[i,j,columns]); must be replaced by inverting i and j variables. With the original line I had problem using flood options with the spatialbox interface

Best regards

F.Toselli

2016-02-24 20:51 GMT+01:00 Alex Liberzon <ale...@eng.tau.ac.il>:

can I use your images for an additional test in the openpiv? if yes, what should be the proper reference to those? 

On Wed, Feb 24, 2016 at 7:03 PM, francesco toselli <francesc...@edu.unito.it> wrote:

Sorry,

Here are the two folders

Bests

F.Toselli

2016-02-24 18:01 GMT+01:00 francesco toselli <francesc...@edu.unito.it>:

Dear Alex,

I never used github before so I will try to fork and pull request in the next days.

I start to send two folders of the same couple of images analysis:

A) Difference fill_holes is a folder with 2 processing having a different line in fill_holes script:

number 4 is the original folder having the line: vector = inpaint_nans(real(vector)) + i*inpaint_nans(imag(vector));
number 5 is mine that has this line: vector = inpaint_nans(vector);

Conclusion A: I dont know why, but I found that also this line replaced gave reasonable and better interpolated vectors.

B) Difference analysis is a folder with 2 processing that show the problem that I explained in the post on google group:

Because of conclusion A in both runs I used the fill_hole line: vector = inpaint_nans(vector);

number 2 is my processing (with the new lines of reshape that I will try to send you very soon with github)

numeber 3 is the original one

As you can see, there is a differences between the filtered/interpolated vectors of the two runs due to the different matrix reshape

Best regards

Francesco Toselli

Questa e-mail è stata inviata da un computer privo di virus protetto da Avast. www.avast.com

[francesc...@edu.unito.it]

Dear all,
I'm using the MATLAB 2014a version of OPENPIV by testing and try to understood how to choose the best parameters from the openpivgui.

Sometimes, while computing images by choosing S/N=peak2secondpeak, it happens to have this matlab error:

Operands to the || and && operators must be convertible to logical scalar values.

Error in find_displacement_rect (line 47)
        if x2 > 1 && y2 > 1 && x2 < NfftHeight && y2 < NfftWidth

Error in openpivgui>pushbutton_start_Callback (line 532)
                        [peak1,peak2,pixi,pixj] = find_displacement_rect(c,s2ntype);

Error in gui_mainfcn (line 95)
        feval(varargin{:});

Error in openpivgui (line 42)
    gui_mainfcn(gui_State, varargin{:});

Error in @(hObject,eventdata)openpivgui('pushbutton_start_Callback',hObject,eventdata,guidata(hObject))

 
Error while evaluating uicontrol Callback

Is clearly a problem of second peak computed on boundary of the interrogation window isn't?

How can avoid this problem without changing my resolution parameters? (Because even if I change the size of the interrogation window this problem can appear in other different images of my sequence)

Thanks a lot for your help!

Francesco

[Alex Liberzon]

Can you please add a small example when it happens? 

[Alex Liberzon]

I cannot reproduce this error 

On Wednesday, April 6, 2016 at 6:05:11 PM UTC+3, Alex Liberzon wrote:

Can you please add a small example when it happens? 

[Alex Liberzon]

See below: 

Moreover I would like to know 1 more thing about lsgradient/gradient function in matlab..

1) As it computes the matrices dudx dudy etc with the same dimension of u and v, how it works at the last row? I mean, how is calculated the last row of the derivates matrices if we don't have a subsequent row to calculate the difference U(N+1,:)-U(N,:)?

The edges are estimated using forward/backward difference. Even in lsgradient.m only 2 row/column is used

% Take forward differences on left and right edges if n > 1 g(1,:) = (f(2,:) - f(1,:))/(h(2)-h(1)); g(2,:) = (f(3,:) - f(2,:))/(h(3)-h(2)); g(n-1,:) = (f(n-1,:) - f(n-2,:))/(h(end-1)-h(end-2)); g(n,:) = (f(n,:) - f(n-1,:))/(h(end)-h(end-1)); end

 

2) I've read on the help that, instead of scalars dx and dy, I can give also the vectors x,y with coordinates. --> [dudx,dudy]=lsgradient(u,x',y);

But in the particular case I have x and y as matrices with a little change of dx,dy in all rows/columns, how can I calculate the gradient in the best way? I tried to pass x and y as matrices but the function takes always only the first row/column of x/y and calculate all the gradient with that vector.

if you try an example in lsgradient.m:

        [x,y] = meshgrid(-2:.2:2, -2:.2:2);

        z = x .* exp(-x.^2 - y.^2);

        [px,py] = lsgradient(z,.2,.2);

        contour(z),hold on, quiver(px,py), hold off

and replace the call by using x,y matrices, it also works: 

[x,y] = meshgrid(-2:.2:2, -2:.2:2);

z = x .* exp(-x.^2 - y.^2);

[px,py] = lsgradient(z,x,y); % <---- Note the call with the matrices. 

contour(z),hold on, quiver(px,py), hold off

 

Thanks