source: trunk/src/series/check_peaklock.m @ 846

Last change on this file since 846 was 810, checked in by g7moreau, 10 years ago
  • Add license
File size: 26.7 KB
RevLine 
[590]1% 'check_peaklocking': estimte peaklocking error in a civ field series TODO: UPDATE
2%------------------------------------------------------------------------
3% function ParamOut=check_peaklocking(Param)
4%
5%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
6%
7% This function is used in four modes by the GUI series:
8%           1) config GUI: with no input argument, the function determine the suitable GUI configuration
9%           2) interactive input: the function is used to interactively introduce input parameters, and then stops
10%           3) RUN: the function itself runs, when an appropriate input  structure Param has been introduced.
11%           4) BATCH: the function itself proceeds in BATCH mode, using an xml file 'Param' as input.
12%
13% This function is used in four modes by the GUI series:
14%           1) config GUI: with no input argument, the function determine the suitable GUI configuration
15%           2) interactive input: the function is used to interactively introduce input parameters, and then stops
16%           3) RUN: the function itself runs, when an appropriate input  structure Param has been introduced.
17%           4) BATCH: the function itself proceeds in BATCH mode, using an xml file 'Param' as input.
18%
19%OUTPUT
20% GUI_input=list of options in the GUI series.fig needed for the function
21%
22%INPUT:
23% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
24% In batch mode, Param is the name of the corresponding xml file containing the same information
25% In the absence of input (as activated when the current Action is selected
26% in series), the function ouput GUI_input set the activation of the needed GUI elements
27%
28% Param contains the elements:(use the menu bar command 'export/GUI config' in series to see the current structure Param)
29%    .InputTable: cell of input file names, (several lines for multiple input)
30%                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}
31%    .OutputSubDir: name of the subdirectory for data outputs
32%    .OutputDirExt: directory extension for data outputs
33%    .Action: .ActionName: name of the current activated function
34%             .ActionPath:   path of the current activated function
35%    .IndexRange: set the file or frame indices on which the action must be performed
36%    .FieldTransform: .TransformName: name of the selected transform function
37%                     .TransformPath:   path  of the selected transform function
38%                     .TransformHandle: corresponding function handle
39%    .InputFields: sub structure describing the input fields withfields
40%              .FieldName: name of the field
41%              .VelType: velocity type
42%              .FieldName_1: name of the second field in case of two input series
43%              .VelType_1: velocity type of the second field in case of two input series
44%    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
45%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[810]46
47%=======================================================================
48% Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France
49%   http://www.legi.grenoble-inp.fr
50%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
51%
52%     This file is part of the toolbox UVMAT.
53%
54%     UVMAT is free software; you can redistribute it and/or modify
55%     it under the terms of the GNU General Public License as published
56%     by the Free Software Foundation; either version 2 of the license,
57%     or (at your option) any later version.
58%
59%     UVMAT is distributed in the hope that it will be useful,
60%     but WITHOUT ANY WARRANTY; without even the implied warranty of
61%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
62%     GNU General Public License (see LICENSE.txt) for more details.
63%=======================================================================
64
[590]65function  ParamOut=check_peaklocking(Param)
66
67%% set the input elements needed on the GUI series when the action is selected in the menu ActionName
68if ~exist('Param','var') % case with no input parameter
69    ParamOut={'AllowInputSort';'off';...% allow alphabetic sorting of the list of input files (options 'off'/'on', 'off' by default)
70        'WholeIndexRange';'off';...% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
71        'NbSlice';'on'; ...%nbre of slices ('off' by default)
72        'VelType';'two';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
73        'FieldName';'off';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
74        'FieldTransform'; 'off';...%can use a transform function
75        'ProjObject';'on';...%can use projection object(option 'off'/'on',
76        'Mask';'off';...%can use mask option   (option 'off'/'on', 'off' by default)
77        'OutputDirExt';'.pklock';...%set the output dir extension
78               ''};
79        return
80end
81
82%%%%%%%%%%%% STANDARD PART  %%%%%%%%%%%%
83%% select different modes,  RUN, parameter input, BATCH
84% BATCH  case: read the xml file for batch case
85if ischar(Param)
86        Param=xml2struct(Param);
87        checkrun=0;
88% RUN case: parameters introduced as the input structure Param
89else
90    hseries=guidata(Param.hseries);%handles of the GUI series
91    if isfield(Param,'Specific')&& strcmp(Param.Specific,'?')
92        checkrun=1;% will only search interactive input parameters (preparation of BATCH mode)
93    else
94        checkrun=2; % indicate the RUN option is used
95    end
96end
97ParamOut=Param; %default output
98OutputDir=[Param.OutputSubDir Param.OutputDirExt];
99
100%% root input file(s) and type
101RootPath=Param.InputTable(:,1);
102RootFile=Param.InputTable(:,3);
103SubDir=Param.InputTable(:,2);
104NomType=Param.InputTable(:,4);
105FileExt=Param.InputTable(:,5);
106[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
107%%%%%%%%%%%%
108% The cell array filecell is the list of input file names, while
109% filecell{iview,fileindex}:
110%        iview: line in the table corresponding to a given file series
111%        fileindex: file index within  the file series,
112% 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
113% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
114%%%%%%%%%%%%
115NbSlice=1;%default
116if isfield(Param.IndexRange,'NbSlice')&&~isempty(Param.IndexRange.NbSlice)
117    NbSlice=Param.IndexRange.NbSlice;
118end
119nbview=1;%number of input file series (lines in InputTable)
120nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
121nbfield_i=size(i1_series{1},2); %nb of fields for the i index
122nbfield=nbfield_j*nbfield_i; %total number of fields
123nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices)
124nbfield=nbfield_i*NbSlice; %total number of fields after adjustement
125
126%determine the file type on each line from the first input file
127ImageTypeOptions={'image','multimage','mmreader','video'};
128NcTypeOptions={'netcdf','civx','civdata'};
129for iview=1:nbview
130    if ~exist(filecell{iview,1}','file')
[668]131        disp_uvmat('ERROR',['the first input file ' filecell{iview,1} ' does not exist'],checkrun)
[590]132        return
133    end
[784]134    [FileInfo{iview},MovieObject{iview}]=get_file_info(filecell{iview,1});
[783]135    FileType{iview}=FileInfo{iview}.FileType;
[590]136    CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images
137    CheckNc{iview}=~isempty(find(strcmp(FileType{iview},NcTypeOptions)));% =1 for netcdf files
138    if ~isempty(j1_series{iview})
139        frame_index{iview}=j1_series{iview};
140    else
141        frame_index{iview}=i1_series{iview};
142    end
143end
144
145%% calibration data and timing: read the ImaDoc files
146%none
147
148%% coordinate transform or other user defined transform
149% none
150
151%%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
152 % EDIT FROM HERE
153
154%% check the validity of  ctinput file types
155%none
156
157%% Set field names and velocity types
158InputFields{1}=[];%default (case of images)
159if isfield(Param,'InputFields')
160    InputFields{1}=Param.InputFields;
161end
162% only one input fieldseries
163
164%% Initiate output fields
165%initiate the output structure as a copy of the first input one (reproduce fields)
166[DataOut,tild,errormsg] = read_field(filecell{1,1},FileType{1},InputFields{1},1);
167if ~isempty(errormsg)
[668]168    disp_uvmat('ERROR',['error reading ' filecell{1,1} ': ' errormsg],checkrun)
[590]169    return
170end
171time_1=[];
172if isfield(DataOut,'Time')
173    time_1=DataOut.Time(1);
174end
175if CheckNc{iview}
176    if isempty(strcmp('Conventions',DataOut.ListGlobalAttribute))
177        DataOut.ListGlobalAttribute=['Conventions' DataOut.ListGlobalAttribute];
178    end
179    DataOut.Conventions='uvmat';
180    DataOut.ListGlobalAttribute=[DataOut.ListGlobalAttribute {Param.Action}];
181    ActionKey='Action';
182    while isfield(DataOut,ActionKey)
183        ActionKey=[ActionKey '_1'];
184    end
185    DataOut.(ActionKey)=Param.Action;
186    DataOut.ListGlobalAttribute=[DataOut.ListGlobalAttribute {ActionKey}];
187    if isfield(DataOut,'Time')
188        DataOut.ListGlobalAttribute=[DataOut.ListGlobalAttribute {'Time','Time_end'}];
189    end
190end
191
192%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
193index_slice=1:nbfield;% select the file indices
194for index=index_slice
195    if checkrun
196        update_waitbar(hseries.Waitbar,index/(nbfield))
197        stopstate=get(hseries.RUN,'BusyAction');
198    else
199        stopstate='queue';
200    end
201    if isequal(stopstate,'queue')% enable STOP command
202        Data=cell(1,nbview);%initiate the set Data;
203        nbtime=0;
204        dt=[];
205        %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
206        for iview=1:nbview
207            % reading input file(s)
208            [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
209            if ~isempty(errormsg)
210                errormsg=['time_series / read_field / ' errormsg];
211                display(errormsg)
212                break
213            end
214        end
215        if isempty(errormsg)
216            Field=Data{1}; % default input field structure
217            % coordinate transform (or other user defined transform)
218            % none
219           
220            %field projection on an object
221            if Param.CheckObject
222                [Field,errormsg]=proj_field(Field,Param.ProjObject);
223                if ~isempty(errormsg)
224                    msgbox_uvmat('ERROR',['time_series / proj_field / ' errormsg])
225                    return
226                end
227            end
228           
229            % initiate the time series at the first iteration
[592]230            if index==1
[590]231                % stop program if the first field reading is in error
232                if ~isempty(errormsg)
[668]233                    disp_uvmat('ERROR',['time_series / sub_field / ' errormsg],checkrun)
[590]234                    return
235                end
236                DataOut=Field;%default
237                DataOut.NbDim=Field.NbDim+1; %add the time dimension for plots
238                nbvar=length(Field.ListVarName);
239                if nbvar==0
[668]240                    disp_uvmat('ERROR','no input variable selected',checkrun)
[590]241                    return
242                end
243                testsum=2*ones(1,nbvar);%initiate flag for action on each variable
244                if isfield(Field,'VarAttribute') % look for coordinate and flag variables
245                    for ivar=1:nbvar
246                        if length(Field.VarAttribute)>=ivar && isfield(Field.VarAttribute{ivar},'Role')
247                            var_role=Field.VarAttribute{ivar}.Role;%'role' of the variable
248                            if isequal(var_role,'errorflag')
[668]249                                disp_uvmat('ERROR','do not handle error flags in time series',checkrun)
[590]250                                return
251                            end
252                            if isequal(var_role,'warnflag')
253                                testsum(ivar)=0;  % not recorded variable
254                                eval(['DataOut=rmfield(DataOut,''' Field.ListVarName{ivar} ''');']);%remove variable
255                            end
256                            if isequal(var_role,'coord_x')| isequal(var_role,'coord_y')|...
257                                    isequal(var_role,'coord_z')|isequal(var_role,'coord')
258                                testsum(ivar)=1; %constant coordinates, record without time evolution
259                            end
260                        end
261                        % check whether the variable ivar is a dimension variable
262                        DimCell=Field.VarDimName{ivar};
263                        if ischar(DimCell)
264                            DimCell={DimCell};
265                        end
266                        if numel(DimCell)==1 && isequal(Field.ListVarName{ivar},DimCell{1})%detect dimension variables
267                            testsum(ivar)=1;
268                        end
269                    end
270                end
271                for ivar=1:nbvar
272                    if testsum(ivar)==2
273                        eval(['DataOut.' Field.ListVarName{ivar} '=[];'])
274                    end
275                end
276                DataOut.ListVarName=[{'Time'} DataOut.ListVarName];
277            end
278           
279            % add data to the current field
280            for ivar=1:length(Field.ListVarName)
281                VarName=Field.ListVarName{ivar};
282                VarVal=Field.(VarName);
283                if testsum(ivar)==2% test for recorded variable
284                    if isempty(errormsg)
285                        if isequal(Param.ProjObject.ProjMode,'inside')% take the average in the domain for 'inside' mode
286                            if isempty(VarVal)
[668]287                                disp_uvmat('ERROR',['empty result at frame index ' num2str(i1_series{iview}(index))],checkrun)
[590]288                                return
289                            end
290                            VarVal=mean(VarVal,1);
291                        end
292                        VarVal=shiftdim(VarVal,-1); %shift dimension
293                        DataOut.(VarName)=cat(1,DataOut.(VarName),VarVal);%concanete the current field to the time series
294                    else
295                        DataOut.(VarName)=cat(1,DataOut.(VarName),0);% put each variable to 0 in case of input reading error
296                    end
297                elseif testsum(ivar)==1% variable representing fixed coordinates
298                    VarInit=DataOut.(VarName);
299                    if isempty(errormsg) && ~isequal(VarVal,VarInit)
[668]300                        disp_uvmat('ERROR',['time series requires constant coordinates ' VarName],checkrun)
[590]301                        return
302                    end
303                end
304            end
[592]305
[590]306        end
307    end
308end
309%%%%%%% END OF LOOP WITHIN A SLICE
310
311%remove time for global attributes if exists
312Time_index=find(strcmp('Time',DataOut.ListGlobalAttribute));
313if ~isempty(Time_index)
314    DataOut.ListGlobalAttribute(Time_index)=[];
315end
316DataOut.Conventions='uvmat';
317for ivar=1:numel(DataOut.ListVarName)
318    VarName=DataOut.ListVarName{ivar};
319    eval(['DataOut.' VarName '=squeeze(DataOut.' VarName ');']) %remove singletons
320end
321
322% add time dimension
323for ivar=1:length(Field.ListVarName)
324    DimCell=Field.VarDimName(ivar);
325    if testsum(ivar)==2%variable used as time series
326        DataOut.VarDimName{ivar}=[{'Time'} DimCell];
327    elseif testsum(ivar)==1
328        DataOut.VarDimName{ivar}=DimCell;
329    end
330end
331indexremove=find(~testsum);
332if ~isempty(indexremove)
333    DataOut.ListVarName(1+indexremove)=[];
334    DataOut.VarDimName(indexremove)=[];
335    if isfield(DataOut,'Role') && ~isempty(DataOut.Role{1})%generaliser aus autres attributs
336        DataOut.Role(1+indexremove)=[];
337    end
338end
339
340%shift variable attributes
341if isfield(DataOut,'VarAttribute')
342    DataOut.VarAttribute=[{[]} DataOut.VarAttribute];
343end
344DataOut.VarDimName=[{'Time'} DataOut.VarDimName];
345DataOut.Action=Param.Action;%name of the processing programme
346test_time=diff(DataOut.Time)>0;% test that the readed time is increasing (not constant)
347if ~test_time
348    DataOut.Time=1:filecounter;
349end
350
351% display nbmissing
352if ~isequal(nbmissing,0)
[668]353    disp_uvmat('WARNING',[num2str(nbmissing) ' files skipped: missing files or bad input, see command window display'],checkrun)
[590]354end
355
356%name of result file
357OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(1),i1_series{1}(end),i_slice,[]);
358errormsg=struct2nc(OutputFile,DataOut); %save result file
359if isempty(errormsg)
360    display([OutputFile ' written'])
361else
[668]362    disp_uvmat('ERROR',['error in Series/struct2nc: ' errormsg],checkrun)
[590]363end
364
365return
366
367%%%%%%%%%%%%%%%%%%  END%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[583]368%evaluation of peacklocking errors
369%use splinhist: give spline coeff cc for a smooth histo (call spline4)
370%use histsmooth(x,cc): calculate the smooth histo for any value x
371%use histder(x,cc): calculate the derivative of the smooth histo
372global hfig1 hfig2 hfig3
373global nbb Uval Vval Uhist Vhist % nbb resolution of the histogram nbb=10: 10 values in unity interval
374global xval xerror yval yerror
375
376set(handles.vector_y,'Value',1)% trigger the option Uhist on the interface
377set(handles.Vhist_input,'Value',1)
378set(handles.cm_switch,'Value',0) % put the switch to 'pixel'
379
380%adjust the extremal values of the histogram in U with respect to integer
381%values
382minimU=round(min(Uval)-0.5)+0.5; %first value of the histogram with integer bins
383maximU=round(max(Uval)-0.5)+0.5;
384minim_fin=(minimU-0.5+1/(2*nbb)); % first bin valueat the beginning of an integer interval
385maxim_fin=(maximU+0.5-1/(2*nbb)); % last integer value
386nb_bin_min= round(-(minim_fin - min(Uval))*nbb); % nbre of bins added below
387nb_bin_max=round((maxim_fin -max(Uval))*nbb); %nbre of bins added above
388Uval=[minim_fin:(1/nbb):maxim_fin];
389histu_min=zeros(nb_bin_min,1);
390histu_max=zeros(nb_bin_max,1);
391Uhist=[histu_min; Uhist ;histu_max]; % column vector
392
393%adjust the extremal values of the histogram in V
394minimV=round(min(Vval-0.5)+0.5);
395maximV=round(max(Vval-0.5)+0.5);
396minim_fin=minimV-0.5+1/(2*nbb); % first bin valueat the beginning of an integer interval
397maxim_fin=maximV+0.5-1/(2*nbb); % last integer value
398nb_bin_min=round((min(Vval) - minim_fin)*nbb); % nbre of bins added below
399nb_bin_max=round((maxim_fin -max(Vval))*nbb);
400Vval=[minim_fin:(1/nbb):maxim_fin];
401histu_min=zeros(nb_bin_min,1);
402histu_max=zeros(nb_bin_max,1);
403Vhist=[histu_min; Vhist ;histu_max]; % column vector
404
405% RUN_histo_Callback(hObject, eventdata, handles)
406% %adjust the histogram to integer values:
407
408%histoU and V
409[Uhistinter,xval,xerror]=peaklock(nbb,minimU,maximU,Uhist);
410[Vhistinter,yval,yerror]=peaklock(nbb,minimV,maximV,Vhist);
411
412% selection of value ranges such that histo>=10 (enough statistics)
413Uval_ind=find(Uhist>=10);
414ind_min=min(Uval_ind);
415ind_max=max(Uval_ind);
416U_min=Uval(ind_min);% minimum allowed value
417U_max=Uval(ind_max);%maximum allowed value
418
419% selection of value ranges such that histo>=10 (enough statistics)
420Vval_ind=find(Vhist>=10);
421ind_min=min(Vval_ind);
422ind_max=max(Vval_ind);
423V_min=Vval(ind_min);% minimum allowed value
424V_max=Vval(ind_max);%maximum allowed value
425
426figure(4)% plot U histogram with smoothed one
427plot(Uval,Uhist,'b')
428grid on
429hold on
430plot(Uval,Uhistinter,'r');
431hold off
432
433figure(5)% plot V histogram with smoothed one
434plot(Vval,Vhist,'b')
435grid on
436hold on
437plot(Vval,Vhistinter,'r');
438hold off
439
440figure(6)% plot pixel error in two subplots
441hfig4=subplot(2,1,1);
442hfig5=subplot(2,1,2);
443axes(hfig4)
444plot(xval,xerror)
445axis([U_min U_max -0.4 0.4])
446xlabel('velocity u (pix)')
447ylabel('peaklocking error (pix)')
448grid on
449axes(hfig5)
450plot(yval,yerror)
451axis([V_min V_max -0.4 0.4]);
452xlabel('velocity v (pix)')
453ylabel('peaklocking error (pix)')
454grid on
455
456
457
458
459
460
461
462
463
464%'peaklock': determines peacklocking errors from velocity histograms.
465%-------------------------------------------------------
466%first smooth the input histogram 'histu' in such a way that the integral over
467%n-n+1 is preserved, then deduce the peaklocking 'error' function of the pixcel displacement 'x'.
468%
469% [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
470%OUTPUT:
471%histinter: smoothed interpolated histogram
472% x: vector of displacement values.
473% error: vector of estimated errors corresponding to x
474%INPUT:
475%histu=vector representing the values of histogram  of measured velocity ;
476%minim, maxim: extremal values of the measured velocity (absica for histu)
477%nbb: number of bins inside each integer interval for the histograms
478%SUBROUTINES INCLUDED:
479%spline4.m% spline interpolation at 4th order
480%splinhist.m: give spline coeff cc for a smooth histo (call spline4)
481%histsmooth.m(x,cc): calculate the smooth histo for any value x
482%histder.m(x,cc): calculate the derivative of the smooth histo
483function [histinter,x,error]=peaklock(nbb,minim,maxim,histu)
484
485nint=maxim-minim+1
486xfin=[minim-0.5+1/(2*nbb):(1/nbb):maxim+0.5-(1/(2*nbb))];
487histo=(reshape(histu,nbb,nint));%extract values with x between integer -1/2 integer +1/2
488Integ=sum(histo)/nbb; %integral of the pdf on each integer bin
489[histinter,cc]=splinhist(Integ,minim,nbb);
490histx=reshape(histinter,nbb,nint);
491xint=[minim:1:maxim];
492x=zeros(nbb,nint);
493%determination of the displacement x(j,:)
494%j=1
495delx=histo(1,:)./histsmooth(-0.5*ones(1,nint),cc)/nbb;
496%del(1,:)=delx;
497x(1,:)=-0.5+delx-(delx.*delx/2).*histder(-0.5*ones(1,nint),cc);
498%histx(1,:)=histsmooth(x(j-1,:),cc);
499for j=2:nbb
500    delx=histo(j,:)./histsmooth(x(j-1,:),cc)/nbb;
501    %delx=delx.*(delx<3*ones(1,nint)/nbb)+3*ones(1,nint)/nbb.*~(delx <3*ones(1,nint)/nbb)
502    x(j,:)=x(j-1,:)+delx-(delx.*delx/2).*histder(x(j-1,:),cc);
503end
504%reshape
505xint=ones(nbb,1)*xint;
506x=x+xint;
507x=reshape(x,1,nbb*nint);
508error=xfin+1/(2*nbb)-x;
509
510%-------------------------------------------------------
511% --- determine the spline coefficients cc for the interpolated histogram.
512%-------------------------------------------------
513function [histsmooth,cc]= splinhist(Integ,mini,nbb)
514% provides a smooth histogramm histmooth, which remains always positive,
515% and is such that its sum over each integer bin [i-1/2 i+1/2] is equal to
516% Integ(i). The function determines histmooth as the exponential of a 4th
517% order spline function and adjust the cefficients by a Newton method to
518% fit the integral conditions Integ
519% histmooth is determined at the abscissa
520% xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)
521%cc(1-5,i) provides the spline coefficients
522
523% order 0
524siz=size(Integ);
525nint=siz(2);
526izero=find(Integ==0); %indices of zero elements
527inonzero=find(Integ);
528Integ(izero)=min(Integ(inonzero));
529aa=log(Integ);%initial guess for a coeff
530spli=spline4(aa,mini,nbb);  %appel à la fonction spline4
531histsmooth=exp(spli);
532
533S=(sum(reshape(histsmooth,nbb,nint)))/nbb;% integral of the fit histsmooth on ]i-1/2 i+1/2[
534epsilon=max(abs(Integ-S));
535iter=0;
536while epsilon > 0.000001 & iter<10
537ident=eye(nint);
538dSda=ones(nint);
539for j=1:nint% determination of the jacobian matrix dSda
540dhistda=spline4(ident(j,:),mini,nbb);
541expdhistda=dhistda.*histsmooth;
542dSda(j,:)=(sum(reshape(expdhistda,nbb,nint)))/nbb;
543end
544aa=aa+(Integ-S)*inv(dSda);%new estimate of coefficients aa by linear interpolation
545[spli,bb]=spline4(aa,mini,nbb);% new fit histsmooth
546histsmooth=exp(spli);
547S=(sum(reshape(histsmooth,nbb,nint)))/nbb;% integral of the fit histsmooth on ]i-1/2 i+1/2[
548epsilon=max(abs(Integ-S));
549iter=iter+1;
550end
551if iter==10, errordlg('splinhist did not converge after 10 iterations'),end
552cc(1,:)=aa;
553cc(2,:)=bb(1,:);
554cc(3,:)=bb(2,:);
555cc(4,:)=bb(3,:);
556cc(5,:)=bb(4,:);
557
558%-------------------------------------------------------
559% --- determine the 4th order spline coefficients from the function values aa.
560%-------------------------------------------------
561function [histsmooth,bb]= spline4(aa,mini,n)
562% spline interpolation at 4th order
563%aa=vector of values of a function at integer abscissa, starting at mini
564%n=number of subdivisions for the interpolated function
565% histmooth =interpolated values at absissa
566% xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))] (maxi=mini+size(aa)-1)
567%bb=[b(i);c(i);d(i); e(i)] matrix of spline coeff
568L1=[1/2 1/4 1/8 1/16;1 1 3/4 1/2;0 2 3 3;0 0 6 12];
569L2=[-1/2 1/4 -1/8 1/16;1 -1 3/4 -1/2;0 2 -3 3;0 0 6 -12];
570M=inv(L2)*L1;
571[V,D]=eig(M);
572F=-inv(V)*inv(L2)*[1 ;0 ;0;0];
573a1rev=[1 -1/D(1,1)];
574b1rev=[F(1)/D(1,1)];
575a2rev=[1 -1/D(2,2)];
576b2rev=[F(2)/D(2,2)];
577a3=[1 -D(3,3)];
578b3=[F(3)];
579a4=[1 -D(4,4)];
580b4=[F(4)];
581
582%data
583% n=10;% résolution de la pdf: nbre de points par unite de u
584% mini=-10.0;%general mini=uint16(min(values)-1 CHOOSE maxi-mini+1 EVEN
585% maxi=9.0; % general maxi=uint16(max(values))+1
586%nint=double(maxi-mini+1); % nombre d'intervals entiers EVEN!
587siz=size(aa);
588nint=siz(2);
589maxi=mini+nint-1;
590npdf=nint*n;% nbre total d'intervals à introduire dans la pdf: hist(u,npdf)
591%simulation de pdf
592xfin=[mini-0.5+1/(2*n):(1/n):maxi+0.5-(1/(2*n))];% valeurs d'interpolation: we take n values in each integer interval
593%histolin=exp(-(xfin-1).*(xfin-1)).*(2+cos(10*(xfin-1)));% simulation d'une pdf
594%histo=log(histolin);
595%histo=sin(2*pi*xfin);
596%histextract=(reshape(histo,n,nint));
597%aa=sum(histextract)/n %integral of the pdf on each integer bin
598IP=[0 diff(aa)];
599Irev=zeros(size(aa));
600for i=1:nint
601    Irev(i)=aa(end-i+1);
602end
603IPrev=[0 diff(Irev)];
604
605%get the spline coelfficients a_d, using filter on the eigen vectors A,B,C
606Arev=filter(b1rev,a1rev,IPrev);
607Brev=filter(b2rev,a2rev,IPrev);
608C=filter(b3,a3,IP);
609D=filter(b4,a4,IP);
610A=zeros(size(Arev));
611B=zeros(size(Brev));
612for i=1:nint
613    A(i)=Arev(end-i+1);
614    B(i)=Brev(end-i+1);
615end
616%Matr=V*[A;B;C;D];
617bb=V*[A;B;C;D];
618%b=Matr(1,:);
619%c=Matr(2,:);
620%d=Matr(3,:);
621%e=Matr(4,:);
622%a=aa;
623
624%calculate the interpolation using the spline coefficients a-d
625%xextract=(reshape(xfin,n,nint));%
626chi=xfin+1/(2*n)-min(xfin)-double(int16(xfin+(1/(2*n))-min(xfin)))-0.5;% decimal part
627chi2=chi.*chi;
628chi3=chi2.*chi;
629chi4=chi3.*chi;
630avec=reshape(ones(n,1)*aa,1,n*nint);
631bvec=reshape(ones(n,1)*bb(1,:),1,n*nint);
632cvec=reshape(ones(n,1)*bb(2,:),1,n*nint);
633dvec=reshape(ones(n,1)*bb(3,:),1,n*nint);
634evec=reshape(ones(n,1)*bb(4,:),1,n*nint);
635histsmooth=avec+bvec.*chi+cvec.*chi2+dvec.*chi3+evec.*chi4;
636
637%-------------------------------------------------------
638% --- determine the interpolated histogram at points chi from the spline ceff cc.
639%-------------------------------------------------
640function histx= histsmooth(chi,cc)
641% provides the value of the interpolated histogram at values chi=x-i
642%(difference with the mnearest integer)
643% cc(5,size(chi)) is the set of spline coefficients obtained by splinhist
644chi2=chi.*chi;
645chi3=chi2.*chi;
646chi4=chi3.*chi;
647histx=exp(cc(1,:)+cc(2,:).*chi+cc(3,:).*chi2+cc(4,:).*chi3+cc(5,:).*chi4);
648
649%-------------------------------------------------------
650% --- determine the derivative p'/p of the interpolated histogram at points chi from the spline ceff cc.
651%-------------------------------------------------
652function histder= histder(chi,cc)
653% provides the logarithmique derivative p'/p of the interpolated histogram
654%at values chi=x-i
655%(difference with the nearest integer)
656% cc(5,size(chi)) is the set of spline coefficients obtained by splinhist
657chi2=chi.*chi;
658chi3=chi2.*chi;
659chi4=chi3.*chi;
[810]660histder=cc(2,:)+2*cc(3,:).*chi+3*cc(4,:).*chi2+4*cc(5,:).*chi3;
Note: See TracBrowser for help on using the repository browser.