Changeset 576 for trunk/src/civ_matlab.m


Ignore:
Timestamp:
Mar 4, 2013, 8:13:53 AM (8 years ago)
Author:
sommeria
Message:

grid improved for civ: computation done closer to the edge. set_grid improved.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/civ_matlab.m

    r575 r576  
    5151check_civ1=0;%default
    5252check_patch1=0;%default
    53 %         ImageFileA=Param.Civ1.ImageA;
    54 %         ImageFileB=Param.Civ1.ImageB;
    55 %         if isfield(Param,'Civ2')
    56 %             ImageFileA_civ2=Param.Civ2.ImageA;
    57 %              ImageFileB_civ2=Param.Civ2.ImageB;
    58 %         end
     53
    5954% case of input Param set by an xml file (batch mode)
    6055if ischar(Param)
    6156    Param=xml2struct(Param); %if Param is the name of an xml file, read this file as a Matlab structure
    62 %     if isfield(Param,'Civ1')
    63 %         if strcmp(Param.Civ1.FileTypeA,'video')
    64 %             Param.Civ1.ImageA=VideoReader(regexprep(Param.Civ1.ImageA,'''','\'));% remove spurious ' appearing in Windows;
    65 %         elseif strcmp(Param.Civ1.FileTypeA,'mmreader')
    66 %             Param.Civ1.ImageA=mmreader(regexprep(Param.Civ1.ImageA,'''','\'));% remove spurious ' appearing in Windows
    67 %         end
    68 %         if strcmp(Param.Civ1.FileTypeB,'video')
    69 %             Param.Civ1.ImageB=VideoReader(regexprep(Param.Civ1.ImageB,'''','\'));% remove spurious ' appearing in Windows
    70 %         elseif strcmp(Param.Civ1.FileTypeB,'mmreader')
    71 %             Param.Civ1.ImageB=mmreader(regexprep(Param.Civ1.ImageB,'''','\'));% remove spurious ' appearing in Windows
    72 %         end
    73 %     end
    74 %     if isfield(Param,'Civ2')
    75 %         if strcmp(Param.Civ2.FileTypeA,'video')
    76 %             Param.Civ2.ImageA=VideoReader(regexprep(Param.Civ2.ImageA,'''','\'));% remove spurious ' appearing in Windows;
    77 %         elseif strcmp(Param.Civ2.FileTypeA,'mmreader')
    78 %             Param.Civ2.ImageA=mmreader(regexprep(Param.Civ2.ImageA,'''','\'));% remove spurious ' appearing in Windows
    79 %             if strcmp(Param.Civ2.FileTypeB,'video')
    80 %                 Param.Civ2.ImageB=VideoReader(regexprep(Param.Civ2.ImageB,'''','\'));% remove spurious ' appearing in Windows
    81 %             elseif strcmp(Param.Civ2.FileTypeB,'mmreader')
    82 %                 Param.Civ2.ImageB=mmreader(regexprep(Param.Civ2.ImageB,'''','\'));% remove spurious ' appearing in Windows
    83 %             end;
    84 %         end
    85 %     end
    8657end
    8758
     
    11687    % caluclate velocity data (y and v in indices, reverse to y component)
    11788    [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
    118    
    119     % to try the reverse_pair method, uncomment below
    120     %     [xtable1 ytable1 utable1 vtable1 ctable1 F1 result_conv1 errormsg1] = civ (Param.Civ1);
    121     %     Param.Civ1.reverse_pair=1;
    122     %     [xtable2 ytable2 utable2 vtable2 ctable2 F2 result_conv2 errormsg2] = civ (Param.Civ1);
    123     %     xtable=[xtable1; xtable2];
    124     %     ytable=[ytable1; ytable2];
    125     %     utable=[utable1; -utable2];
    126     %     vtable=[vtable1; -vtable2];
    127     %     ctable=[ctable1; ctable2];
    128     %     F=[F1; F2];
    129     %     result_conv=[result_conv1; result_conv2];
    130     %     errormsg=[errormsg1; errormsg2];
    13189    if ~isempty(errormsg)
    13290        return
     
    271229        end
    272230    end
    273     ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
    274     iby2=ceil(par_civ2.CorrBoxSize(2)/2);
    275     isx2=ibx2+4;% search ara +-4 pixels around the guess
    276     isy2=iby2+4;
    277     % shift from par_civ2.filename_nc1
    278     % shiftx=velocity interpolated at position
    279     miniy=max(1+isy2,1+iby2);
    280     minix=max(1+isx2,1+ibx2);
    281     maxiy=min(size(par_civ2.ImageA,1)-isy2,size(par_civ2.ImageA,1)-iby2);
    282     maxix=min(size(par_civ2.ImageA,2)-isx2,size(par_civ2.ImageA,2)-ibx2);
    283     [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
    284     GridX=reshape(GridX,[],1);
    285     GridY=reshape(GridY,[],1);
    286     Shiftx=zeros(size(GridX));% shift expected from civ1 data
    287     Shifty=zeros(size(GridX));
    288     nbval=zeros(size(GridX));
     231    %     ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
     232    %     iby2=ceil(par_civ2.CorrBoxSize(2)/2);
     233    %     isx2=ibx2+4;% search ara +-4 pixels around the guess
     234    %     isy2=iby2+4;
     235    %     % shift from par_civ2.filename_nc1
     236    %     % shiftx=velocity interpolated at position
     237    %     miniy=max(1+isy2,1+iby2);
     238    %     minix=max(1+isx2,1+ibx2);
     239    %     maxiy=min(size(par_civ2.ImageA,1)-isy2,size(par_civ2.ImageA,1)-iby2);
     240    %     maxix=min(size(par_civ2.ImageA,2)-isx2,size(par_civ2.ImageA,2)-ibx2);
     241    %     [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
     242    %     GridX=reshape(GridX,[],1);
     243    %     GridY=reshape(GridY,[],1);
     244    if isfield(par_civ2,'Grid')% grid points set as input file
     245        if ischar(par_civ2.Grid)%read the grid file if the input is a file name
     246            par_civ2.Grid=dlmread(par_civ2.Grid);
     247            par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     248        end
     249    else% automatic grid
     250        minix=floor(par_civ2.Dx/2)-0.5;
     251        maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
     252        miniy=floor(par_civ2.Dy/2)-0.5;
     253        maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
     254        [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
     255        par_civ2.Grid(:,1)=reshape(GridX,[],1);
     256        par_civ2.Grid(:,2)=reshape(GridY,[],1);
     257    end
     258    Shiftx=zeros(size(par_civ2.Grid,1));% shift expected from civ1 data
     259    Shifty=zeros(size(par_civ2.Grid,1));
     260    nbval=zeros(size(par_civ2.Grid,1));
    289261    if par_civ2.CheckDeformation
    290         DUDX=zeros(size(GridX));
    291         DUDY=zeros(size(GridX));
    292         DVDX=zeros(size(GridX));
    293         DVDY=zeros(size(GridX));
     262        DUDX=zeros(size(par_civ2.Grid,1));
     263        DUDY=zeros(size(par_civ2.Grid,1));
     264        DVDX=zeros(size(par_civ2.Grid,1));
     265        DVDY=zeros(size(par_civ2.Grid,1));
    294266    end
    295267    NbSubDomain=size(Data.Civ1_SubRange,3);
     
    316288        mask=imread(par_civ2.Mask);
    317289    end
    318     par_civ2.SearchBoxSize(1)=2*isx2+1;
    319     par_civ2.SearchBoxSize(2)=2*isy2+1;
     290    ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
     291    iby2=ceil(par_civ2.CorrBoxSize(2)/2);
     292    %     isx2=ibx2+4;% search ara +-4 pixels around the guess
     293    %     isy2=iby2+4;
     294    par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess
     295    par_civ2.SearchBoxSize(2)=2*iby2+9;
     296    %par_civ2.SearchBoxSize(1)=2*isx2+1;
     297    %par_civ2.SearchBoxSize(2)=2*isy2+1;
    320298    par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
    321 %     par_civ2.SearchBoxShift(2)=Shifty(nbval>=1)./nbval(nbval>=1);
    322299    par_civ2.Grid=[GridX(nbval>=1)-par_civ2.SearchBoxShift(:,1)/2 GridY(nbval>=1)-par_civ2.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors
    323300    if par_civ2.CheckDeformation
     
    329306    % caluclate velocity data (y and v in indices, reverse to y component)
    330307    [xtable ytable utable vtable ctable F] = civ (par_civ2);
    331     %     diff_squared=(utable-par_civ2.Shiftx).*(utable-par_civ2.Shiftx)+(vtable+par_civ2.Shifty).*(vtable+par_civ2.Shifty);
    332     %     F(diff_squared>=4)=4; %flag vectors whose distance to the guess exceeds 2 pixels
    333308    list_param=(fieldnames(Param.Civ2))';
    334309    list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'};
     
    350325    end
    351326    Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param {'Civ2_Time','Civ2_Dt'}];
    352     %     Data.Civ2_Time=par_civ2.Time;
    353     %     Data.Civ2_Dt=par_civ2.Dt;
    354327    nbvar=numel(Data.ListVarName);
    355328    Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C'}];%  cell array containing the names of the fields to record
     
    475448%methods perform better and will maybe be implemented in the future.
    476449
    477 %% prepare grid
     450%% prepare measurement grid
     451if isfield(par_civ,'Grid')% grid points set as input
     452    if ischar(par_civ.Grid)%read the drid file if the input is a file name
     453        par_civ.Grid=dlmread(par_civ.Grid);
     454        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     455    end
     456else% automatic grid
     457    minix=floor(par_civ.Dx/2)-0.5;
     458    maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
     459    miniy=floor(par_civ.Dy/2)-0.5;
     460    maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
     461    [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     462    par_civ.Grid(:,1)=reshape(GridX,[],1);
     463    par_civ.Grid(:,2)=reshape(GridY,[],1);
     464end
     465nbvec=size(par_civ.Grid,1);
     466
     467%% prepare correlation and search boxes
    478468ibx2=ceil(par_civ.CorrBoxSize(1)/2);
    479469iby2=ceil(par_civ.CorrBoxSize(2)/2);
     
    482472shiftx=round(par_civ.SearchBoxShift(:,1));
    483473shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
    484 if isfield(par_civ,'Grid')
    485     if ischar(par_civ.Grid)%read the drid file if the input is a file name
    486         par_civ.Grid=dlmread(par_civ.Grid);
    487         par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
    488     end
    489 else% automatic measurement grid
    490 %     ibx2=ceil(par_civ.Bx/2);
    491 %     iby2=ceil(par_civ.By/2);
    492 %     isx2=ceil(par_civ.Searchx/2);
    493 %     isy2=ceil(par_civ.Searchy/2);
    494     miniy=max(1+isy2+shifty,1+iby2);
    495     minix=max(1+isx2-shiftx,1+ibx2);
    496     maxiy=min(par_civ.ImageHeight-isy2+shifty,par_civ.ImageHeight-iby2);
    497     maxix=min(par_civ.ImageWidth-isx2-shiftx,par_civ.ImageWidth-ibx2);
    498     [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
    499     par_civ.Grid(:,1)=reshape(GridX,[],1);
    500     par_civ.Grid(:,2)=reshape(GridY,[],1);
    501 end
    502 nbvec=size(par_civ.Grid,1);
    503 if numel(shiftx)==1
     474if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
    504475    shiftx=shiftx*ones(nbvec,1);
    505476    shifty=shifty*ones(nbvec,1);
    506477end
     478
    507479%% Default output
    508480xtable=par_civ.Grid(:,1);
     
    561533    %  20>=mask: velocity=0
    562534checkmask=0;
     535MinA=min(min(par_civ.ImageA));
     536MinB=min(min(par_civ.ImageB));
    563537if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
    564538   checkmask=1;
     
    569543  %  check_noflux=(par_civ.Mask<100) ;%TODO: to implement
    570544    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
    571     par_civ.ImageA(check_undefined)=min(min(par_civ.ImageA));% put image A to zero (i.e. the min image value) in the undefined  area
    572     par_civ.ImageB(check_undefined)=min(min(par_civ.ImageB));% put image B to zero (i.e. the min image value) in the undefined  area
     545    par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined  area
     546    par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined  area
    573547end
    574548
     
    587561    jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
    588562    if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
    589         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
    590               jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
    591             F(ivec)=3;
    592         else
    593             image1_crop=par_civ.ImageA(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2);%extract a subimage (correlation box) from image A
    594             image2_crop=par_civ.ImageB(jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2,iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2);%extract a larger subimage (search box) from image B
     563%         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
     564%               jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
     565%             F(ivec)=3;
     566%         else
     567            F(ivec)=0;
     568            subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
     569            subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
     570            subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
     571            subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
     572            image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
     573            image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
     574            check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
     575            check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
     576            check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
     577            check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
     578            image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A
     579            image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B
    595580            image1_mean=mean(mean(image1_crop));
    596581            image2_mean=mean(mean(image2_crop));
     
    603588                F(ivec)=3;
    604589            end
    605         end     
     590%         end     
    606591        if F(ivec)~=3
    607592            image1_crop=image1_crop-image1_mean;%substract the mean
Note: See TracChangeset for help on using the changeset viewer.