Changeset 950 for trunk/src


Ignore:
Timestamp:
Jun 11, 2016, 9:32:29 PM (9 years ago)
Author:
sommeria
Message:

various updates

Location:
trunk/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/imadoc2struct.m

    r924 r950  
    126126        end
    127127    end
     128    if ~isfield(s.Camera,'FirstFrameIndexI')
     129        s.Camera.FirstFrameIndexI=1; %first index qssumed equl to 1 by default
     130    end
    128131    s.Time=[zeros(size(s.Time,1),1) s.Time]; %insert a vertical line of zeros (to deal with zero file indices)
    129     s.Time=[zeros(1,size(s.Time,2)); s.Time]; %insert a horizontal line of zeros
     132    if s.Camera.FirstFrameIndexI~=0
     133    s.Time=[zeros(s.Camera.FirstFrameIndexI,size(s.Time,2)); s.Time]; %insert a horizontal line of zeros
     134    end
    130135end
    131136

  • trunk/src/script_readlvm.m

    r893 r950  
    11%% get the input file
    2 fileinput=uigetfile_uvmat('pick an input file','/fsnet/project/coriolis/2015/15STRATJET/Probes');
     2fileinput=uigetfile_uvmat('pick an input file','/fsnet/project/coriolis/2016/16MILESTONE/Data');
    33[Path,Name,Ext]=fileparts(fileinput);
    44

  • trunk/src/series.m

    r942 r950  
    990990        Time=XmlData.Time;
    991991        TimeName='xml';
     992                        if XmlData.Time(1,:)==XmlData.Time(2,:)% case starting with index 1
     993                    sizDti=size(XmlData.Time,1)-1;%size of the time vector explicitly defined in the xml file
     994                    ind_start=1;
     995                else
     996                    sizDti=size(XmlData.Time,1);% case starting with index 0
     997                    ind_start=0;
     998                end
     999        % complement the input if the whole time series is not defined
     1000            if size(i1_series,3)>size(XmlData.Time,1)-ind_start %only the first time interval is defined, extrapolate to the whole series
     1001                Dti_total=XmlData.Time(end)-XmlData.Time(1);%total time interval covered by the time vector
     1002                missing_indices=sizDti+1+ind_start:size(i1_series,3)+1;% remaining set of frame indices for which time needs to be found
     1003                repeat_nbre=1+floor((missing_indices-sizDti-ind_start)/(sizDti-1));% number of repetitions of Dti
     1004                time_indices=1+mod(missing_indices-sizDti-1,sizDti-1);
     1005                for j=1:size(XmlData.Time,2)
     1006                Time(missing_indices,j)=XmlData.Time(time_indices,j)+repeat_nbre'*Dti_total;
     1007                end
     1008                % update the xml file with NbDti
     1009                t=xmltree(XmlFileName);
     1010                uid_NbDti=find(t,'ImaDoc/Camera/BurstTiming/NbDti')
     1011                if isempty(uid_NbDti)
     1012                    uid_BurstTiming=find(t,'ImaDoc/Camera/BurstTiming')
     1013                    [t,uid_NbDti]=add(t,uid_BurstTiming,'element','NbDti');
     1014                end
     1015                [t,uid_NbDti]=add(t,uid_NbDti,'chardata',num2str(repeat_nbre(end)-1));
     1016                save(t,XmlFileName)
     1017            end
    9921018    end
    9931019    if isfield(XmlData,'Camera')

  • trunk/src/series/civ_series.m

    r949 r950  
    283283if iview_A~=0
    284284    XmlFileName=find_imadoc(RootPath_A,SubDir_A,RootFile_A,FileExt_A);
    285     time=[];
     285    Time=[];
    286286    if ~isempty(XmlFileName)
    287287        XmlData=imadoc2struct(XmlFileName);
    288288        if isfield(XmlData,'Time')
    289             time=XmlData.Time;
     289            Time=XmlData.Time;
    290290            TimeSource='xml';
    291291        end
    … …  
    302302        end
    303303    end
    304     if isempty(time) && ~isempty(find(strcmp(FileType_A,{'mmreader','video'})))% case of video input
    305         time=zeros(FileInfo_A.NumberOfFrames+1,2);
    306         time(:,2)=(0:1/FileInfo_A.FrameRate:(FileInfo_A.NumberOfFrames)/FileInfo_A.FrameRate)';
     304    if isempty(Time) && ~isempty(find(strcmp(FileType_A,{'mmreader','video'})))% case of video input
     305        Time=zeros(FileInfo_A.NumberOfFrames+1,2);
     306        Time(:,2)=(0:1/FileInfo_A.FrameRate:(FileInfo_A.NumberOfFrames)/FileInfo_A.FrameRate)';
    307307        TimeSource='video';
    308308        ColorType='truecolor';
    309309    end
    310     if isempty(time)% time = index i +0.001 index j by default
     310    if isempty(Time)% Time = index i +0.001 index j by default
    311311        %MinIndex_i=min(i1_series_Civ1);
    312312        MaxIndex_i=max(i2_series_Civ1);
    313313        %MinIndex_j=min(j1_series_Civ1);
    314314        MaxIndex_j=max(j2_series_Civ1);
    315         time=(1:MaxIndex_i)'*ones(1,MaxIndex_j);
    316         time=time+0.001*ones(MaxIndex_i,1)*(1:MaxIndex_j);
    317         time=[zeros(1,MaxIndex_j);time];% insert a first line of zeros
    318         time=[zeros(MaxIndex_i+1,1) time];% insert a first column of zeros
     315        Time=(1:MaxIndex_i)'*ones(1,MaxIndex_j);
     316        Time=Time+0.001*ones(MaxIndex_i,1)*(1:MaxIndex_j);
     317        Time=[zeros(1,MaxIndex_j);Time];% insert a first line of zeros
     318        Time=[zeros(MaxIndex_i+1,1) Time];% insert a first column of zeros
    319319    end
    320320   
    … …  
    438438            end
    439439            if strcmp(Param.ActionInput.ListCompareMode,'displacement')
    440                 Data.Civ1_Time=time(i2+1,j2+1);% the time is the time of the secodn image
    441                 Data.Civ1_Dt=1;% time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading
     440                Data.Civ1_Time=Time(i2+1,j2+1);% the Time is the Time of the secodn image
     441                Data.Civ1_Dt=1;% Time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading
    442442            else
    443             Data.Civ1_Time=(time(i2+1,j2+1)+time(i1+1,j1+1))/2;% the time is the time at the middle of the image pair
    444             Data.Civ1_Dt=time(i2+1,j2+1)-time(i1+1,j1+1);
     443            Data.Civ1_Time=(Time(i2+1,j2+1)+Time(i1+1,j1+1))/2;% the Time is the Time at the middle of the image pair
     444            Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
    445445            end
    446446            for ilist=1:length(list_param)
    … …  
    736736                Civ2_Dt=1;
    737737            else
    738                 Civ2_Dt=time(i2_civ2+1,j2_civ2+1)-time(i1_civ2+1,j1_civ2+1);
     738                Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1);
    739739            end
    740740        end
    … …  
    761761            Data.Civ2_ImageB=ImageName_B;
    762762             if strcmp(Param.ActionInput.ListCompareMode,'displacement')
    763                 Data.Civ2_Time=time(i2_civ2+1,j2_civ2+1);% the time is the time of the secodn image
    764                 Data.Civ2_Dt=1;% time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading
     763                Data.Civ2_Time=Time(i2_civ2+1,j2_civ2+1);% the Time is the Time of the secodn image
     764                Data.Civ2_Dt=1;% Time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading
    765765             else
    766             Data.Civ2_Time=(time(i2_civ2+1,j2_civ2+1)+time(i1_civ2+1,j1_civ2+1))/2;
     766            Data.Civ2_Time=(Time(i2_civ2+1,j2_civ2+1)+Time(i1_civ2+1,j1_civ2+1))/2;
    767767            Data.Civ2_Dt=Civ2_Dt;
    768768             end

  • trunk/src/series/sub_background.m

    r924 r950  
    340340        if Param.ActionInput.CheckLevelTransform
    341341            C=levels(Acor);
    342             imwrite(C,newname,'BitDepth',8); % save the new image
     342            imwrite(C,newname,'BitDepth',16); % save the new image
    343343        else
    344344            if isequal(FileInfo{1}.BitDepth,16)
    … …  
    393393                if Param.ActionInput.CheckLevelTransform
    394394                    C=levels(Acor);
    395                     imwrite(C,newname,'BitDepth',8); % save the new image
     395                    imwrite(C,newname,'BitDepth',16); % save the new image
    396396                else
    397397                    if isequal(FileInfo{1}.BitDepth,16)
    … …  
    421421        if Param.ActionInput.CheckLevelTransform
    422422            C=levels(Acor);
    423             imwrite(C,newname,'BitDepth',8); % save the new image
     423            imwrite(C,newname,'BitDepth',16); % save the new image
    424424        else
    425425            if isequal(FileInfo{1}.BitDepth,16)
    … …  
    435435end
    436436
    437 
    438437function C=levels(A)
    439 %whos A;
    440 B=double(A(:,:,1));
    441 windowsize=round(min(size(B,1),size(B,2))/20);
    442 windowsize=floor(windowsize/2)*2+1;
    443 ix=1/2-windowsize/2:-1/2+windowsize/2;%
    444 %del=np/3;
    445 %fct=exp(-(ix/del).^2);
    446 fct2=cos(ix/(windowsize-1)/2*pi/2);
    447 %Mfiltre=(ones(5,5)/5^2);
    448 %Mfiltre=fct2';
    449 Mfiltre=fct2'*fct2;
    450 Mfiltre=Mfiltre/(sum(sum(Mfiltre)));
    451 
    452 C=filter2(Mfiltre,B);
    453 C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize);
    454 C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize);
    455 C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:);
    456 C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:);
    457 C=tanh(B./(2*C));
    458 [n,c]=hist(reshape(C,1,[]),100);
    459 % figure;plot(c,n);
    460 
    461 [m,i]=max(n);
    462 c_max=c(i);
    463 [dummy,index]=sort(abs(c-c(i)));
    464 n=n(index);
    465 c=c(index);
    466 i_select = find(cumsum(n)<0.95*sum(n));
    467 if isempty(i_select)
    468     i_select = 1:length(c);
    469 end
    470 c_select=c(i_select);
    471 n_select=n(i_select);
    472 cmin=min(c_select);
    473 cmax=max(c_select);
    474 C=(C-cmin)/(cmax-cmin)*256;
    475 C=uint8(C);
     438
     439nblock_y=100;%2*Param.TransformInput.BlockSize;
     440nblock_x=100;%2*Param.TransformInput.BlockSize;
     441[npy,npx]=size(A);
     442[X,Y]=meshgrid(1:npx,1:npy);
     443
     444%Backg=zeros(size(A));
     445%Aflagmin=sparse(imregionalmin(A));%Amin=1 for local image minima
     446%Amin=A.*Aflagmin;%values of A at local minima
     447% local background: find all the local minima in image subblocks
     448fctblock= inline('median(x(:))');
     449Backg=blkproc(A,[nblock_y nblock_x],fctblock);% take the median in  blocks
     450fctblock= inline('mean(x(:))');
     451B=imresize(Backg,size(A),'bilinear');% interpolate to the initial size image
     452A=(A-B);%substract background
     453AMean=blkproc(A,[nblock_y nblock_x],fctblock);% take the mean in  blocks
     454fctblock= inline('var(x(:))');
     455AVar=blkproc(A,[nblock_y nblock_x],fctblock);% take the mean in  blocks
     456Avalue=AVar./AMean% typical value of particle luminosity
     457Avalue=imresize(Avalue,size(A),'bilinear');% interpolate to the initial size image
     458C=uint16(1000*tanh(A./(2*Avalue)));
     459%Bmin=blkproc(Aflagmin,[nblock_y nblock_x],sumblock);% find the number of minima in blocks
     460%Backg=Backg./Bmin; % find the average of minima in blocks
     461% function C=levels(A)
     462% %whos A;
     463% B=double(A(:,:,1));
     464% windowsize=round(min(size(B,1),size(B,2))/20);
     465% windowsize=floor(windowsize/2)*2+1;
     466% ix=1/2-windowsize/2:-1/2+windowsize/2;%
     467% %del=np/3;
     468% %fct=exp(-(ix/del).^2);
     469% fct2=cos(ix/(windowsize-1)/2*pi/2);
     470% %Mfiltre=(ones(5,5)/5^2);
     471% %Mfiltre=fct2';
     472% Mfiltre=fct2'*fct2;
     473% Mfiltre=Mfiltre/(sum(sum(Mfiltre)));
     474%
     475% C=filter2(Mfiltre,B);
     476% C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize);
     477% C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize);
     478% C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:);
     479% C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:);
     480% C=tanh(B./(2*C));
     481% [n,c]=hist(reshape(C,1,[]),100);
     482% % figure;plot(c,n);
     483%
     484% [m,i]=max(n);
     485% c_max=c(i);
     486% [dummy,index]=sort(abs(c-c(i)));
     487% n=n(index);
     488% c=c(index);
     489% i_select = find(cumsum(n)<0.95*sum(n));
     490% if isempty(i_select)
     491%     i_select = 1:length(c);
     492% end
     493% c_select=c(i_select);
     494% n_select=n(i_select);
     495% cmin=min(c_select);
     496% cmax=max(c_select);
     497% C=(C-cmin)/(cmax-cmin)*256;
     498% C=uint8(C);

  • trunk/src/transform_field/ima_levels.m

    r924 r950  
     1% 'ima_remove_background': removes backgound from an image (using the local minimum)
     2% requires the Matlab image processing toolbox
     3%------------------------------------------------------------------------
     4%%%%  Use the general syntax for transform fields with a single input %%%%
     5% OUTPUT:
     6% DataOut:   output field structure
     7%
     8%INPUT:
     9% DataIn:  first input field structure
     10
    111%=======================================================================
    212% Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
    … …  
    1727%=======================================================================
    1828
    19 function DataOut=im_levels(DataIn)
    20 %% set GUI config: no action defined
    21 DataOut=[];  %default  output field
    22 if strcmp(DataIn,'*')
     29function DataOut=ima_remove_background_blocks(DataIn,Param)
     30%------------------------------------------------------------------------
     31%% request input parameters
     32if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
     33    prompt = {'block size(pixels)'};
     34    dlg_title = 'get the block size (in pixels) used to calculate the local statistics';
     35    num_lines= 1;
     36    def     = { '100'};
     37    if isfield(Param,'TransformInput')&&isfield(Param.TransformInput,'BlockSize')
     38        def={num2str(Param.TransformInput.BlockSize)};
     39    end
     40    answer = inputdlg(prompt,dlg_title,num_lines,def);
     41    DataOut.TransformInput.BlockSize=str2num(answer{1});
    2342    return
    2443end
    25 %-----------------------------------------------
    26 %parameters
    27 np=30
     44if ~isfield(DataIn,'A')
     45    DataOut.Txt='remove_particles only valid for input images';
     46    return
     47end
     48if ~exist('imerode','file');
     49        DataOut.Txt='the function imerode from the image processing toolbox is needed';
     50    return
     51end
     52
    2853%---------------------------------------------------------
    2954DataOut=DataIn;%default
     55nblock_y=2*Param.TransformInput.BlockSize;
     56nblock_x=2*Param.TransformInput.BlockSize;
     57[npy,npx]=size(DataIn.A);
     58[X,Y]=meshgrid(1:npx,1:npy);
    3059
    31 B=double(DataIn.A(:,:,1));
    32 windowsize=round(min(size(B,1),size(B,2))/20);
    33 windowsize=floor(windowsize/2)*2+1;
    34 ix=[1/2-windowsize/2:-1/2+windowsize/2];%
    35 %del=np/3;
    36 %fct=exp(-(ix/del).^2);
    37 fct2=cos(ix/((np-1)/2)*pi/2);
    38 %Mfiltre=(ones(5,5)/5^2);
    39 %Mfiltre=fct2';
    40 Mfiltre=fct2'*fct2;
    41 Mfiltre=Mfiltre/(sum(sum(Mfiltre)));
    42 
    43 C=filter2(Mfiltre,B);
    44 C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize);
    45 C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize);
    46 C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:);
    47 C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:);
    48 C=tanh(B./(2*C));
    49 [n,c]=hist(reshape(C,1,[]),100);
    50 % figure;plot(c,n);
    51 
    52 [m,i]=max(n);
    53 c_max=c(i);
    54 [dummy,index]=sort(abs(c-c(i)));
    55 n=n(index);
    56 c=c(index);
    57 i_select = find(cumsum(n)<0.95*sum(n));
    58 if isempty(i_select)
    59     i_select = 1:length(c);
    60 end
    61 c_select=c(i_select);
    62 n_select=n(i_select);
    63 cmin=min(c_select);
    64 cmax=max(c_select);
    65 C=(C-cmin)/(cmax-cmin)*256;
    66 DataOut.AA=uint8(C);
     60%BACKGROUND LEVEL
     61Atype=class(DataIn.A);
     62A=double(DataIn.A);
     63%Backg=zeros(size(A));
     64%Aflagmin=sparse(imregionalmin(A));%Amin=1 for local image minima
     65%Amin=A.*Aflagmin;%values of A at local minima
     66% local background: find all the local minima in image subblocks
     67fctblock= inline('median(x(:))');
     68Backg=blkproc(A,[nblock_y nblock_x],fctblock);% take the median in  blocks
     69fctblock= inline('mean(x(:))');
     70B=imresize(Backg,size(A),'bilinear');% interpolate to the initial size image
     71A=(A-B);%substract background
     72AMean=blkproc(A,[nblock_y nblock_x],fctblock);% take the mean in  blocks
     73fctblock= inline('var(x(:))');
     74AVar=blkproc(A,[nblock_y nblock_x],fctblock);% take the mean in  blocks
     75Avalue=AVar./AMean% typical value of particle luminosity
     76Avalue=imresize(Avalue,size(A),'bilinear');% interpolate to the initial size image
     77DataOut.A=uint16(1000*tanh(A./(2*Avalue)));
     78%Bmin=blkproc(Aflagmin,[nblock_y nblock_x],sumblock);% find the number of minima in blocks
     79%Backg=Backg./Bmin; % find the average of minima in blocks
    6780
    6881
    69 
    70 DataOut.A=uint8(C);

  • trunk/src/uvmat.m

    r946 r950  
    19701970UvData.j1_series{index}=j1_series;
    19711971UvData.j2_series{index}=j2_series;
     1972nbfield=max(max(max(i2_series)));% total number of fields (i index)
     1973if isempty(nbfield)
     1974    nbfield=max(max(max(i1_series)));
     1975end
     1976nbfield_j=max(max(max(j2_series)));% number of fields along j index
     1977if isempty(nbfield_j)
     1978    nbfield_j=max(max(max(j1_series)));
     1979end
    19721980
    19731981%% read timing and total frame number from the current file (e.g. movie files)
    … …  
    20202028        if isfield(XmlDataRead,'TimeUnit')&& ~isempty(XmlDataRead.TimeUnit)
    20212029            TimeUnit=XmlDataRead.TimeUnit;
    2022         end
     2030        end   
    20232031        if isfield(XmlDataRead,'Time')&& ~isempty(XmlDataRead.Time)
    20242032            XmlData.Time=XmlDataRead.Time;
     2033             if XmlDataRead.Time(1,:)==XmlDataRead.Time(2,:)% case starting with index 1
     2034            sizDti=size(XmlDataRead.Time,1)-1;%size of the time vector explicitly defined in the xml file
     2035            ind_start=1;
     2036        else
     2037            sizDti=size(XmlDataRead.Time,1);% case starting with index 0
     2038            ind_start=0;
     2039        end
     2040            % complement the input if the whole time series is not defined
     2041            if size(i1_series,3)>size(XmlDataRead.Time,1)-ind_start %only the first time interval is defined, extrapolate to the whole series           
     2042                Dti_total=XmlDataRead.Time(end)-XmlDataRead.Time(1);%total time interval covered by the time vector
     2043                missing_indices=sizDti+1+ind_start:size(i1_series,3)+1;% remaining set of frame indices for which time needs to be found
     2044                repeat_nbre=1+floor((missing_indices-sizDti-ind_start)/(sizDti-1));% number of repetitions of Dti
     2045                time_indices=1+mod(missing_indices-sizDti-1,sizDti-1);
     2046                for j=1:size(XmlDataRead.Time,2)
     2047                    XmlData.Time(missing_indices,j)=XmlDataRead.Time(time_indices,j)+repeat_nbre'*Dti_total;
     2048                end
     2049                % update the xml file with NbDti
     2050                t=xmltree(XmlFileName);
     2051                uid_NbDti=find(t,'ImaDoc/Camera/BurstTiming/NbDti')
     2052                if isempty(uid_NbDti)
     2053                    uid_BurstTiming=find(t,'ImaDoc/Camera/BurstTiming')
     2054                    [t,uid_NbDti]=add(t,uid_BurstTiming,'element','NbDti');
     2055                end
     2056                [t,uid_NbDti]=add(t,uid_NbDti,'chardata',num2str(repeat_nbre(end)-1));
     2057                save(t,XmlFileName)
     2058            end
    20252059        end
    20262060        set(handles.view_xml,'BackgroundColor',[1 1 1])% paint back to white
    … …  
    20702104
    20712105%% store last index in handles.MaxIndex_i and .MaxIndex_j
    2072 nbfield=max(max(max(i2_series)));
    2073 if isempty(nbfield)
    2074     nbfield=max(max(max(i1_series)));
    2075 end
    2076 nbfield_j=max(max(max(j2_series)));
    2077 if isempty(nbfield_j)
    2078     nbfield_j=max(max(max(j1_series)));
    2079 end
    20802106if isfield(XmlData,'Time')&& ~isempty(XmlData.Time)
    20812107    %transform .Time to a column vector if it is a line vector the nomenclature uses a single index
Note: See TracChangeset for help on using the changeset viewer.