Changeset 1115 for trunk/src


Ignore:
Timestamp:
Apr 6, 2022, 8:36:16 PM (3 years ago)
Author:
sommeria
Message:

geometry calib generalised to order_4 and proj_field debugged

Location:
trunk/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/find_field_cells.m

    r1107 r1115  
    162162        cell_counter=cell_counter+1;
    163163        icell=cell_counter;
     164        CellInfo{icell}.FieldName=Data.ListVarName{ind_scalar(iscalar)};
    164165        CellInfo{icell}.VarType='scalar';
    165166        DimCell{icell}=Data.VarDimName{ind_scalar(iscalar)};

  • trunk/src/geometry_calib.m

    r1113 r1115  
    122122
    123123%% set menu of calibration options
    124 set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_extrinsic'})
     124set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_order4';'3D_extrinsic'})
    125125if exist('inputfile','var')&& ~isempty(inputfile)
    126126    [RootPath,SubDir,RootFile,tild,tild,tild,tild,FileExt]=fileparts_uvmat(inputfile);
    … …  
    232232
    233233%% Apply calibration
    234 [GeometryCalib,index,ind_removed]=calibrate(Coord,CalibFcn,Intrinsic);% apply calibration
     234est_dist=[0;0;0;0;0]; %default
     235switch CalibFcn
     236    case 'calib_3D_linear'
     237        CalibFcn='calib_3D';
     238    case 'calib_3D_quadr'
     239        est_dist=[1;0;0;0;0];
     240        CalibFcn='calib_3D'
     241    case 'calib_3D_order4'
     242        est_dist=[1;1;0;0;0];
     243        CalibFcn='calib_3D'
     244end
     245
     246[GeometryCalib,index,ind_removed]=calibrate(Coord,CalibFcn,est_dist,Intrinsic);% apply calibration
    235247if isempty(GeometryCalib)
    236248    return
    … …  
    393405%------------------------------------------------------------------------
    394406% --- activate calibration and store parameters in ouputfile .
    395 function [GeometryCalib,ind_max,ind_removed]=calibrate(Coord,CalibFcn,Intrinsic)
     407function [GeometryCalib,ind_max,ind_removed]=calibrate(Coord,CalibFcn,est_dist,Intrinsic)
    396408%------------------------------------------------------------------------
    397409
    … …  
    403415CoordFlip(:,3)=Coord(:,3);% the calibration function assume a z ccordinate along the camera view, opposite to ours
    404416if ~isempty(Coord)
    405     GeometryCalib=feval(CalibFcn,CoordFlip,Intrinsic);
     417    GeometryCalib=feval(CalibFcn,CoordFlip,est_dist,Intrinsic);
    406418else
    407419    msgbox_uvmat('ERROR','No calibration points, abort')
    … …  
    433445    CoordFlip=Coord;
    434446    CoordFlip(:,3)=Coord(:,3);% the calibration function assume a z ccordinate along the camera view, opposite to ours
    435     GeometryCalib=feval(CalibFcn,CoordFlip,Intrinsic);
     447    GeometryCalib=feval(CalibFcn,CoordFlip,est_dist,Intrinsic);
    436448    X=Coord(:,1);
    437449    Y=Coord(:,2);
    … …  
    448460%------------------------------------------------------------------------
    449461% --- determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
    450 function GeometryCalib=calib_rescale(Coord,Intrinsic)
     462function GeometryCalib=calib_rescale(Coord,est_dist,Intrinsic)
    451463%------------------------------------------------------------------------
    452464X=Coord(:,1);
    … …  
    464476%------------------------------------------------------------------------
    465477% --- determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
    466 function GeometryCalib=calib_linear(Coord,Intrinsic)
     478function GeometryCalib=calib_linear(Coord,est_dist,Intrinsic)
    467479%------------------------------------------------------------------------
    468480X=Coord(:,1);
    … …  
    494506% --- determine the tsai parameters for a view normal to the grid plane
    495507% NOT USED
    496 function GeometryCalib=calib_normal(Coord,Intrinsic)
     508function GeometryCalib=calib_normal(Coord,est_dist,Intrinsic)
    497509%------------------------------------------------------------------------
    498510Calib.fx=str2num(get(handles.fx,'String'));
    … …  
    555567alpha=calib_param(5);
    556568GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
    557 
    558 %------------------------------------------------------------------------
    559 function GeometryCalib=calib_3D_linear(Coord,Intrinsic)
    560 %------------------------------------------------------------------------
    561 coord_files=Intrinsic.coord_files;
    562 if ischar(Intrinsic.coord_files)
    563     coord_files={coord_files};
    564 end
    565 if isempty(coord_files{1}) || isequal(coord_files,{''})
    566     coord_files={};
    567 end
    568 %retrieve the calibration points stored in the files listed in the popup list ListCoordFiles
    569 x_1=Coord(:,4:5)';%px coordinates of the ref points
    570 
    571 nx=Intrinsic.Npx;
    572 ny=Intrinsic.Npy;
    573 x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
    574 X_1=Coord(:,1:3)';%phys coordinates of the ref points
    575 n_ima=numel(coord_files)+1;
    576 if ~isempty(coord_files)
    577     msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
    578     for ifile=1:numel(coord_files)
    579         t=xmltree(coord_files{ifile});
    580         s=convert(t);%convert to matlab structure
    581         if isfield(s,'GeometryCalib')
    582             if isfield(s.GeometryCalib,'SourceCalib')
    583                 if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
    584                     PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
    585                     Coord_file=zeros(length(PointCoord),5);%default
    586                     for i=1:length(PointCoord)
    587                         line=str2num(PointCoord{i});
    588                         Coord_file(i,4:5)=line(4:5);%px x
    589                         Coord_file(i,1:3)=line(1:3);%phys x
    590                     end
    591                     eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
    592                     eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
    593                     eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
    594                 end
    595             end
    596         end
    597     end
    598 end
    599 n_ima=numel(coord_files)+1;
    600 est_dist=[0;0;0;0;0];
    601 est_aspect_ratio=0;
    602 est_fc=[1;1];
    603 center_optim=0;
    604 path_uvmat=which('uvmat');% check the path detected for source file uvmat
    605 path_UVMAT=fileparts(path_uvmat); %path to UVMAT
    606 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));% apply fct 'toolbox_calib/go_calib_optim'
    607 if exist('Rc_1','var')
    608     GeometryCalib.CalibrationType='3D_extrinsic';
    609     GeometryCalib.fx_fy=fc';
    610     GeometryCalib.Cx_Cy=cc';
    611     GeometryCalib.kc=kc(1);
    612     GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    613     GeometryCalib.Tx_Ty_Tz=Tc_1';
    614     GeometryCalib.R=Rc_1;
    615     GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
    616     GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
    617     GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
    618     GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees
    619     GeometryCalib.ErrorRMS=[];
    620     GeometryCalib.ErrorMax=[];
    621 else
    622     msgbox_uvmat('ERROR',['calibration function ' fullfile('toolbox_calib','go_calib_optim') ' did not converge: use multiple views or option 3D_extrinsic'])
    623     GeometryCalib=[];
    624 end
    625 
    626 %------------------------------------------------------------------------
    627 function GeometryCalib=calib_3D_quadr(Coord,Intrinsic)
     569%
     570% %------------------------------------------------------------------------
     571% function GeometryCalib=calib_3D_linear(Coord,Intrinsic)
     572% %------------------------------------------------------------------------
     573% coord_files=Intrinsic.coord_files;
     574% if ischar(Intrinsic.coord_files)
     575%     coord_files={coord_files};
     576% end
     577% if isempty(coord_files{1}) || isequal(coord_files,{''})
     578%     coord_files={};
     579% end
     580% %retrieve the calibration points stored in the files listed in the popup list ListCoordFiles
     581% x_1=Coord(:,4:5)';%px coordinates of the ref points
     582%
     583% nx=Intrinsic.Npx;
     584% ny=Intrinsic.Npy;
     585% x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
     586% X_1=Coord(:,1:3)';%phys coordinates of the ref points
     587% n_ima=numel(coord_files)+1;
     588% if ~isempty(coord_files)
     589%     msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
     590%     for ifile=1:numel(coord_files)
     591%         t=xmltree(coord_files{ifile});
     592%         s=convert(t);%convert to matlab structure
     593%         if isfield(s,'GeometryCalib')
     594%             if isfield(s.GeometryCalib,'SourceCalib')
     595%                 if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
     596%                     PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
     597%                     Coord_file=zeros(length(PointCoord),5);%default
     598%                     for i=1:length(PointCoord)
     599%                         line=str2num(PointCoord{i});
     600%                         Coord_file(i,4:5)=line(4:5);%px x
     601%                         Coord_file(i,1:3)=line(1:3);%phys x
     602%                     end
     603%                     eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
     604%                     eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
     605%                     eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
     606%                 end
     607%             end
     608%         end
     609%     end
     610% end
     611% n_ima=numel(coord_files)+1;
     612% est_dist=[0;0;0;0;0];
     613% est_aspect_ratio=0;
     614% est_fc=[1;1];
     615% center_optim=0;
     616% path_uvmat=which('uvmat');% check the path detected for source file uvmat
     617% path_UVMAT=fileparts(path_uvmat); %path to UVMAT
     618% run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));% apply fct 'toolbox_calib/go_calib_optim'
     619% if exist('Rc_1','var')
     620%     GeometryCalib.CalibrationType='3D_extrinsic';
     621%     GeometryCalib.fx_fy=fc';
     622%     GeometryCalib.Cx_Cy=cc';
     623%     GeometryCalib.kc=kc(1);
     624%     GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     625%     GeometryCalib.Tx_Ty_Tz=Tc_1';
     626%     GeometryCalib.R=Rc_1;
     627%     GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
     628%     GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
     629%     GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
     630%     GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees
     631%     GeometryCalib.ErrorRMS=[];
     632%     GeometryCalib.ErrorMax=[];
     633% else
     634%     msgbox_uvmat('ERROR',['calibration function ' fullfile('toolbox_calib','go_calib_optim') ' did not converge: use multiple views or option 3D_extrinsic'])
     635%     GeometryCalib=[];
     636% end
     637
     638%------------------------------------------------------------------------
     639function GeometryCalib=calib_3D(Coord,est_dist,Intrinsic)
    628640%------------------------------------------------------------------
    629641
    … …  
    678690end
    679691n_ima=numel(coord_files)+1;
    680 est_dist=[1;0;0;0;0];
     692%est_dist=[1;0;0;0;0];
    681693est_aspect_ratio=1;
    682694center_optim=0;
    … …  
    687699    GeometryCalib.fx_fy=fc';
    688700    GeometryCalib.Cx_Cy=cc';
    689     GeometryCalib.kc=kc(1);
     701    ind_kc=find(kc);
     702    if ~isempty(ind_kc)
     703    GeometryCalib.kc=(kc(ind_kc))';%include cases quadr (one value kc(1)) and order_4 (2 values for kc)
     704    end
    690705    GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    691706    GeometryCalib.Tx_Ty_Tz=Tc_1';
    … …  
    703718
    704719%------------------------------------------------------------------------
    705 function GeometryCalib=calib_3D_extrinsic(Coord, Intrinsic)
     720function GeometryCalib=calib_3D_extrinsic(Coord,est_dist, Intrinsic)
    706721%------------------------------------------------------------------
    707722path_uvmat=which('geometry_calib');% check the path detected for source file uvmat
    … …  
    747762addpath(fct_path)
    748763GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%reverse Cx_Cy(2) for calibration (inversion of px ordinate)
     764kc=zeros(1,5);
     765kc(1:numel(GeometryCalib.kc))=GeometryCalib.kc;
    749766[omc,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
    750     (GeometryCalib.fx_fy)',GeometryCalib.Cx_Cy',[GeometryCalib.kc 0 0 0 0]);
     767    (GeometryCalib.fx_fy)',GeometryCalib.Cx_Cy',kc);
    751768rmpath(fct_path);
    752 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     769GeometryCalib.CoordUnit='';% default value, to be updated by the calling function
    753770GeometryCalib.Tx_Ty_Tz=Tc1';
    754771%inversion of z axis
    … …  
    10861103Tmod(:,2)=(TIndex(:,2)-1)*Dy+Rangy(1);
    10871104%Tmod=T(:,(1:2))+Delta;% 'phys' coordinates of the detected points
    1088 [Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));% image coordinates of the detected points
     1105[Xpx,Ypx]=px_XYZ(GeometryCalib,[],Tmod(:,1),Tmod(:,2));% image coordinates of the detected points
    10891106Coord=[T Xpx Ypx zeros(size(T,1),1)];
    10901107set(handles.ListCoord,'Data',Coord)
    … …  
    13261343set(handles.Cx,'String',num2str(Cx_Cy(1),'%1.1f'))
    13271344set(handles.Cy,'String',num2str(Cx_Cy(2),'%1.1f'))
    1328 set(handles.kc,'String',num2str(kc,'%1.4f'))
     1345%set(handles.kc,'String',num2str(kc,'%1.4f'))
     1346set(handles.kc,'String',num2str(kc,4))
    13291347
    13301348%------------------------------------------------------------------------

  • trunk/src/phys_XYZ.m

    r1113 r1115  
    7474    Calib.Cx_Cy=[0 0];
    7575end
    76 if ~isfield(Calib,'kc')
    77     Calib.kc=0;
     76kc=[0 0];
     77if isfield(Calib,'kc')
     78    kc(1:numel(Calib.kc))=Calib.kc;
    7879end
     80% if ~isfield(Calib,'kc2')
     81%     Calib.kc2=0;
     82% end
    7983if isfield(Calib,'R')
    8084    R=(Calib.R)';
    81 %     R(3)=-R(3);
    82 %     R(6)=-R(6);
    83 %     R(9)=-R(9);
    8485    c=Z0virt;
    8586    cvirt=Z0virt;
    … …  
    126127    Xd=(X-Calib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates
    127128    Yd=(Y-Calib.Cx_Cy(2))/Calib.fx_fy(2);
    128     dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd);% distortion factor, first approximation Xu,Yu=Xd,Yd
     129    dist_fact=1+kc(1)*(Xd.*Xd+Yd.*Yd);% distortion factor, first approximation Xu,Yu=Xd,Yd
    129130    test=0;
    130131    niter=0;
    … …  
    133134        Xu=Xd./dist_fact;%undistorted sensor coordinates, second iteration
    134135        Yu=Yd./dist_fact;
    135         dist_fact=1+Calib.kc*(Xu.*Xu+Yu.*Yu);% distortion factor,next approximation
     136        R2=Xu.*Xu+Yu.*Yu;
     137        dist_fact=1+kc(1)*R2+kc(2)*R2.*R2;% distortion factor,next approximation
    136138        test=max(max(abs(dist_fact-dist_fact_old)))<0.00001; % reducing the relative error to 10^-5 forthe inversion of the quadraticcorrection
    137139        niter=niter+1;

  • trunk/src/proj_field.m

    r1114 r1115  
    10951095            end
    10961096        else
    1097             ProjData.ProjObjectAngle=FieldData.PlaneAngle;
     1097            ProjData.ProjObjectAngle=FieldData.PlaneAngle;ProjMode
    10981098        end
    10991099    end
    … …  
    11251125icell_grid=[];% field cell index which defines the grid
    11261126icell_scattered=[];% field cell index which defines fields with scattered coordinates
    1127 for icell=1:numel(CellInfo)
    1128     if strcmp(ObjectData.ProjMode,'interp_lin')&& ~strcmp(ProjMode{icell},'interp_tps')
    1129         errormsg='ProjMode interp_tps needed ';
    1130         return
    1131     end
    1132 end
     1127% for icell=1:numel(CellInfo)
     1128%     if strcmp(ObjectData.ProjMode,'interp_lin')&& ~strcmp(ProjMode{icell},'interp_tps')
     1129%         errormsg='ProjMode interp_tps needed ';
     1130%         return
     1131%     end
     1132% end
    11331133if strcmp(ObjectData.ProjMode,'projection')
    11341134    %% case of a grid requested by the input field

  • trunk/src/px_XYZ.m

    r1113 r1115  
    6969    if ~isfield(Calib,'kc')
    7070        r2=1; %no quadratic distortion
     71    elseif numel(Calib.kc)==1
     72        r2=1+Calib.kc*(Xu.*Xu+Yu.*Yu);
    7173    else
    72         r2=1+Calib.kc*(Xu.*Xu+Yu.*Yu);
     74        R2=Xu.*Xu+Yu.*Yu;
     75        r2=1+Calib.kc(1)*R2+Calib.kc(2)*R2.*R2;
    7376    end
    7477   

  • trunk/src/series/civ_series.m

    r1107 r1115  
    628628            ind_good=1:numel(Data.Civ1_X);
    629629        end
    630        
     630        if isempty(ind_good)
     631                        disp_uvmat('ERROR','all vectors of civ1 are bad, check input parameters' ,checkrun)
     632                        return
     633        end
     634
    631635        % perform Patch calculation using the UVMAT fct 'filter_tps'
    632636        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
    … …  
    919923        else
    920924            ind_good=1:numel(Data.Civ2_X);
     925        end
     926                if isempty(ind_good)
     927                        disp_uvmat('ERROR','all vectors of civ2 are bad, check input parameters' ,checkrun)
     928                        return
    921929        end
    922930        [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=...

  • trunk/src/series/merge_proj.m

    r1114 r1115  
    8686    return
    8787end
    88 
     88if 0==1
     89    phys; % used to include phys when compiling is done
     90end
    8991%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
    9092ParamOut=[]; %default output

  • trunk/src/series/turb_stat.m

    r1107 r1115  
    204204% for i_slice=1:Param.IndexRange.NbSlice
    205205%     i_slice
    206     ind_first=Param.IndexRange.first_i
     206    ind_first=Param.IndexRange.first_i;
    207207    for index_i=ind_first:Param.IndexRange.NbSlice:Param.IndexRange.last_i
    208208        if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue')
    … …  
    210210            break
    211211        end
    212         for index_j=Param.IndexRange.first_j:Param.IndexRange.last_j
     212        for index_j=first_j:last_j
    213213            InputFile=fullfile_uvmat(RootPath{1},SubDir{1},RootFile{1},FileExt{1},NomType{1},index_i,index_i,index_j,index_j);
    214214            [Field,tild,errormsg] = read_field(InputFile,FileType{iview},InputFields{iview});
    … …  
    217217           
    218218            %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    219             if index_i==ind_first && index_j==Param.IndexRange.first_j %initiate the output data structure in the first field
     219            if index_i==ind_first && index_j==first_j %initiate the output data structure in the first field
    220220                [CellInfo,NbDim,errormsg]=find_field_cells(Field);
    221221                YName='coord_y';%default
Note: See TracChangeset for help on using the changeset viewer.