Changeset 1184 for trunk/src/series/bed_scan.m

Timestamp:

Nov 6, 2025, 7:03:02 PM (4 weeks ago)

Author:

sommeria

Message:

bed-scan updated and many updates

File:

• trunk/src/series/bed_scan.m (modified) (7 diffs)

trunk/src/series/bed_scan.m

@@ -61 +61 @@
 %=======================================================================
 
-function ParamOut=bed_scan (Param)
+function ParamOut=bed_scan(Param)
 
 %% set the input elements needed on the GUI series when the action is selected in the menu ActionName or InputTable refreshed
@@ -67 +67 @@
     ParamOut.NbViewMax=1;% max nbre of input file series (default , no limitation)
     ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
-    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
-    ParamOut.NbSlice=1; %nbre of slices ('off' by default)
+    ParamOut.WholeIndexRange='on';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
+    ParamOut.NbSlice='off'; %nbre of slices ('off' by default)
     ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
     ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
@@ -75 +75 @@
     ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
     ParamOut.OutputDirExt='.bed';%set the output dir extension
-    ParamOut.OutputFileMode='NbSlice';% ='=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice
-    %check the type of the existence and type of the first input file:
-    Param.IndexRange.last_i=Param.IndexRange.first_i;%keep only the first index in the series
-    if isfield(Param.IndexRange,'first_j')
-    Param.IndexRange.last_j=Param.IndexRange.first_j;
-    end
-    filecell=get_file_series(Param);
-    if ~exist(filecell{1,1},'file')
-        msgbox_uvmat('WARNING','the first input file does not exist')
-    else
-        FileInfo=get_file_info(filecell{1,1});
-        FileType=FileInfo.FileType;
-        if isempty(find(strcmp(FileType,{'image','multimage','mmreader','video'})))% =1 for images
-            msgbox_uvmat('ERROR',['bad input file type for ' mfilename ': an image is needed'])
-        end
-    end
-return
+    ParamOut.OutputFileMode='NbInput_i';% ='=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice   
+    return
 end
 
@@ -111 +96 @@
 %% root input file names and nomenclature type (cell arrays with one element)
 RootPath=Param.InputTable{1,1};
-
+SubDir=Param.InputTable{1,2};
+RootFile=Param.InputTable{1,3};
+NomType=Param.InputTable{1,4};
+FileExt=Param.InputTable{1,5};
+i_series=Param.IndexRange.first_i:Param.IndexRange.incr_i:Param.IndexRange.last_i;
+j_series=Param.IndexRange.first_j:Param.IndexRange.incr_j:Param.IndexRange.last_j;
+nbfield_i=numel(i_series);
+nbfield_j=numel(j_series);
 
 %% directory for output files
@@ -117 +109 @@
 
 %% get the set of input file names (cell array filecell), and file indices
-[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
+%[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
 % filecell{iview,fileindex}: cell array representing the list of file names
 %        iview: line in the table corresponding to a given file series
-%        fileindex: file index within  the file series, 
-% i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j 
+%        fileindex: file index within  the file series,
+% i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j
 % i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
-nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
-nbfield_i=size(i1_series{1},2); %nb of fields for the i index
-
-%% set of y positions  
-ycalib=[-51 -1 49];% calibration planes
-y_scan=-51+(100/400)*(i1_series{1}-1);% transverse position given by the translating system: first view at y=-51, view 400 at y=+49
+%nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
+% nbfield_i=size(i1_series{1},2); %nb of fields for the i index
+% 
+% nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
+CheckVirtual=false;
+if isfield(Param,'FileSeries')% virtual file indexing used (e.g. multitif images)
+    CheckVirtual=true;
+end
+
+%%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%fullfile(
+% EDIT FROM HERE
+%% load the initial scan
+[RootRoot,CamName]=fileparts(RootPath);
+RootRoot=fileparts(RootRoot);
+CalibFolder=fullfile(RootRoot,'EXP_INIT',CamName);
+File_init=fullfile(CalibFolder,'images.png.bed','Z_init.nc');
+Data_init=nc2struct(File_init);
+%% set of y positions
+
+nb_scan=400;% nbre of planes for a scan
+if nbfield_j<400
+    nb_scan=nbfield_j;
+end
+%ycalib=[-51 -1 49];% calibration planes
+y_scan=-51+0.25*(1:nb_scan);% transverse position given by the translating system: first view at y=-51, view 400 at y=+49
+coord_x=0.25:0.25:450;%% coord x in phys coordinates for final projection
+
 Mfiltre=ones(2,10)/20;%filter matrix for imnages
 
 %% calibration data and timing: read the ImaDoc files
-XmlData_A=xml2struct(fullfile(RootPath,'planeA.xml'));
-XmlData_B=xml2struct(fullfile(RootPath,'planeB.xml'));
-XmlData_C=xml2struct(fullfile(RootPath,'planeC.xml'));
-ycalib=[-51 -1 49];% the three y positions for calibration
-fx(1)=XmlData_A.GeometryCalib.fx_fy(1);
+
+XmlData_A=xml2struct(fullfile(CalibFolder,'planeA.xml'));
+XmlData_B=xml2struct(fullfile(CalibFolder,'planeB.xml'));
+XmlData_C=xml2struct(fullfile(CalibFolder,'planeC.xml'));
+ycalib=[-51 -1 49];% the three y positions for calibration=
+fx(1)=XmlData_C.GeometryCalib.fx_fy(1);
 fx(2)=XmlData_B.GeometryCalib.fx_fy(1);
-fx(3)=XmlData_C.GeometryCalib.fx_fy(1);
-fy(1)=XmlData_A.GeometryCalib.fx_fy(2);
+fx(3)=XmlData_A.GeometryCalib.fx_fy(1);
+fy(1)=XmlData_C.GeometryCalib.fx_fy(2);
 fy(2)=XmlData_B.GeometryCalib.fx_fy(2);
-fy(3)=XmlData_C.GeometryCalib.fx_fy(2);
-Tx(1)=XmlData_A.GeometryCalib.Tx_Ty_Tz(1);
+fy(3)=XmlData_A.GeometryCalib.fx_fy(2);
+Tx(1)=XmlData_C.GeometryCalib.Tx_Ty_Tz(1);
 Tx(2)=XmlData_B.GeometryCalib.Tx_Ty_Tz(1);
-Tx(3)=XmlData_C.GeometryCalib.Tx_Ty_Tz(1);
-Ty(1)=XmlData_A.GeometryCalib.Tx_Ty_Tz(2);
+Tx(3)=XmlData_A.GeometryCalib.Tx_Ty_Tz(1);
+Ty(1)=XmlData_C.GeometryCalib.Tx_Ty_Tz(2);
 Ty(2)=XmlData_B.GeometryCalib.Tx_Ty_Tz(2);
-Ty(3)=XmlData_C.GeometryCalib.Tx_Ty_Tz(2);
-R11(1)=XmlData_A.GeometryCalib.R(1,1);
+Ty(3)=XmlData_A.GeometryCalib.Tx_Ty_Tz(2);
+R11(1)=XmlData_C.GeometryCalib.R(1,1);
 R11(2)=XmlData_B.GeometryCalib.R(1,1);
-R11(3)=XmlData_C.GeometryCalib.R(1,1);
-R12(1)=XmlData_A.GeometryCalib.R(1,2);
+R11(3)=XmlData_A.GeometryCalib.R(1,1);
+R12(1)=XmlData_C.GeometryCalib.R(1,2);
 R12(2)=XmlData_B.GeometryCalib.R(1,2);
-R12(3)=XmlData_C.GeometryCalib.R(1,2);
-R21(1)=XmlData_A.GeometryCalib.R(2,1);
+R12(3)=XmlData_A.GeometryCalib.R(1,2);
+R21(1)=XmlData_C.GeometryCalib.R(2,1);
 R21(2)=XmlData_B.GeometryCalib.R(2,1);
-R21(3)=XmlData_C.GeometryCalib.R(2,1);
-R22(1)=XmlData_A.GeometryCalib.R(2,2);
+R21(3)=XmlData_A.GeometryCalib.R(2,1);
+R22(1)=XmlData_C.GeometryCalib.R(2,2);
 R22(2)=XmlData_B.GeometryCalib.R(2,2);
-R22(3)=XmlData_C.GeometryCalib.R(2,2);
-pfx=polyfit(ycalib,fx,1);%get the linear interpolation of each parameter of the three calibrations
+R22(3)=XmlData_A.GeometryCalib.R(2,2);
+pfx=polyfit(ycalib,fx,1);%get thfield_ie linear interpolation of each parameter of the three calibrations
 pfy=polyfit(ycalib,fy,1);
 pTx=polyfit(ycalib,Tx,1);
@@ -169 +183 @@
 p22=polyfit(ycalib,R22,1);
 %get the calibration parameters at each position y by interpolation of the 3 calibration parameters
-for img=1:nbfield_i
+for img=1:nb_scan
     Calib(img).fx_fy(1)=pfx(1)*y_scan(img)+pfx(2);
     Calib(img).fx_fy(2)=pfy(1)*y_scan(img)+pfy(2);
@@ -184 +198 @@
 end
 
-%% check coincdence in time for several input file series
-%not relevant for this function
-
-%% coordinate transform or other user defined transform
-%not relevant for this function
-
-%%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
- % EDIT FROM HERE
+
+
 
 %% Load the init bed scan
 tic
-nb_scan=10;
-for img=1:nb_scan
-     img
-    a=flipud(imread(filecell{1,img}));%image of the initial bed
-    if img==1
-        x=1:size(a,2);%image absissa in pixel coordinates
+%filecell=reshape(filecell,nbfield_j,nbfield_i)
+% main loop
+for ifield=1:nbfield_i
+    for img=1:nbfield_j% loop on positions
+        if CheckVirtual
+            [FileName,FrameIndex]=index2filename(Param.FileSeries,i_series(ifield),j_series(img),nbfield_j);
+            InputFile=fullfile(RootPath,SubDir,FileName);
+        else
+            InputFile=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,FileIndex,ifield,[],img);
+            FrameIndex=1;
+        end
+        a=flipud(read_image(InputFile,'multimage',[],FrameIndex));
+        %a=flipud(imread(InputFile));%image of the initial bed  [X_b_new(img,:),Z_b_new(img,:)]=phys_XYZ(Calib(img),[],x,Z_sb(img,:))
+        if img==1
+            [nby,nbx]=size(a);
+            x_ima=1:nbx;%image absissa in pixel coordinates
+            X_phys=zeros(nb_scan,nbx);
+            Z_phys=zeros(nb_scan,nbx);
+        end
+        % filtering
+            a=filter2(Mfiltre,a);%smoothed image
+        amean=mean(a,2);
+     [~,ind_max]=max(amean);% get the max of the image averaged along x, to restrict the search region
+    ind_range=max(1,ind_max-30):min(nby,ind_max+30);% search band to find the line
+    z_ima=get_max(a(ind_range,:))+ind_range(1)-1;% get the max in the search band and shift to express it in indices of the original image
+
+ %       [~,iy]=max(a);% find the max along the first coordinate     Z_s_new=zeros(nb_scan,size(Z_s,2)); y, max values imax and the corresponding  y index iy along the first coordinate y
+        z_ima=smooth(z_ima,20,'rloess');%smooth Z, the image index of max luminosity (dependning on x)
+        [X_phys(img,:),Z_phys(img,:)]=phys_XYZ(Calib(img),[],x_ima,z_ima');
     end
-    % filtering
-    a=filter2(Mfiltre,a);%smoothed image
-    [imax,iy]=max(a);% find the max along the first coordinate y, max values imax and the corresponding  y index iy along the first coordinate y
-    Z_s(img,:)=smooth(iy,40,'rloess');%smooth Z, the image index of max luminosity (dependning on x)
-    Yima=y_scan(img)*ones(size(x));%positions Y transformed into a vector
-    X_new(img,:)=phys_XYZ(Calib(img),x,Yima,1);
-    %X_new(:,img)=phys_scan(x,y(img));
-end
-
-toc
-
-[X,Y]=meshgrid(x,y_scan);
-
-
-%% Load the current bed scan
-for img=1:nb_scan
-   b=flipud(imread(filecell{2,img}));%image of the current bed
-    b=filter2(Mfiltre,b); % filtering
-    [imaxb,iyb]=max(b);
-    Z_sb(img,:)=smooth(iyb,20,'rloess');
-end
-
-%% bed change
-dZ=Z_s-Z_sb;% displacement between current position and initial
-dZ_new=zeros(nb_scan,size(dZ,2));
-for img=1:nb_scan  
-    Yima=y_scan(img)*ones(1,size(dZ,2));
-    [~,dZ_new(img,:)]=phys_XYZ(Calib(img),dZ(img,:),Yima,1);
-   % dZ_new(:,img)=phys_scanz(dZ(:,img),y(img));
-end
-
-
-%% PLOTS
-coord_x=X_new(end,1):0.1:X_new(end,end);
-[Y_m,X_m]=meshgrid(y_scan,coord_x);
-%Y_new=Y';
-dZ_mesh=griddata(X_new,Y,dZ_new,X_m,Y_m);
-
-if checkrun
-    figure(1)
-    hold on
-    plot(x,Z_s+700)
-    % xlim([0 4096])
-    % ylim([0 3000])
-    
-    figure(2)
-    hold on
-    plot(x,Z_sb+700)
-    xlim([0 4096])
-    ylim([0 3000])
-    
-    figure(3)
-    surfc(X_m,Y_m,dZ_mesh)
-    shading interp;
-    colorbar;
-    caxis([0 3]);
-    
-    figure
-    pcolor(X_m,Y_m,dZ_mesh);
-    colormap;
-    set(gca,'Xdir','reverse');
-    caxis([0 3]);
-    shading flat
-    hold on
-    colorbar
-    title('Dz')
-end
-
-%save(fullfile(DirOut,'18OS_f.mat'),'dZ','dZ_new','X','Y','Z_s','Z_sb','y')
-
-% save netcdf
-Data.ListVarName={'coord_x','coord_y','dZ'};
-Data.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'}};
-Data.VarAttribute{1}.Role='coord_x';
-Data.VarAttribute{1}.unit='cm';
-Data.VarAttribute{2}.Role='coord_y';
-Data.VarAttribute{2}.unit='cm';
-Data.VarAttribute{3}.Role='scalar';
-Data.VarAttribute{3}.unit='cm';
-Data.coord_x=[coord_x(1) coord_x(end)];
-Data.coord_y=[y(1) y(end)];
-Data.dZ=dZ_mesh';
-struct2nc(fullfile(DirOut,'dZ.nc'),Data)
-
-%% gives the physical position x from the image position X and the physical position y of the laser plane
-function F=phys_scan(X,y)
-% linear fct of X whose coefficient depend on y in a quadratic way
-F=(9.4*10^(-7)*y.^2-3.09*10^(-4)*y+0.07).*X +(-0.001023*y.^2+0.469*y+186.9);
-
-%% gives the physical position z from the image position Z and the physical position y of the laser plane
-function Fz=phys_scanz(Z,y)
-% scale factor applied to Z depending on the physical position y of the laser plane
-Fz=(-1.4587*10^(-5)*y.^2 + 0.001072*y+0.0833).*Z; %+(-2.1*10^(-6)*x.^2+5.1*10^(-4)*x+0.0735).*Z;
-
- 
+    disp(['last file of ' num2str(ifield)])
+     disp(FileName)
+     disp(FrameIndex)
+
+    %% interpolate on a regular grid
+    %coord_x=X_phys(end,1):0.1:X_phys(end,end);%% coord x in phys coordinates based in the last view plane (the last)
+    [X_m,Y_m]=meshgrid(coord_x,y_scan);
+    Y=y_scan'*ones(1,nbx);%initialisation of X, Y final topography map
+
+    Data.Z=griddata(X_phys,Y,Z_phys,X_m,Y_m);% dZ interpolated on the regular ph1ys grid X_m,Y_m
+    size(Data.Z)
+    size(Data_init.Z_init)
+    Data.dZ=Data.Z-Data_init.Z_init;
+
+    toc
+
+    % save netcdf
+    Data.ListVarName={'coord_x','y_scan','Z','dZ'};
+    Data.VarDimName={'coord_x','y_scan',{'y_scan','coord_x'},{'y_scan','coord_x'}};
+    Data.VarAttribute{1}.Role='coord_x';
+    Data.VarAttribute{1}.unit='cm';
+    Data.VarAttribute{2}.Role='coord_y';
+    Data.VarAttribute{2}.unit='cm';
+    Data.VarAttribute{3}.Role='scalar';
+    Data.VarAttribute{3}.unit='cm';
+        Data.VarAttribute{4}.Role='scalar';
+    Data.VarAttribute{4}.unit='cm';
+    Data.coord_x=coord_x;
+    Data.y_scan=y_scan;
+    struct2nc(fullfile(DirOut,['dZ_' num2str(ifield) '.nc']),Data)
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function iy=get_max(a)% get the max with sub picel resolution
+[a_max,iy]=max(a);
+[Nby,Nbx]=size(a);
+for ind_x=1:Nbx
+    if iy(ind_x)>1 && iy(ind_x)<Nby
+        a_plus=a(iy(ind_x)+1,ind_x);
+        a_min=a(iy(ind_x)-1,ind_x);
+        denom=2*a_max(ind_x)-a_plus-a_min;
+        if denom >0
+            iy(ind_x)=iy(ind_x)+0.5*(a_plus-a_min)/denom;%adjust the position of the max with a quadratic fit of the three points around the max
+        end
+    end
+end