Changeset 214 for trunk


Ignore:
Timestamp:
Mar 6, 2011, 10:22:33 PM (13 years ago)
Author:
sommeria
Message:

relabel_i_j and sub_background improved. Possibility of using levels

Location:
trunk/src/series
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/series/relabel_i_j.m

    r184 r214  
    1 % relabel_i_j: relabel an image series with two indices, according to the time matrix given by ImaDoc
    2 % (specific to RDvision system)
     1% relabel_i_j: relabel an image series with two indices, and correct errors from the RDvision transfer program
    32%----------------------------------------------------------------------
    43function GUI_input=relabel_i_j(num_i1,num_i2,num_j1,num_j2,Series)
    54%requests for the visibility of input windows in the GUI series  (activated directly by the selection in the menu ACTION)
    6 if ~exist('num_i1','var')
    7     GUI_input={};
    8     return %exit the function
    9 end
     5
     6GUI_input={};
    107
    118%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%enable waitbar
    12 hseries=guidata(Series.hseries);%handles of the GUI series
     9hGUI=findobj(allchild(0),'name','series');
     10hseries=guidata(hGUI);%handles of the GUI series
    1311WaitbarPos=get(hseries.waitbar_frame,'Position');
    1412%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    1715display('RDvision system')
    1816first_label=0 %image numbers start from 0
    19 errorfactor=2 %correct a factor of 2 in NbDk+1
    20 
    21 %% create dir of the new images
    22 basename=fullfile(Series.RootPath,Series.RootFile) ;
    23 [dir_images,namebase]=fileparts(basename);
    24 [path,subdir_ima]=fileparts(dir_images);
    25 newdir=fullfile(path,[subdir_ima '_ij']);
    26 mkdir(newdir);
    27 [xx,msg2] = fileattrib(newdir,'+w','g'); %yield writing access (+w) to user group (g)
    28 if ~strcmp(msg2,'')
    29     msgbox_uvmat('ERROR',['pb of permission for ' subdir_ima ': ' msg2])%error message for directory creation
     17%errorfactor=1 %correct a factor of 2 in NbDk+1
     18
     19%% read imadoc
     20RootPath=get(hseries.RootPath,'String');
     21RootFile=get(hseries.RootFile,'String');
     22basename=fullfile(RootPath{1},RootFile{1});
     23SeriesData=get(hGUI,'UserData');
     24if ~strcmp(SeriesData.NomType,'_000001')
     25    msgbox_uvmat('WARNING','the input is not a file from RDvision: this function relabel_i_j has no action')%error message for directory creation
    3026    return
    31 end
    32 display(['relabelled images in the directory ' newdir])
    33 basename_new=fullfile(newdir,namebase);
    34 
    35 %% read imadoc
    36 [XmlData,warntext]=imadoc2struct([basename '.xml']);
    37 nbfield1=size(XmlData.Time,1)/errorfactor;
     27else
     28    msgbox_uvmat('CONFIRMATION','this function will relabel the file series from RDvision and correct the xml file')%error message for directory creation
     29end
     30[XmlData,warntext]=imadoc2struct([basename '.xml'])% read the xml file appended to the present function (containing bug corrections)
     31if ~isempty(warntext)
     32    msgbox_uvmat('ERROR',warntext)%error message for xml file reading
     33end
     34nbfield1=size(XmlData.Time,1);
    3835nbfield2=size(XmlData.Time,2);
     36
    3937set(hseries.first_i,'String',num2str(first_label))% display the first image in the process
    4038set(hseries.last_i,'String',num2str(nbfield1*nbfield2-1+first_label))% display the last image in the process
     39set(hseries.nb_field,'String',{num2str(nbfield1*nbfield2-1+first_label)})% display the total nbre of images
     40
     41%% stop program ther when it is selected in the menu (no run action)
     42if ~exist('num_i1','var')
     43    return
     44end
     45
     46answer=msgbox_uvmat('CONFIRMATION',[num2str(nbfield1) ' bursts containing ' num2str(nbfield2) ' images each'])%error message for directory creation
    4147
    4248%% apply the image rescaling function 'level' (avoid bright particles)
    43 answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m');
    44 test_level=isequal(answer,'Yes');
     49% answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m');
     50% test_level=isequal(answer,'Yes');
    4551
    4652%% copy and adapt the xml file
    4753if exist([basename '.xml'],'file')
    48     copyfile([basename '.xml'],[basename_new '.xml']);% copy the .civ file
    49     t=xmltree([basename_new '.xml']);
     54    try
     55    copyfile([basename '.xml'],[basename '.xml~']);% backup the xml file
     56    catch
     57    errormsg=lasterr
     58            msgbox_uvmat('ERROR',errormsg);
     59            return
     60    end
     61    t=xmltree([basename '.xml']);
    5062   
    5163    %update information on the first image name in the series
     
    5567    end   
    5668    uid_ImageName=find(t,'ImaDoc/Heading/ImageName');
    57     ImageName=name_generator(basename_new,1,1,'.png','_i_j');
     69    ImageName=name_generator(basename,1,1,'.png','_i_j');
    5870    [pth,ImageName]=fileparts(ImageName);
    5971    ImageName=[ImageName '.png']
     
    6375    uid_value=children(t,uid_ImageName);
    6476    if isempty(uid_value)
    65         t=add(t,uid_ImageName,'chardata',ImageName)%indicate  name of the first image, with ;png extension
     77        t=add(t,uid_ImageName,'chardata',ImageName);%indicate  name of the first image, with ;png extension
    6678    else
    67         t=set(t,uid_value(1),'value',ImageName)%indicate  name of the first image, with ;png extension
    68     end 
    69 
    70 %     %add information about image transform
    71 %     [t,new_uid]=add(t,1,'element','ImageTransform');
    72 %     [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
    73 %     [t]=add(t,NameFunction_uid,'chardata','sub_background');     
    74 %     [t,NbSlice_uid]=add(t,new_uid,'element','NbSlice');
    75 %     [t]=add(t,new_uid,'chardata',num2str(nbslice_i));
    76 %     [t,NbSlidingImages_uid]=add(t,new_uid,'element','NbSlidingImages');
    77 %     [t]=add(t,NbSlidingImages_uid,'chardata',num2str(nbaver));
    78 %     [t,LuminosityRank_uid]=add(t,new_uid,'element','RankBackground');
    79 %     [t]=add(t,LuminosityRank_uid,'chardata',num2str(rank));% luminosity rank almong the nbaver sliding images
    80     save(t,[basename_new '.xml'])
     79        t=set(t,uid_value(1),'value',ImageName);%indicate  name of the first image, with ;png extension
     80    end
     81   
     82    %%%% correction RDvision %%%%
     83    uid_NbDtj=find(t,'ImaDoc/Camera/BurstTiming/NbDtj');
     84    t=set(t,uid_NbDtj,'value',num2str(s.NbDtj))
     85    uid_NbDtk=find(t,'ImaDoc/Camera/BurstTiming/NbDtk');
     86    t=set(t,uid_NbDtk,'value',num2str(s.NbDtk))
     87    %%%
     88   
     89    save(t,[basename '.xml'])
    8190end
    8291
     
    8796    num_j=mod(ifile-1+first_label,nbfield2)+1;
    8897    num_i=floor((ifile-1+first_label)/nbfield2)+1;
    89     filename_new=name_generator(basename_new,num_i,num_j,'.png','_i_j');
    90     if test_level
    91         A=imread(filename);
    92         C=levels(A);
    93         imwrite(C,filename_new)
    94     else
    95         copyfile(filename,filename_new);
    96     end   
    97 end
    98 
     98    filename_new=name_generator(basename,num_i,num_j,'.png','_i_j');
     99    try
     100        movefile(filename,filename_new);
     101        [s,errormsg] = fileattrib(filename_new,'-w','a'); %set images to read only '-w' for all users ('a')
     102        if ~s
     103            msgbox_uvmat('ERROR',errormsg);
     104            return
     105        end
     106    catch
     107        errormsg=lasterr
     108        msgbox_uvmat('ERROR',errormsg);
     109        return
     110    end
     111    %     if test_level
     112    %         A=imread(filename);
     113    %         C=levels(A);
     114    %         imwrite(C,filename_new)
     115    %     else
     116    %         try
     117    %             copyfile(filename,filename_new);
     118    %         catch
     119    %             errormsg=lasterr
     120    %             msgbox_uvmat('ERROR',errormsg);
     121    %             return
     122    %         end
     123    %     end
     124end
    99125
    100126
     
    137163C=(C-cmin)/(cmax-cmin)*256;
    138164C=uint8(C);
     165
     166%'imadoc2struct': reads the xml file for image documentation
     167%------------------------------------------------------------------------
     168% function [s,errormsg]=imadoc2struct(ImaDoc,option)
     169%
     170% OUTPUT:
     171% s: structure representing ImaDoc
     172%   s.Heading: information about the data hierarchical structure
     173%   s.Time: matrix of times
     174%   s.TimeUnit
     175%  s.GeometryCalib: substructure containing the parameters for geometric calibration
     176% errormsg: error message
     177%
     178% INPUT:
     179% ImaDoc: full name of the xml input file with head key ImaDoc
     180% option: ='GeometryCalib': read  the data of GeometryCalib, including source point coordinates
     181
     182function [s,errormsg]=imadoc2struct(ImaDoc,option)
     183%% default input and output
     184if ~exist('option','var')
     185    option='*';
     186end
     187errormsg=[];%default
     188s.Heading=[];%default
     189s.Time=[]; %default
     190s.TimeUnit=[]; %default
     191s.GeometryCalib=[];
     192tsai=[];%default
     193
     194%% opening the xml file
     195if exist(ImaDoc,'file')~=2, errormsg=[ ImaDoc ' does not exist']; return;end;%input file does not exist
     196try
     197    t=xmltree(ImaDoc);
     198catch
     199    errormsg={[ImaDoc ' is not a valid xml file']; lasterr};
     200    display(errormsg);
     201    return
     202end
     203uid_root=find(t,'/ImaDoc');
     204if isempty(uid_root), errormsg=[ImaDoc ' is not an image documentation file ImaDoc']; return; end;%not an ImaDoc .xml file
     205
     206
     207%% Heading
     208uid_Heading=find(t,'/ImaDoc/Heading');
     209if ~isempty(uid_Heading),
     210    uid_Campaign=find(t,'/ImaDoc/Heading/Campaign');
     211    uid_Exp=find(t,'/ImaDoc/Heading/Experiment');
     212    uid_Device=find(t,'/ImaDoc/Heading/Device');
     213    uid_Record=find(t,'/ImaDoc/Heading/Record');
     214    uid_FirstImage=find(t,'/ImaDoc/Heading/ImageName');
     215    s.Heading.Campaign=get(t,children(t,uid_Campaign),'value');
     216    s.Heading.Experiment=get(t,children(t,uid_Exp),'value');
     217    s.Heading.Device=get(t,children(t,uid_Device),'value');
     218    if ~isempty(uid_Record)
     219        s.Heading.Record=get(t,children(t,uid_Record),'value');
     220    end
     221    s.Heading.ImageName=get(t,children(t,uid_FirstImage),'value');
     222end
     223
     224%% Camera  and timing
     225if strcmp(option,'*') || strcmp(option,'Camera')
     226    uid_Camera=find(t,'/ImaDoc/Camera');
     227    if ~isempty(uid_Camera)
     228        uid_ImageSize=find(t,'/ImaDoc/Camera/ImageSize');
     229        if ~isempty(uid_ImageSize);
     230            ImageSize=get(t,children(t,uid_ImageSize),'value');
     231            xindex=findstr(ImageSize,'x');
     232            if length(xindex)>=2
     233                s.Npx=str2double(ImageSize(1:xindex(1)-1));
     234                s.Npy=str2double(ImageSize(xindex(1)+1:xindex(2)-1));
     235            end
     236        end
     237        uid_TimeUnit=find(t,'/ImaDoc/Camera/TimeUnit');
     238        if ~isempty(uid_TimeUnit)
     239            s.TimeUnit=get(t,children(t,uid_TimeUnit),'value');
     240        end
     241        uid_BurstTiming=find(t,'/ImaDoc/Camera/BurstTiming');
     242        if ~isempty(uid_BurstTiming)
     243            for k=1:length(uid_BurstTiming)
     244                subt=branch(t,uid_BurstTiming(k));%subtree under BurstTiming
     245                % reading Dtk
     246                Frequency=get_value(subt,'/BurstTiming/FrameFrequency',1);
     247                Dtj=get_value(subt,'/BurstTiming/Dtj',[]);
     248                Dtj=Dtj/Frequency;%Dtj converted from frame unit to TimeUnit (e.g. 's')
     249                NbDtj=get_value(subt,'/BurstTiming/NbDtj',1);
     250                %%%% correction RDvision %%%%
     251                NbDtj=NbDtj/numel(Dtj);
     252                s.NbDtj=NbDtj;
     253                %%%%
     254                Dti=get_value(subt,'/BurstTiming/Dti',[]);
     255                Dti=Dti/Frequency;%Dtj converted from frame unit to TimeUnit (e.g. 's')
     256                NbDti=get_value(subt,'/BurstTiming/NbDti',1);
     257                Time_val=get_value(subt,'/BurstTiming/Time',0);%time in TimeUnit
     258                if ~isempty(Dti)
     259                    Dti=reshape(Dti'*ones(1,NbDti),NbDti*numel(Dti),1); %concatene Dti vector NbDti times
     260                    Time_val=[Time_val;Time_val(end)+cumsum(Dti)];%append the times defined by the intervals  Dti
     261                end
     262                if ~isempty(Dtj)
     263                    Dtj=reshape(Dtj'*ones(1,NbDtj),1,NbDtj*numel(Dtj)); %concatene Dtj vector NbDtj times
     264                    Dtj=[0 Dtj];
     265                    Time_val=Time_val*ones(1,numel(Dtj))+ones(numel(Time_val),1)*cumsum(Dtj);% produce a time matrix with Dtj
     266                end
     267                % reading Dtk
     268                Dtk=get_value(subt,'/BurstTiming/Dtk',[]);
     269                NbDtk=get_value(subt,'/BurstTiming/NbDtk',1);
     270                %%%% correction RDvision %%%%
     271                NbDtk=-1+(NbDtk+1)/(NbDti+1)
     272                s.NbDtk=NbDtk;
     273                %%%%%
     274                if isempty(Dtk)
     275                    s.Time=[s.Time;Time_val];
     276                else
     277                    for kblock=1:NbDtk+1
     278                        Time_val_k=Time_val+(kblock-1)*Dtk;
     279                        s.Time=[s.Time;Time_val_k];
     280                    end
     281                end
     282            end
     283        end
     284    end
     285end
     286
     287%% motor
     288if strcmp(option,'*') || strcmp(option,'GeometryCalib')
     289    uid_subtree=find(t,'/ImaDoc/TranslationMotor');
     290    if length(uid_subtree)==1
     291        subt=branch(t,uid_subtree);%subtree under GeometryCalib
     292       [s.TranslationMotor,errormsg]=read_subtree(subt,{'Nbslice','ZStart','ZEnd'},[1 1 1],[1 1 1]);
     293    end
     294end
     295%%  geometric calibration
     296if strcmp(option,'*') || strcmp(option,'GeometryCalib')
     297    uid_GeometryCalib=find(t,'/ImaDoc/GeometryCalib');
     298    if ~isempty(uid_GeometryCalib)
     299        if length(uid_GeometryCalib)>1
     300            errormsg=['More than one GeometryCalib in ' filecivxml];
     301            return
     302        end
     303        subt=branch(t,uid_GeometryCalib);%subtree under GeometryCalib
     304        cont=get(subt,1,'contents');
     305        if ~isempty(cont)
     306            uid_CalibrationType=find(subt,'/GeometryCalib/CalibrationType');
     307            if isequal(length(uid_CalibrationType),1)
     308                tsai.CalibrationType=get(subt,children(subt,uid_CalibrationType),'value');
     309            end
     310            uid_CoordUnit=find(subt,'/GeometryCalib/CoordUnit');
     311            if isequal(length(uid_CoordUnit),1)
     312                tsai.CoordUnit=get(subt,children(subt,uid_CoordUnit),'value');
     313            end
     314            uid_fx_fy=find(subt,'/GeometryCalib/fx_fy');
     315            focal=[];%default fro old convention (Reg Wilson)
     316            if isequal(length(uid_fx_fy),1)
     317                tsai.fx_fy=str2num(get(subt,children(subt,uid_fx_fy),'value'));
     318            else %old convention (Reg Wilson)
     319                uid_focal=find(subt,'/GeometryCalib/focal');
     320                uid_dpx_dpy=find(subt,'/GeometryCalib/dpx_dpy');
     321                uid_sx=find(subt,'/GeometryCalib/sx');
     322                if ~isempty(uid_focal) && ~isempty(uid_dpx_dpy) && ~isempty(uid_sx)
     323                    dpx_dpy=str2num(get(subt,children(subt,uid_dpx_dpy),'value'));
     324                    sx=str2num(get(subt,children(subt,uid_sx),'value'));
     325                    focal=str2num(get(subt,children(subt,uid_focal),'value'));
     326                    tsai.fx_fy(1)=sx*focal/dpx_dpy(1);
     327                    tsai.fx_fy(2)=focal/dpx_dpy(2);
     328                end
     329            end
     330            uid_Cx_Cy=find(subt,'/GeometryCalib/Cx_Cy');
     331            if ~isempty(uid_Cx_Cy)
     332                tsai.Cx_Cy=str2num(get(subt,children(subt,uid_Cx_Cy),'value'));
     333            end
     334            uid_kc=find(subt,'/GeometryCalib/kc');
     335            if ~isempty(uid_kc)
     336                tsai.kc=str2double(get(subt,children(subt,uid_kc),'value'));
     337            else %old convention (Reg Wilson)
     338                uid_kappa1=find(subt,'/GeometryCalib/kappa1');
     339                if ~isempty(uid_kappa1)&& ~isempty(focal)
     340                    kappa1=str2double(get(subt,children(subt,uid_kappa1),'value'));
     341                    tsai.kc=-kappa1*focal*focal;
     342                end
     343            end
     344            uid_Tx_Ty_Tz=find(subt,'/GeometryCalib/Tx_Ty_Tz');
     345            if ~isempty(uid_Tx_Ty_Tz)
     346                tsai.Tx_Ty_Tz=str2num(get(subt,children(subt,uid_Tx_Ty_Tz),'value'));
     347            end
     348            uid_R=find(subt,'/GeometryCalib/R');
     349            if ~isempty(uid_R)
     350                RR=get(subt,children(subt,uid_R),'value');
     351                if length(RR)==3
     352                    tsai.R=[str2num(RR{1});str2num(RR{2});str2num(RR{3})];
     353                end
     354            end
     355           
     356            %look for laser plane definitions
     357            uid_Angle=find(subt,'/GeometryCalib/PlaneAngle');
     358            uid_Pos=find(subt,'/GeometryCalib/SliceCoord');
     359            if isempty(uid_Pos)
     360                uid_Pos=find(subt,'/GeometryCalib/PlanePos');%old convention
     361            end
     362            if ~isempty(uid_Angle)
     363                tsai.PlaneAngle=str2num(get(subt,children(subt,uid_Angle),'value'));
     364            end
     365            if ~isempty(uid_Pos)
     366                for j=1:length(uid_Pos)
     367                    tsai.SliceCoord(j,:)=str2num(get(subt,children(subt,uid_Pos(j)),'value'));
     368                end
     369                uid_DZ=find(subt,'/GeometryCalib/SliceDZ');
     370                uid_NbSlice=find(subt,'/GeometryCalib/NbSlice');
     371                if ~isempty(uid_DZ) && ~isempty(uid_NbSlice)
     372                    DZ=str2double(get(subt,children(subt,uid_DZ),'value'));
     373                    NbSlice=get(subt,children(subt,uid_NbSlice),'value');
     374                    if isequal(NbSlice,'volume')
     375                        tsai.NbSlice='volume';
     376                        NbSlice=NbDtj+1;
     377                    else
     378                        tsai.NbSlice=str2double(NbSlice);
     379                    end
     380                    tsai.SliceCoord=ones(NbSlice,1)*tsai.SliceCoord+DZ*(0:NbSlice-1)'*[0 0 1];
     381                end
     382            end   
     383            tsai.SliceAngle=get_value(subt,'/GeometryCalib/SliceAngle',[0 0 0]);
     384            tsai.VolumeScan=get_value(subt,'/GeometryCalib/VolumeScan','n');
     385            tsai.InterfaceCoord=get_value(subt,'/GeometryCalib/InterfaceCoord',[0 0 0]);
     386            tsai.RefractionIndex=get_value(subt,'/GeometryCalib/RefractionIndex',1);
     387           
     388            if strcmp(option,'GeometryCalib')
     389                tsai.PointCoord=get_value(subt,'/GeometryCalib/SourceCalib/PointCoord',[0 0 0 0 0]);
     390            end
     391            s.GeometryCalib=tsai;
     392        end
     393    end
     394end
     395
     396%--------------------------------------------------
     397%  read a subtree
     398% INPUT:
     399% t: xltree
     400% head_element: head elelemnt of the subtree
     401% Data, structure containing
     402%    .Key: element name
     403%    .Type: type of element ('charg', 'float'....)
     404%    .NbOccur: nbre of occurrence, NaN for un specified number
     405function [s,errormsg]=read_subtree(subt,Data,NbOccur,NumTest)
     406%--------------------------------------------------
     407s=[];%default
     408errormsg='';
     409head_element=get(subt,1,'name');
     410    cont=get(subt,1,'contents');
     411    if ~isempty(cont)
     412        for ilist=1:length(Data)
     413            uid_key=find(subt,[head_element '/' Data{ilist}]);
     414            if ~isequal(length(uid_key),NbOccur(ilist))
     415                errormsg=['wrong number of occurence for ' Data{ilist}];
     416                return
     417            end
     418            for ival=1:length(uid_key)
     419                val=get(subt,children(subt,uid_key(ival)),'value');
     420                if ~NumTest(ilist)
     421                    eval(['s.' Data{ilist} '=val;']);
     422                else
     423                    eval(['s.' Data{ilist} '=str2double(val);'])
     424                end
     425            end
     426        end
     427    end
     428
     429
     430%--------------------------------------------------
     431%  read an xml element
     432function val=get_value(t,label,default)
     433%--------------------------------------------------
     434val=default;
     435uid=find(t,label);%find the element iud(s)
     436if ~isempty(uid) %if the element named label exists
     437   uid_child=children(t,uid);%find the children
     438   if ~isempty(uid_child)
     439       data=get(t,uid_child,'type');%get the type of child
     440       if iscell(data)% case of multiple element
     441           for icell=1:numel(data)
     442               val_read=str2num(get(t,uid_child(icell),'value'));
     443               if ~isempty(val_read)
     444                   val(icell,:)=val_read;
     445               end
     446           end
     447%           val=val';
     448       else % case of unique element value
     449           val_read=str2num(get(t,uid_child,'value'));
     450           if ~isempty(val_read)
     451               val=val_read;
     452           else
     453              val=get(t,uid_child,'value');%char string data
     454           end
     455       end
     456   end
     457end
  • trunk/src/series/sub_background.m

    r169 r214  
    5050        %'GetObject';...;%can use projection object
    5151        %'GetMask';...;%can use mask option 
    52         'PARAMETER';'NbSliding';...
    53         'PARAMETER';'VolumeScan';...
    54         'PARAMETER';'RankBrightness';...
     52        %'PARAMETER';'NbSliding';...
     53        %'PARAMETER';'VolumeScan';...
     54        %'PARAMETER';'RankBrightness';...
    5555               ''};
    5656    return %exit the function
     
    5858
    5959%----------------------------------------------------------------
    60 % initiate the waitbar
     60%% initiate the waitbar
    6161hseries=guidata(Series.hseries);%handles of the GUI series
    6262WaitbarPos=get(hseries.waitbar_frame,'Position');
     
    6767end
    6868
    69 %determine input image type
     69%% determine input image type
    7070FileType=[];%default
    7171MovieObject=[];
     
    101101nbaver_init=23;%approximate number of images used for the sliding background: to be adjusted later to include an integer number of bursts
    102102
    103 %adjust the proposed number of images in the sliding average to include an integer number of bursts
     103
     104%% apply the image rescaling function 'level' (avoid the blinking effects of bright particles)
     105answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m after sub_background');
     106test_level=isequal(answer,'Yes');
     107
     108%% adjust the proposed number of images in the sliding average to include an integer number of bursts
    104109if siz(2)~=1
    105110    nbaver=floor(nbaver_init/siz(1)); % number of bursts used for the sliding background,
     
    114119nom_type=Series.NomType;
    115120
    116 %create dir of the new images
    117 [dir_images,namebase]=fileparts(filebase);
    118 [path,subdir_ima]=fileparts(dir_images);
    119 curdir=pwd;
    120 cd(path);
    121 mkdir([subdir_ima '_b']);
    122 [xx,msg2] = fileattrib([subdir_ima '_b'],'+w','g'); %yield writing access (+w) to user group (g)
     121%% create dir of the new images
     122% [dir_images,namebase]=fileparts(filebase);
     123if test_level
     124    term='_b_levels';
     125else
     126    term='_b';
     127end
     128[pp,subdir_ima]=fileparts(Series.RootPath);
     129try
     130    mkdir([dir_images term]);
     131catch
     132    errormsg=lasterr
     133            msgbox_uvmat('ERROR',errormsg);
     134            return
     135end
     136[xx,msg2] = fileattrib([dir_images term],'+w','g'); %yield writing access (+w) to user group (g)
    123137if ~strcmp(msg2,'')
    124     msgbox_uvmat('ERROR',['pb of permission for ' subdir_ima '_b' ': ' msg2])%error message for directory creation
    125     cd(curdir)
     138    msgbox_uvmat('ERROR',['pb of permission for ' subdir_ima term ': ' msg2])%error message for directory creation
    126139    return
    127140end
    128 cd(curdir);
    129 filebase_b=fullfile(path,[subdir_ima '_b'],namebase);
    130 
     141filebase_b=fullfile([dir_images term],Series.RootFile);
     142
     143%% set processing parameters
    131144prompt = {'Number of images for the sliding background';'The number of positions (laser slices)';'volume scan mode (Yes/No)';...
    132145        'the luminosity rank chosen to define the background (0.1=for dense particle seeding, 0.5 (median) for sparse particles'};
    133 dlg_title = ['get (slice by slice) a sliding background and substract to each image, result in subdir ' subdir_ima '_b'];
     146dlg_title = ['get (slice by slice) a sliding background and substract to each image, result in subdir ' subdir_ima term];
    134147num_lines= 3;
    135148def     = { num2str(nbaver_init);num2str(nbslice_i);'No';'0.1'};
     
    139152
    140153nbaver_ima=str2num(answer{1});%number of images for the sliding background
    141 nbaver=ceil(nbaver_ima/siz(1))%number of bursts for the sliding background
     154nbaver=ceil(nbaver_ima/siz(1));%number of bursts for the sliding background
    142155if isequal(floor(nbaver/2),nbaver)
    143156   nbaver=nbaver+1%put the number of burst to an odd number (so the middle burst is defined)
    144157end
    145 % if isequal(nbaver,round(nbaver))
    146158step=siz(1);%case of bursts: the sliding background is shifted by one burst
    147 % else
    148 %    step=1;
    149 % end
    150159vol_test=answer{3};
    151160if isequal(vol_test,'Yes')
     
    164173    end
    165174end
    166 % nbaver=floor(nbaver_init/nbfield2); % number of bursts used for the sliding background,
    167 % if isequal(floor(nbaver/2),nbaver)
    168 %     nbaver=nbaver+1;%put the number of burst to an odd number
    169 % end
    170 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    171 %nbaver_ima=nbaver*nbfield2;% adjust the number of sliding images A  REMETRE
    172 % if ~isequal(nbaver_ima,nbaver_init)
    173 %     hwarn=warndlg(['number of images in the sliding average adjusted to ' num2str(nbaver_ima)]);
    174 %     set(hwarn,'Units','normalized')
    175 %     set(hwarn,'Position',[0.3 0.3 0.4 0.1])
    176 % end
    177 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    178175rank=floor(str2num(answer{4})*nbaver_ima);
    179176if rank==0
     
    183180nbfield=floor(lengthtot/(nbfield2*nbslice_i));%total number of i indexes (adjusted to an integer number of slices)
    184181nbfield_slice=nbfield*nbfield2;% number of fields per slice
    185 % test_plot=isequal(answer{5},'Yes'); %=1 to display background images
    186182if nbaver_ima > nbfield*nbfield2
    187183    errordlg('number of images in a slice smaller than the proposed number of images for the sliding average')
     
    190186nbfirst=(ceil(nbaver/2))*step;
    191187if nbfirst>nbaver_ima
    192     nbfirst=ceil(nbaver_ima/2)
     188    nbfirst=ceil(nbaver_ima/2);
    193189    step=1;
    194190    nbaver=nbaver_ima;
    195191end
    196192
    197 %copy the xml file
     193%% copy the xml file
    198194if exist([filebase '.xml'],'file')
    199195    copyfile([filebase '.xml'],[filebase_b '.xml']);% copy the .civ file
     
    223219    [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
    224220    [t]=add(t,NameFunction_uid,'chardata','sub_background');     
     221    if test_levels
     222            [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
     223            [t]=add(t,NameFunction_uid,'chardata','levels');
     224    end
    225225    [t,NbSlice_uid]=add(t,new_uid,'element','NbSlice');
    226226    [t]=add(t,new_uid,'chardata',num2str(nbslice_i));
     
    240240%MAIN LOOP ON SLICES
    241241for islice=1:nbslice_i
    242     %select the series of image indices at the level islice
     242    %% select the series of image indices at the level islice
    243243    for ifield=1:nbfield
    244244        for iburst=1:nbfield2
     
    246246        end
    247247    end 
    248     %read the first series of nbaver_ima images and sort by luminosity at each pixel
     248   
     249    %% read the first series of nbaver_ima images and sort by luminosity at each pixel
    249250    for ifield = 1:nbaver_ima
    250251        ifile=indselect(ifield);
     
    255256    Asort=sort(Ak,3);%sort the luminosity of images at each point
    256257    B=Asort(:,:,rank);%background image
    257 %
    258 %     namemean=name_generator([filebase '_back'],islice,[],'.png','_i');% makes the file name
    259 %     imwrite(B,namemean,'BitDepth',16); % save the first background image
    260 %     ['background image for slice ' num2str(islice) ' saved in ' namemean]
    261     %substract the background from each of the first images and save the new images
    262 %     for ifield=1:floor(nbaver_ima/2)+1
    263     'first background image will be substracted'
    264 
     258   display( 'first background image will be substracted')
    265259    for ifield=1:nbfirst
    266260            Acor=double(Ak(:,:,ifield))-double(B);%substract background to the current image
     
    269263            ifile=indselect(ifield);
    270264            newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name
    271             imwrite(C,newname,'BitDepth',16); % save the new image
    272     end
    273     %repeat the operation on a sliding series of nbaver*nbfield2 images
    274     'sliding background image will be substracted'
     265            if test_level
     266                C=levels(C);
     267                 imwrite(C,newname,'BitDepth',8); % save the new image
     268            else
     269                 imwrite(C,newname,'BitDepth',16); % save the new image
     270            end
     271    end
     272   
     273    %% repeat the operation on a sliding series of nbaver*nbfield2 images
     274    display('sliding background image will be substracted')
    275275    if nbfield_slice > nbaver_ima
    276 %         for ifield = floor(nbaver_ima/2)+2:step:nbfield_slice-floor(nbaver_ima/2)
    277276        for ifield = step*ceil(nbaver/2)+1:step:nbfield_slice-step*floor(nbaver/2)
    278277            stopstate=get(hseries.RUN,'BusyAction');
     
    291290                B=Asort(:,:,rank);%current background image
    292291                for iburst=1:step
    293 %                     Acor=double(Ak(:,:,floor(nbaver_ima/2)+iburst-1))-double(B);
    294292                    index=step*floor(nbaver/2)+iburst;
    295293                    Acor=double(Ak(:,:,index))-double(B);
     
    298296                    ifile=indselect(ifield+iburst-1);
    299297                    [newname]=...
    300                        name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',Series.NomType)
    301 %                     newname=name_generator(filebase_b,num_i1(indselect(ifield+iburst-1)),num_j1(indselect(ifield+iburst-1)),'.png',nom_type)% makes the new file name
    302                     imwrite(C,newname,'BitDepth',16); % save the new image
    303                 end 
     298                        name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',Series.NomType) % makes the new file name
     299                    if test_level
     300                        C=levels(C);
     301                        imwrite(C,newname,'BitDepth',8); % save the new image
     302                    else
     303                        imwrite(C,newname,'BitDepth',16); % save the new image
     304                    end
     305                end
    304306            else
    305307                return
     
    308310    end
    309311
    310 %substract the background from the last images
    311 %     for ifield=nbfield_slice-floor(nbaver_ima/2)+1:nbfield_slice
    312      'last background image will be substracted'
     312%% substract the background from the last images
     313    display('last background image will be substracted')
    313314     ifield=nbfield_slice-(step*ceil(nbaver/2))+1:nbfield_slice;
    314     for ifield=nbfield_slice-(step*floor(nbaver/2))+1:nbfield_slice 
    315         index=ifield-nbfield_slice+step*(2*floor(nbaver/2)+1);
    316         Acor=double(Ak(:,:,index))-double(B);
    317         Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
    318         C=uint16(Acor);
    319         ifile=indselect(ifield);
    320         newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name
    321         imwrite(C,newname,'BitDepth',16); % save the new image
    322     end
     315     for ifield=nbfield_slice-(step*floor(nbaver/2))+1:nbfield_slice
     316         index=ifield-nbfield_slice+step*(2*floor(nbaver/2)+1);
     317         Acor=double(Ak(:,:,index))-double(B);
     318         Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
     319         C=uint16(Acor);
     320         ifile=indselect(ifield);
     321         newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name
     322         if test_level
     323             C=levels(C);
     324             imwrite(C,newname,'BitDepth',8); % save the new image
     325         else
     326             imwrite(C,newname,'BitDepth',16); % save the new image
     327         end
     328     end
    323329end
    324330
     
    348354end
    349355   
     356
     357function C=levels(A)
     358%whos A;
     359B=double(A(:,:,1));
     360windowsize=round(min(size(B,1),size(B,2))/20);
     361windowsize=floor(windowsize/2)*2+1;
     362ix=[1/2-windowsize/2:-1/2+windowsize/2];%
     363%del=np/3;
     364%fct=exp(-(ix/del).^2);
     365fct2=cos(ix/(windowsize-1)/2*pi/2);
     366%Mfiltre=(ones(5,5)/5^2);
     367%Mfiltre=fct2';
     368Mfiltre=fct2'*fct2;
     369Mfiltre=Mfiltre/(sum(sum(Mfiltre)));
     370
     371C=filter2(Mfiltre,B);
     372C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize);
     373C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize);
     374C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:);
     375C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:);
     376C=tanh(B./(2*C));
     377[n,c]=hist(reshape(C,1,[]),100);
     378% figure;plot(c,n);
     379
     380[m,i]=max(n);
     381c_max=c(i);
     382[dummy,index]=sort(abs(c-c(i)));
     383n=n(index);
     384c=c(index);
     385i_select = find(cumsum(n)<0.95*sum(n));
     386if isempty(i_select)
     387    i_select = 1:length(c);
     388end
     389c_select=c(i_select);
     390n_select=n(i_select);
     391cmin=min(c_select);
     392cmax=max(c_select);
     393C=(C-cmin)/(cmax-cmin)*256;
     394C=uint8(C);
Note: See TracChangeset for help on using the changeset viewer.