Index: trunk/src/find_field_cells.m =================================================================== --- trunk/src/find_field_cells.m (revision 1114) +++ trunk/src/find_field_cells.m (revision 1115) @@ -162,4 +162,5 @@ cell_counter=cell_counter+1; icell=cell_counter; + CellInfo{icell}.FieldName=Data.ListVarName{ind_scalar(iscalar)}; CellInfo{icell}.VarType='scalar'; DimCell{icell}=Data.VarDimName{ind_scalar(iscalar)}; Index: trunk/src/geometry_calib.m =================================================================== --- trunk/src/geometry_calib.m (revision 1114) +++ trunk/src/geometry_calib.m (revision 1115) @@ -122,5 +122,5 @@ %% set menu of calibration options -set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_extrinsic'}) +set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_order4';'3D_extrinsic'}) if exist('inputfile','var')&& ~isempty(inputfile) [RootPath,SubDir,RootFile,tild,tild,tild,tild,FileExt]=fileparts_uvmat(inputfile); @@ -232,5 +232,17 @@ %% Apply calibration -[GeometryCalib,index,ind_removed]=calibrate(Coord,CalibFcn,Intrinsic);% apply calibration +est_dist=[0;0;0;0;0]; %default +switch CalibFcn + case 'calib_3D_linear' + CalibFcn='calib_3D'; + case 'calib_3D_quadr' + est_dist=[1;0;0;0;0]; + CalibFcn='calib_3D' + case 'calib_3D_order4' + est_dist=[1;1;0;0;0]; + CalibFcn='calib_3D' +end + +[GeometryCalib,index,ind_removed]=calibrate(Coord,CalibFcn,est_dist,Intrinsic);% apply calibration if isempty(GeometryCalib) return @@ -393,5 +405,5 @@ %------------------------------------------------------------------------ % --- activate calibration and store parameters in ouputfile . -function [GeometryCalib,ind_max,ind_removed]=calibrate(Coord,CalibFcn,Intrinsic) +function [GeometryCalib,ind_max,ind_removed]=calibrate(Coord,CalibFcn,est_dist,Intrinsic) %------------------------------------------------------------------------ @@ -403,5 +415,5 @@ CoordFlip(:,3)=Coord(:,3);% the calibration function assume a z ccordinate along the camera view, opposite to ours if ~isempty(Coord) - GeometryCalib=feval(CalibFcn,CoordFlip,Intrinsic); + GeometryCalib=feval(CalibFcn,CoordFlip,est_dist,Intrinsic); else msgbox_uvmat('ERROR','No calibration points, abort') @@ -433,5 +445,5 @@ CoordFlip=Coord; CoordFlip(:,3)=Coord(:,3);% the calibration function assume a z ccordinate along the camera view, opposite to ours - GeometryCalib=feval(CalibFcn,CoordFlip,Intrinsic); + GeometryCalib=feval(CalibFcn,CoordFlip,est_dist,Intrinsic); X=Coord(:,1); Y=Coord(:,2); @@ -448,5 +460,5 @@ %------------------------------------------------------------------------ % --- determine the parameters for a calibration by an affine function (rescaling and offset, no rotation) -function GeometryCalib=calib_rescale(Coord,Intrinsic) +function GeometryCalib=calib_rescale(Coord,est_dist,Intrinsic) %------------------------------------------------------------------------ X=Coord(:,1); @@ -464,5 +476,5 @@ %------------------------------------------------------------------------ % --- determine the parameters for a calibration by a linear transform matrix (rescale and rotation) -function GeometryCalib=calib_linear(Coord,Intrinsic) +function GeometryCalib=calib_linear(Coord,est_dist,Intrinsic) %------------------------------------------------------------------------ X=Coord(:,1); @@ -494,5 +506,5 @@ % --- determine the tsai parameters for a view normal to the grid plane % NOT USED -function GeometryCalib=calib_normal(Coord,Intrinsic) +function GeometryCalib=calib_normal(Coord,est_dist,Intrinsic) %------------------------------------------------------------------------ Calib.fx=str2num(get(handles.fx,'String')); @@ -555,75 +567,75 @@ alpha=calib_param(5); GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1]; - -%------------------------------------------------------------------------ -function GeometryCalib=calib_3D_linear(Coord,Intrinsic) -%------------------------------------------------------------------------ -coord_files=Intrinsic.coord_files; -if ischar(Intrinsic.coord_files) - coord_files={coord_files}; -end -if isempty(coord_files{1}) || isequal(coord_files,{''}) - coord_files={}; -end -%retrieve the calibration points stored in the files listed in the popup list ListCoordFiles -x_1=Coord(:,4:5)';%px coordinates of the ref points - -nx=Intrinsic.Npx; -ny=Intrinsic.Npy; -x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates -X_1=Coord(:,1:3)';%phys coordinates of the ref points -n_ima=numel(coord_files)+1; -if ~isempty(coord_files) - msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used']) - for ifile=1:numel(coord_files) - t=xmltree(coord_files{ifile}); - s=convert(t);%convert to matlab structure - if isfield(s,'GeometryCalib') - if isfield(s.GeometryCalib,'SourceCalib') - if isfield(s.GeometryCalib.SourceCalib,'PointCoord') - PointCoord=s.GeometryCalib.SourceCalib.PointCoord; - Coord_file=zeros(length(PointCoord),5);%default - for i=1:length(PointCoord) - line=str2num(PointCoord{i}); - Coord_file(i,4:5)=line(4:5);%px x - Coord_file(i,1:3)=line(1:3);%phys x - end - eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']); - eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]); - eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']); - end - end - end - end -end -n_ima=numel(coord_files)+1; -est_dist=[0;0;0;0;0]; -est_aspect_ratio=0; -est_fc=[1;1]; -center_optim=0; -path_uvmat=which('uvmat');% check the path detected for source file uvmat -path_UVMAT=fileparts(path_uvmat); %path to UVMAT -run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));% apply fct 'toolbox_calib/go_calib_optim' -if exist('Rc_1','var') - GeometryCalib.CalibrationType='3D_extrinsic'; - GeometryCalib.fx_fy=fc'; - GeometryCalib.Cx_Cy=cc'; - GeometryCalib.kc=kc(1); - GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function - GeometryCalib.Tx_Ty_Tz=Tc_1'; - GeometryCalib.R=Rc_1; - GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate - GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate - GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate - GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees - GeometryCalib.ErrorRMS=[]; - GeometryCalib.ErrorMax=[]; -else - msgbox_uvmat('ERROR',['calibration function ' fullfile('toolbox_calib','go_calib_optim') ' did not converge: use multiple views or option 3D_extrinsic']) - GeometryCalib=[]; -end - -%------------------------------------------------------------------------ -function GeometryCalib=calib_3D_quadr(Coord,Intrinsic) +% +% %------------------------------------------------------------------------ +% function GeometryCalib=calib_3D_linear(Coord,Intrinsic) +% %------------------------------------------------------------------------ +% coord_files=Intrinsic.coord_files; +% if ischar(Intrinsic.coord_files) +% coord_files={coord_files}; +% end +% if isempty(coord_files{1}) || isequal(coord_files,{''}) +% coord_files={}; +% end +% %retrieve the calibration points stored in the files listed in the popup list ListCoordFiles +% x_1=Coord(:,4:5)';%px coordinates of the ref points +% +% nx=Intrinsic.Npx; +% ny=Intrinsic.Npy; +% x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates +% X_1=Coord(:,1:3)';%phys coordinates of the ref points +% n_ima=numel(coord_files)+1; +% if ~isempty(coord_files) +% msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used']) +% for ifile=1:numel(coord_files) +% t=xmltree(coord_files{ifile}); +% s=convert(t);%convert to matlab structure +% if isfield(s,'GeometryCalib') +% if isfield(s.GeometryCalib,'SourceCalib') +% if isfield(s.GeometryCalib.SourceCalib,'PointCoord') +% PointCoord=s.GeometryCalib.SourceCalib.PointCoord; +% Coord_file=zeros(length(PointCoord),5);%default +% for i=1:length(PointCoord) +% line=str2num(PointCoord{i}); +% Coord_file(i,4:5)=line(4:5);%px x +% Coord_file(i,1:3)=line(1:3);%phys x +% end +% eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']); +% eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]); +% eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']); +% end +% end +% end +% end +% end +% n_ima=numel(coord_files)+1; +% est_dist=[0;0;0;0;0]; +% est_aspect_ratio=0; +% est_fc=[1;1]; +% center_optim=0; +% path_uvmat=which('uvmat');% check the path detected for source file uvmat +% path_UVMAT=fileparts(path_uvmat); %path to UVMAT +% run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));% apply fct 'toolbox_calib/go_calib_optim' +% if exist('Rc_1','var') +% GeometryCalib.CalibrationType='3D_extrinsic'; +% GeometryCalib.fx_fy=fc'; +% GeometryCalib.Cx_Cy=cc'; +% GeometryCalib.kc=kc(1); +% GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function +% GeometryCalib.Tx_Ty_Tz=Tc_1'; +% GeometryCalib.R=Rc_1; +% GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate +% GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate +% GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate +% GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees +% GeometryCalib.ErrorRMS=[]; +% GeometryCalib.ErrorMax=[]; +% else +% msgbox_uvmat('ERROR',['calibration function ' fullfile('toolbox_calib','go_calib_optim') ' did not converge: use multiple views or option 3D_extrinsic']) +% GeometryCalib=[]; +% end + +%------------------------------------------------------------------------ +function GeometryCalib=calib_3D(Coord,est_dist,Intrinsic) %------------------------------------------------------------------ @@ -678,5 +690,5 @@ end n_ima=numel(coord_files)+1; -est_dist=[1;0;0;0;0]; +%est_dist=[1;0;0;0;0]; est_aspect_ratio=1; center_optim=0; @@ -687,5 +699,8 @@ GeometryCalib.fx_fy=fc'; GeometryCalib.Cx_Cy=cc'; - GeometryCalib.kc=kc(1); + ind_kc=find(kc); + if ~isempty(ind_kc) + GeometryCalib.kc=(kc(ind_kc))';%include cases quadr (one value kc(1)) and order_4 (2 values for kc) + end GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=Tc_1'; @@ -703,5 +718,5 @@ %------------------------------------------------------------------------ -function GeometryCalib=calib_3D_extrinsic(Coord, Intrinsic) +function GeometryCalib=calib_3D_extrinsic(Coord,est_dist, Intrinsic) %------------------------------------------------------------------ path_uvmat=which('geometry_calib');% check the path detected for source file uvmat @@ -747,8 +762,10 @@ addpath(fct_path) GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%reverse Cx_Cy(2) for calibration (inversion of px ordinate) +kc=zeros(1,5); +kc(1:numel(GeometryCalib.kc))=GeometryCalib.kc; [omc,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,... - (GeometryCalib.fx_fy)',GeometryCalib.Cx_Cy',[GeometryCalib.kc 0 0 0 0]); + (GeometryCalib.fx_fy)',GeometryCalib.Cx_Cy',kc); rmpath(fct_path); -GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function +GeometryCalib.CoordUnit='';% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=Tc1'; %inversion of z axis @@ -1086,5 +1103,5 @@ Tmod(:,2)=(TIndex(:,2)-1)*Dy+Rangy(1); %Tmod=T(:,(1:2))+Delta;% 'phys' coordinates of the detected points -[Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));% image coordinates of the detected points +[Xpx,Ypx]=px_XYZ(GeometryCalib,[],Tmod(:,1),Tmod(:,2));% image coordinates of the detected points Coord=[T Xpx Ypx zeros(size(T,1),1)]; set(handles.ListCoord,'Data',Coord) @@ -1326,5 +1343,6 @@ set(handles.Cx,'String',num2str(Cx_Cy(1),'%1.1f')) set(handles.Cy,'String',num2str(Cx_Cy(2),'%1.1f')) -set(handles.kc,'String',num2str(kc,'%1.4f')) +%set(handles.kc,'String',num2str(kc,'%1.4f')) +set(handles.kc,'String',num2str(kc,4)) %------------------------------------------------------------------------ Index: trunk/src/phys_XYZ.m =================================================================== --- trunk/src/phys_XYZ.m (revision 1114) +++ trunk/src/phys_XYZ.m (revision 1115) @@ -74,12 +74,13 @@ Calib.Cx_Cy=[0 0]; end -if ~isfield(Calib,'kc') - Calib.kc=0; +kc=[0 0]; +if isfield(Calib,'kc') + kc(1:numel(Calib.kc))=Calib.kc; end +% if ~isfield(Calib,'kc2') +% Calib.kc2=0; +% end if isfield(Calib,'R') R=(Calib.R)'; -% R(3)=-R(3); -% R(6)=-R(6); -% R(9)=-R(9); c=Z0virt; cvirt=Z0virt; @@ -126,5 +127,5 @@ Xd=(X-Calib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates Yd=(Y-Calib.Cx_Cy(2))/Calib.fx_fy(2); - dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd);% distortion factor, first approximation Xu,Yu=Xd,Yd + dist_fact=1+kc(1)*(Xd.*Xd+Yd.*Yd);% distortion factor, first approximation Xu,Yu=Xd,Yd test=0; niter=0; @@ -133,5 +134,6 @@ Xu=Xd./dist_fact;%undistorted sensor coordinates, second iteration Yu=Yd./dist_fact; - dist_fact=1+Calib.kc*(Xu.*Xu+Yu.*Yu);% distortion factor,next approximation + R2=Xu.*Xu+Yu.*Yu; + dist_fact=1+kc(1)*R2+kc(2)*R2.*R2;% distortion factor,next approximation test=max(max(abs(dist_fact-dist_fact_old)))<0.00001; % reducing the relative error to 10^-5 forthe inversion of the quadraticcorrection niter=niter+1; Index: trunk/src/proj_field.m =================================================================== --- trunk/src/proj_field.m (revision 1114) +++ trunk/src/proj_field.m (revision 1115) @@ -1095,5 +1095,5 @@ end else - ProjData.ProjObjectAngle=FieldData.PlaneAngle; + ProjData.ProjObjectAngle=FieldData.PlaneAngle;ProjMode end end @@ -1125,10 +1125,10 @@ icell_grid=[];% field cell index which defines the grid icell_scattered=[];% field cell index which defines fields with scattered coordinates -for icell=1:numel(CellInfo) - if strcmp(ObjectData.ProjMode,'interp_lin')&& ~strcmp(ProjMode{icell},'interp_tps') - errormsg='ProjMode interp_tps needed '; - return - end -end +% for icell=1:numel(CellInfo) +% if strcmp(ObjectData.ProjMode,'interp_lin')&& ~strcmp(ProjMode{icell},'interp_tps') +% errormsg='ProjMode interp_tps needed '; +% return +% end +% end if strcmp(ObjectData.ProjMode,'projection') %% case of a grid requested by the input field Index: trunk/src/px_XYZ.m =================================================================== --- trunk/src/px_XYZ.m (revision 1114) +++ trunk/src/px_XYZ.m (revision 1115) @@ -69,6 +69,9 @@ if ~isfield(Calib,'kc') r2=1; %no quadratic distortion + elseif numel(Calib.kc)==1 + r2=1+Calib.kc*(Xu.*Xu+Yu.*Yu); else - r2=1+Calib.kc*(Xu.*Xu+Yu.*Yu); + R2=Xu.*Xu+Yu.*Yu; + r2=1+Calib.kc(1)*R2+Calib.kc(2)*R2.*R2; end Index: trunk/src/series/civ_series.m =================================================================== --- trunk/src/series/civ_series.m (revision 1114) +++ trunk/src/series/civ_series.m (revision 1115) @@ -628,5 +628,9 @@ ind_good=1:numel(Data.Civ1_X); end - + if isempty(ind_good) + disp_uvmat('ERROR','all vectors of civ1 are bad, check input parameters' ,checkrun) + return + end + % perform Patch calculation using the UVMAT fct 'filter_tps' [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=... @@ -919,4 +923,8 @@ else ind_good=1:numel(Data.Civ2_X); + end + if isempty(ind_good) + disp_uvmat('ERROR','all vectors of civ2 are bad, check input parameters' ,checkrun) + return end [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=... Index: trunk/src/series/merge_proj.m =================================================================== --- trunk/src/series/merge_proj.m (revision 1114) +++ trunk/src/series/merge_proj.m (revision 1115) @@ -86,5 +86,7 @@ return end - +if 0==1 + phys; % used to include phys when compiling is done +end %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%% ParamOut=[]; %default output Index: trunk/src/series/turb_stat.m =================================================================== --- trunk/src/series/turb_stat.m (revision 1114) +++ trunk/src/series/turb_stat.m (revision 1115) @@ -204,5 +204,5 @@ % for i_slice=1:Param.IndexRange.NbSlice % i_slice - ind_first=Param.IndexRange.first_i + ind_first=Param.IndexRange.first_i; for index_i=ind_first:Param.IndexRange.NbSlice:Param.IndexRange.last_i if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue') @@ -210,5 +210,5 @@ break end - for index_j=Param.IndexRange.first_j:Param.IndexRange.last_j + for index_j=first_j:last_j InputFile=fullfile_uvmat(RootPath{1},SubDir{1},RootFile{1},FileExt{1},NomType{1},index_i,index_i,index_j,index_j); [Field,tild,errormsg] = read_field(InputFile,FileType{iview},InputFields{iview}); @@ -217,5 +217,5 @@ %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% - if index_i==ind_first && index_j==Param.IndexRange.first_j %initiate the output data structure in the first field + if index_i==ind_first && index_j==first_j %initiate the output data structure in the first field [CellInfo,NbDim,errormsg]=find_field_cells(Field); YName='coord_y';%default