image thumbnail

mpiv

by

 

20 Sep 2002 (Updated )

PIV method in MATLAB

piv_mrs( im1, im2, ...
function [xi_2,yi_2, iu_2,iv_2] = piv_mrs( im1, im2, ... 
				     xi_1, yi_1, iu_1, iv_1, ...
				     nx_pixel, ny_pixel, n_scale, ...
				     overlap_x, overlap_y, ...
				     iu_max, iv_max, ...
				     i_filter, vec_std, i_interp, ...
				     ir, n_recur )
%========================================================================
%
% version 0.10
%
%
% 	piv_mrs.m
%
%
% Description:
%
%	This program is a part of
%	calculation of PIV displacement using MQD method 
%	with hierarchical (iterative) scheme.
%	This program is to be called by the main program 'piv_mqr.m'
%	  - simple MQD algorism
%	  - hierarchical/recursive algorism
%	  - search close to image edge
%	  - Gaussian subpixel fit
%
% Requires:
%	  - piv_mqd.m
%	  - piv_mqr.m
%	  - func_findpeak2.m
%	  - piv_windowsize.m
%	  - vector_check.m
%	  - vector_filter_linear.m
%	  - vector_filter_median.m
%	  - vector_filter_vecstd.m
%	  - vector_exterp_linear.m
%	  - vector_interp.m
%	  - vector_interp_kriging.m
%	  - vector_interp_kriging_local.m
%	  - vector_interp_linear.m
%	  - vector_interp_spline.m
%
% Requires toolboxes:
%
%   - This program requires DACE, Kriging Toolbox, developed by
%	S.N. Lophaven, H.B. Nielsen and J. Sondergaard
%	at Technical University of Denmark
%
% Variables:
%
%  Input:
%	im1 and 2	image files (double precision)
%	nx_pixel	subwindow size in x
%			(should be larger than 20, typical 32 or 64)
%	ny_pixel	subwindow size in y
%			(should be larger than 20, typical 32 or 64)
%	overlap_x	overlap ratio of adjacent subwindows in x
%	overlap_y	overlap ratio of adjacent subwindows in y
%	iu_max		maximum displacement in x (unit: pixel)
%	iv_max		maximum displacement in y (unit: pixel)
%			-> iu_max and iv_max give search area
%	i_opt		= -1: no subpixel fit
%			=  1: with Gaussian subpixel fit
%
%  Output:
%	xi, yi		location of velocity vector
%	iu, iv		velocity vector
%	D		maximum value of image intensity difference
%
%======================================================================
%
% Terms:
%
%       Distributed under the terms of the terms of the BSD License
%
% Copyright:
%
%       Nobuhito Mori
%           Disaster Prevention Research Institue
%           Kyoto University, JAPAN
%           mori@oceanwave.jp
%
%======================================================================
%
% Update:
%       0.11    2009/07/01 BSD License applied
%	0.10	2003/07/02 MMR has been used to eliminate error vector
%	0.01	2003/07/01 First version
%
%========================================================================

%
% <===== Input for professional use				<=====
%      can be replaced with other appropriate values

% i_plot = 1        check plotting
%        = other    no plotting
i_plot = 9;

% SNR: signal to noise ratio to find peak
r_SNR = 3.0;

% r_peak: ratio of maximum and mean
r_MMR = 1.10;

% r_peak: ratio of 1st peak and 2nd peak
r_PPR = 1.05;

% set area of search
p_search = 1/3; 			% percentage of subwindow

% set minimum subwindow size for interrogation
n_window_min = 9;                       % pixel

% min and max values for MQD
d_min      = 10^(-5);
d_max      = Inf;

% =====> End of input for professional use			======>

%
% --- initialization
%

nx = size(im1,1);
ny = size(im1,2);

% window size
mx_pixel = nx_pixel/n_scale;
my_pixel = ny_pixel/n_scale;

[ xi_2,yi_2,mx_start,my_start,mx_overlap,my_overlap,dx_center, dy_center ] ...
    = func_pivwindowsize( 'mqd', ...
			  nx,ny, mx_pixel,my_pixel, overlap_x,overlap_y );

lx = max(size(xi_1));
ly = max(size(yi_1));
mx = max(size(xi_2));
my = max(size(yi_2));

disp(' ========================================')
c_tmp = strcat( '  Total number of vectors = ', ...
		int2str(mx),' x ',int2str(my) );
disp( c_tmp );
c_tmp = strcat( '  Interrogation window size = ', ...
		int2str(mx_pixel),' x ',int2str(my_pixel) );
disp( c_tmp );
disp(' ========================================')

%
% --- preprocess to interpolate velocity vector
%

[X1 Y1] = meshgrid(xi_1,yi_1);
[X2 Y2] = meshgrid(xi_2,yi_2);

[ ut, vt, i_cond ] = vector_check( iu_1, iv_1, vec_std, i_filter ); 

if i_interp == 3
  [ u ] = vector_interp_kriging_local( ut );
  [ v ] = vector_interp_kriging_local( vt ); 
else
  [ u ] = vector_interp_linear( ut );
  [ v ] = vector_interp_linear( vt ); 
  [ u ] = vector_exterp_linear( ut );
  [ v ] = vector_exterp_linear( vt ); 
end

% remove NaN vector
if i_interp ~= 3
  while max(max(isnan(u))) == 1
    [ u ] = vector_interp_linear( u );
    [ u ] = vector_exterp_linear( u );
    [ u ] = vector_interp_NaN( u );
  end  
  while max(max(isnan(v))) == 1
    [ v ] = vector_interp_linear( v );
    [ v ] = vector_exterp_linear( v );
    [ v ] = vector_interp_NaN( v );
  end
end

%
% --- predicting velocity vector using interp2 interpolation 
%

Ui = interp2( X1,Y1, u', X2, Y2, '*spline' );
Vi = interp2( X1,Y1, v', X2, Y2, '*spline' );
Ui = Ui';
Vi = Vi';

[ Ui, Vi, i_cond ] = vector_check( Ui, Vi, vec_std, i_filter ); 

if i_interp == 3
  [ UI ] = vector_interp_kriging_local( Ui );
  [ VI ] = vector_interp_kriging_local( Vi ); 
else
  [ UI ] = vector_exterp_linear( Ui );
  [ VI ] = vector_exterp_linear( Vi );
end

% test
if i_plot == 1
  figure(2);clf
  quiver(X1,Y1,u',v','r-');
  hold on
  quiver(X2,Y2,UI',VI','b-');
  quiver(X2,Y2,Ui',Vi','g-');
  hold off
  pause
end

%
% --- main routine
%

  for iy = 1:my

    c_proc = strcat( 'Recursive PIV step', ...
		   int2str(ir),'/',int2str(n_recur), ...
		   '.  Process accomplished : ', ...
		   num2str( 100*(iy-1)/(my-1),' %03.0f' ), '/100' );
    disp( c_proc )

    for ix = 1:mx

      % create target window from 1st image
      % location of target window
      ix1 = xi_2(ix) - dx_center;
      % ix1 = mx_start + (mx_pixel - mx_overlap)*(ix-1) + 1;
      ix2 = ix1 + mx_pixel - 1;
      iy1 = yi_2(iy) - dy_center;
      % iy1 = my_start + (my_pixel - my_overlap)*(iy-1) + 1;
      iy2 = iy1 + my_pixel - 1;

      % center of target window
      ix_center = xi_2(ix);
      iy_center = yi_2(iy);

      f1 = im1( ix1:ix2, iy1:iy2 );

      %
      % --- set search area
      %

      % search area must be larger than n_wind_min or 1+p_seach percentage
      mx_move = round( [ (1-p_search), (1+p_search) ]*UI(ix,iy) );
      my_move = round( [ (1-p_search), (1+p_search) ]*VI(ix,iy) );
      mx_d = abs( max(mx_move)-min(mx_move) );
      my_d = abs( max(my_move)-min(my_move) );
      lx = min( fix(mx_pixel/2)+1, n_window_min);
      ly = min( fix(my_pixel/2)+1, n_window_min);
      
      if mx_d > lx
        mx_movemin = min(mx_move);
        mx_movemax = max(mx_move);
      else
        mx_movemin = round(UI(ix,iy)) - ceil(lx/2);
        mx_movemax = round(UI(ix,iy)) + ceil(ly/2);
      end
      
      if my_d > ly
        my_movemin = min(my_move);
        my_movemax = max(my_move);
      else
        my_movemin = round(VI(ix,iy)) - fix(lx/2);
        my_movemax = round(VI(ix,iy)) + fix(ly/2);
      end

      % search area must be odd number
      if rem(mx_movemax-mx_movemin+1,2) == 0
        mx_movemax = mx_movemax + 1;
      end
      
      if rem(my_movemax-my_movemin+1,2) == 0
        my_movemax = my_movemax + 1;
      end

      % check the value of movemax
      if isnan(mx_movemax*mx_movemin*my_movemax*my_movemin) == 1
        mx_movemin = 0;
        mx_movemax = 0;
        my_movemin = 0;
        my_movemax = 0;
      end

      clear C;

      %
      % --- calculation uisng MQD
      %

      isy = 0;
      for jy = my_movemin:my_movemax
        isx = 0;
        isy = isy + 1;

        for jx = mx_movemin:mx_movemax
          isx = isx + 1;

          C(isx,isy) = d_min;

          %
          % --- create search subwindow in 2nd image
          %

          f2 = zeros(mx_pixel,my_pixel);

          kx1 = ix1 + jx;
          kx2 = kx1 + mx_pixel - 1;
          ky1 = iy1 + jy;
          ky2 = ky1 + my_pixel - 1;
          n   = mx_pixel*my_pixel;

          if ( kx1 >= 1 ) & ( kx2 <= nx ) & ( ky1 >= 1 ) & ( ky2 <= ny ) 
     
            f2 = im2( kx1:kx2, ky1:ky2 );   
	     d = sum( sum( abs(f1 - f2) ) )/n;

	    if d ~= 0
	      C(isx,isy) = d;
            else
	      C(isx,isy) = d_min;
            end

	  else
    	     C(isx,isy) = d_max;
          end

        end
      end

     %
     % --- peak finding /  C is odd number maxtrix
     %

     C = 1./C;

     % Peak finding with Gaussian subpixel fit
     [ ip_x, ip_y, SNR, MMR, PPR ] = func_findpeak2( C, -2 );

     ix_peak = mx_movemin + ip_x - 1;
     iy_peak = my_movemin + ip_y - 1;

     if i_plot == 1
       figure(1);clf;colormap(jet)
       %contour( C', 25 );
       surf(C');
       hold on
       C_max = max(max(C));
       plot3( ip_x, ip_y, C_max, 'bo', ...
	         'MarkerSize', 8, ...
	         'MarkerFaceColor', 'b', ...
	         'MarkerEdgeColor', 'k' );
       hold off
       disp( '-> push key' )
       pause
     end

     %
     % --- eliminate vectors with small MMR and PPR
     %

     if MMR < r_MMR | PPR < r_PPR
       ix_peak = NaN;
       iy_peak = NaN;
     end

     %
     % --- eliminate vectors in large displacement with iu_max and iv_max
     %

     if abs(ix_peak) >= iu_max
       ix_peak = NaN;
       iy_peak = NaN;
     end
     if abs(iy_peak) >= iv_max
       ix_peak = NaN;
       iy_peak = NaN;
     end

     %
     % --- end of main loop
     %

     is_x(ix,iy) = ix_peak;
     is_y(ix,iy) = iy_peak;

    end
  end 

%
% --- post process
%

iu_2 = is_x;
iv_2 = is_y;

Contact us