Changeset 246 for trunk/src/civ_uvmat.m


Ignore:
Timestamp:
Apr 28, 2011, 10:52:31 AM (14 years ago)
Author:
sommeria
Message:

thin plate shell (patch) introduced

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/civ_uvmat.m

    r243 r246  
    11% To develop....
    22function [Data,errormsg]= civ_uvmat(Param,ncfile)
     3errormsg='';
    34Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
    45Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    56Data.Program='civ_uvmat';
    67Data.CivStage=0;%default
     8ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};
     9ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'};
    710
    811%% Civ1
    … …  
    2023    shiftx=str2num(par_civ1.shiftx);
    2124    shifty=str2num(par_civ1.shifty);
    22     miniy=max(1+isy2-shifty,1+iby2);
     25    miniy=max(1+isy2+shifty,1+iby2);
    2326    minix=max(1+isx2-shiftx,1+ibx2);
    24     maxiy=min(size(image1,1)-isy2-shifty,size(image1,1)-iby2);
     27    maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2);
    2528    maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2);
    2629    [GridX,GridY]=meshgrid(minix:stepx:maxix,miniy:stepy:maxiy);
    … …  
    2932   
    3033    % caluclate velocity data (y and v in indices, reverse to y component)
    31     [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty,PointCoord,str2num(par_civ1.rho), []);
     34    [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho), []);
    3235    list_param=(fieldnames(par_civ1))';
    3336    list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'};
    … …  
    6770    Data.CivStage=1;
    6871else
    69     Data=nc2struct(ncfile)%read existing netcdf file
     72    if isfield(Param,'Fix1')
     73        Data=nc2struct(ncfile,ListVarCiv1);%read civ1 data in the existing netcdf file
     74    else
     75        Data=nc2struct(ncfile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file
     76    end
     77    if isfield(Data,'Txt')
     78        msgbox_uvmat('ERROR',Data.Txt)
     79        return
     80    end
     81    % read Civx data
    7082    if isfield(Data,'absolut_time_T0')%read civx file
    7183        var={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF';'vec_X','vec_Y','vec_U','vec_V','vec_C','vec_F','vec_FixFlag'};
    … …  
    95107        ParamName=ListFixParam{ilist};
    96108        ListName=['Fix1_' ParamName];
    97         ['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']
    98109        eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
    99110        eval(['Data.' ListName '=Param.Fix1.' ParamName ';'])
    … …  
    107118    Data.VarDimName=[Data.VarDimName {'nbvec1'}];
    108119    nbvar=length(Data.ListVarName);
    109     Data.VarAttribute{nbvar}.Role='errorflag';
    110     [Data.Civ1_FF]=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
     120    Data.VarAttribute{nbvar}.Role='errorflag';   
     121    Data.Civ1_FF=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
    111122    Data.CivStage=2;                               
    112123end   
    … …  
    117128    Data.Patch1_Threshold=str2double(Param.Patch1.Threshold);
    118129    Data.Patch1_SubDomain=str2double(Param.Patch1.SubDomain);
    119     Data.ListVarName=[Data.ListVarName {'Patch1_U','Patch1_V'}];
    120     Data.VarDimName=[Data.VarDimName {'nbvec1','nbvec1'}];
     130    Data.ListVarName=[Data.ListVarName {'Civ1_U_Diff','Civ1_V_Diff','Civ1_X_SubRange','Civ1_Y_SubRange','Civ1_NbSites','Civ1_X_tps','Civ1_Y_tps','Civ1_U_tps','Civ1_V_tps'}];
     131    Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'NbSubDomain1','Two'},{'NbSubDomain1','Two'},'NbSubDomain1',...
     132             {'NbVec1Sub','NbSubDomain1'},{'NbVec1Sub','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}];
    121133    nbvar=length(Data.ListVarName);
    122134    Data.VarAttribute{nbvar-1}.Role='vector_x';
    123135    Data.VarAttribute{nbvar}.Role='vector_y';
    124     Data.Patch1_U=zeros(size(Data.Civ1_X));
    125     Data.Patch1_V=zeros(size(Data.Civ1_X));
     136    Data.Civ1_U_Diff=zeros(size(Data.Civ1_X));
     137    Data.Civ1_V_Diff=zeros(size(Data.Civ1_X));
    126138    if isfield(Data,'Civ1_FF')
    127139        ind_good=find(Data.Civ1_FF==0);
    … …  
    129141        ind_good=1:numel(Data.Civ1_X);
    130142    end
    131     Data.Civ1_X
    132     [Ures, Vres,FFres]=...
     143    [Data.Civ1_X_SubRange,Data.Civ1_Y_SubRange,Data.Civ1_NbSites,FFres,Ures, Vres,Data.Civ1_X_tps,Data.Civ1_Y_tps,Data.Civ1_U_tps,Data.Civ1_V_tps]=...
    133144                            patch_uvmat(Data.Civ1_X(ind_good)',Data.Civ1_Y(ind_good)',Data.Civ1_U(ind_good)',Data.Civ1_V(ind_good)',Data.Patch1_Rho,Data.Patch1_Threshold,Data.Patch1_SubDomain);
    134                         size(Ures)
    135                         size(Vres)
    136                         size(FFres)
    137                         size(ind_good)
    138       Data.Patch1_U(ind_good)=Ures;
    139       Data.Patch1_V(ind_good)=Vres;
    140       Data.Civ1_FF(ind_good)=FFres
    141       Data.CivStage=3;                               
     145      Data.Civ1_U_Diff(ind_good)=Data.Civ1_U(ind_good)-Ures;
     146      Data.Civ1_V_Diff(ind_good)=Data.Civ1_V(ind_good)-Vres;
     147      Data.Civ1_FF(ind_good)=FFres;
     148      Data.CivStage=3;                             
    142149end   
    143150%% write result
     151% 'TESTcalc'
     152% [DataOut,errormsg]=calc_field('velocity',Data)
     153if exist('ncfile','var')
    144154errormsg=struct2nc(ncfile,Data);
    145 
     155end
    146156
    147157
    … …  
    192202    FF=FF==1 | (U.*U+V.*V)>thresh;
    193203end
    194  
     204FF=double(FF);
    195205%
    196206% % criterium on velocity values
    … …  
    270280%------------------------------------------------------------------------
    271281% patch function
    272 function [U_patch,V_patch,FF,SubRangx,SubRangy,X_ctrs,Y_ctrs] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)
     282function [SubRangx,SubRangy,nbpoints,FF,U_smooth,V_smooth,X_tps,Y_tps,U_tps,V_tps] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain)
    273283%subdomain decomposition
    274284warning off
    … …  
    278288Y=reshape(Y,[],1);
    279289nbvec=numel(X);
    280 NbSubDomain=ceil(nbvec/SubDomain)
    281 MinX=min(X)
    282 MinY=min(Y)
    283 MaxX=max(X)
    284 MaxY=max(Y)
     290NbSubDomain=ceil(nbvec/SubDomain);
     291MinX=min(X);
     292MinY=min(Y);
     293MaxX=max(X);
     294MaxY=max(Y);
    285295RangX=MaxX-MinX;
    286296RangY=MaxY-MinY;
    287 AspectRatio=RangY/RangX
    288 NbSubDomainX=ceil(sqrt(NbSubDomain/AspectRatio))
    289 NbSubDomainY=ceil(sqrt(NbSubDomain*AspectRatio))
    290 
     297AspectRatio=RangY/RangX;
     298NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);
     299NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);
     300NbSubDomain=NbSubDomainX*NbSubDomainY;
    291301SizX=RangX/NbSubDomainX;%width of subdomains
    292302SizY=RangY/NbSubDomainY;%height of subdomains
    293303CentreX=linspace(MinX+SizX/2,MaxX-SizX/2,NbSubDomainX);
    294304CentreY=linspace(MinY+SizY/2,MaxY-SizY/2,NbSubDomainY);
    295 SubIndexX=ceil((X-MinX)/SizX);%subdomain index of vectors
    296 SubIndexY=ceil((Y-MinY)/SizY);
    297 rho=RangX*RangY*Rho;%optimum rho increase as teh area of the subdomain
    298 U_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    299 V_tps=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    300 X_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    301 Y_dist=zeros(length(X),NbSubDomainY,NbSubDomainX);%default spline
    302 U_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
    303 V_smooth=zeros(length(X),NbSubDomainY,NbSubDomainX);
    304 dist_ctre=zeros(length(X),NbSubDomainY,NbSubDomainX);
     305[CentreX,CentreY]=meshgrid(CentreX,CentreY);
     306CentreY=reshape(CentreY,1,[]);
     307CentreX=reshape(CentreX,1,[]);
     308rho=SizX*SizY*Rho/1000000;%optimum rho increase as the area of the subdomain (division by 10^6 to reach good values with the default GUI input)
     309U_tps_sub=zeros(length(X),NbSubDomain);%default spline
     310V_tps_sub=zeros(length(X),NbSubDomain);%default spline
     311U_smooth=zeros(length(X),1);
     312V_smooth=zeros(length(X),1);
     313
     314nb_select=zeros(length(X),1);
    305315FF=zeros(length(X),1);
    306 for isubx=1:NbSubDomainX
    307     for isuby=1:NbSubDomainY
    308         SubRangx(isuby,isubx,:)=[CentreX(isubx)-SizX/2 CentreX(isubx)+SizX/2];
    309         SubRangy(isuby,isubx,:)=[CentreY(isubx)-SizY/2 CentreY(isubx)+SizY/2];
    310         for iter=1:3 %increase the subdomain during three iterations at most
    311             ind_sel=find(X>SubRangx(isuby,isubx,1) & X<SubRangx(isuby,isubx,2) & Y>SubRangy(isuby,isubx,1) & Y<SubRangy(isuby,isubx,2)); 
    312             size(ind_sel)
    313             if numel(ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
    314                 SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
    315                 SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4;
    316                 SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4;
    317                 SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
     316test_empty=zeros(1,NbSubDomain);
     317for isub=1:NbSubDomain
     318    SubRangx(isub,:)=[CentreX(isub)-SizX/2 CentreX(isub)+SizX/2];
     319    SubRangy(isub,:)=[CentreY(isub)-SizY/2 CentreY(isub)+SizY/2];
     320    ind_sel_previous=[];
     321    ind_sel=0;
     322    while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most
     323        ind_sel_previous=ind_sel;
     324        ind_sel=find(X>SubRangx(isub,1) & X<SubRangx(isub,2) & Y>SubRangy(isub,1) & Y<SubRangy(isub,2));
     325        % if no vector in the subdomain, skip the subdomain
     326        if isempty(ind_sel)
     327            test_empty(isub)=1;   
     328            U_tps(1,isub)=0;%define U_tps and V_tps by default
     329            V_tps(1,isub)=0;
     330            break
     331            % if too few selected vectors, increase the subrange for next iteration
     332        elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
     333            SubRangx(isub,1)=SubRangx(isub,1)-SizX/4;
     334            SubRangx(isub,2)=SubRangx(isub,2)+SizX/4;
     335            SubRangy(isub,1)=SubRangy(isub,1)-SizY/4;
     336            SubRangy(isub,2)=SubRangy(isub,2)+SizY/4;
     337        else
     338            [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),U(ind_sel),rho);
     339            [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),V(ind_sel),rho);
     340            UDiff=U_smooth_sub-U(ind_sel);
     341            VDiff=V_smooth_sub-V(ind_sel);
     342            NormDiff=UDiff.*UDiff+VDiff.*VDiff;
     343            FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
     344            ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors
     345            % no value exceeds threshold, the result is recorded
     346            if isequal(numel(ind_ind_sel),numel(ind_sel))
     347                U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub;
     348                V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub;
     349                nbpoints(isub)=numel(ind_sel);
     350                X_tps(1:nbpoints(isub),isub)=X(ind_sel);
     351                Y_tps(1:nbpoints(isub),isub)=Y(ind_sel);
     352                U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub;
     353                V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub;         
     354                nb_select(ind_sel)=nb_select(ind_sel)+1;
     355                 display('good')
     356                break
     357                % too few selected vectors, increase the subrange for next iteration
     358            elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous);
     359                SubRangx(isub,1)=SubRangx(isub,1)-SizX/4;
     360                SubRangx(isub,2)=SubRangx(isub,2)+SizX/4;
     361                SubRangy(isub,1)=SubRangy(isub,1)-SizY/4;
     362                SubRangy(isub,2)=SubRangy(isub,2)+SizY/4;
     363%                 display('fewsmooth')
     364                % interpolation-smoothing is done again with the selected vectors
    318365            else
    319                 [U_smooth_sub,U_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),U(ind_sel),rho);
    320                 [V_smooth_sub,V_tps_sub]=tps_uvmat(X(ind_sel),Y(ind_sel),V(ind_sel),rho);
    321                 size(U_smooth_sub)
    322                 size(U(ind_sel))
    323                 UDiff=U_smooth_sub-U(ind_sel);
    324                 VDiff=V_smooth_sub-V(ind_sel);
    325                 NormDiff=UDiff.*UDiff+VDiff.*VDiff;
    326                 FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination
    327                 ind_ind_sel=find(FF(ind_sel)==0);
    328                 if isequal(numel(ind_ind_sel),numel(ind_sel))
    329                     U_smooth(ind_sel,isuby,isubx)=U_smooth_sub;
    330                     V_smooth(ind_sel,isuby,isubx)=V_smooth_sub;
    331                     break
    332                 elseif numel(ind_ind_sel)<SubDomain/4;% too few selected vectors, increase the subrange for next iteration
    333                     SubRangx(isuby,isubx,1)=SubRangx(isuby,isubx,1)-SizX/4;
    334                     SubRangx(isuby,isubx,2)=SubRangx(isuby,isubx,2)+SizX/4;
    335                     SubRangy(isuby,isubx,1)=SubRangy(isuby,isubx,1)-SizY/4;
    336                     SubRangy(isuby,isubx,2)=SubRangy(isuby,isubx,2)+SizY/4;
    337                 else
    338                     [U_smooth(ind_sel(ind_ind_sel),isuby,isubx),U_tps(ind_sel(ind_ind_sel),isuby,isubx),U_tps3(isuby,isubx)]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),U(ind_sel(ind_ind_sel)),rho);
    339                     [V_smooth(ind_sel(ind_ind_sel),isuby,isubx),V_tps(ind_sel(ind_ind_sel),isuby,isubx),V_tps3(isuby,isubx)]=tps_uvmat(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),V(ind_sel(ind_ind_sel)),rho);
    340                      break                 
    341                 end
    342             end       
    343         end       
    344         X_ctrs(isuby,isubx)=mean(X(ind_sel));%gravity centre of selected points
    345         Y_ctrs(isuby,isubx)=mean(Y(ind_sel));%positions of tps sources for the subdomain i,     
    346         dist_ctre(ind_sel,isuby,isubx)=sqrt(abs(((X(ind_sel)-X_ctrs(isuby,isubx)).*(X(ind_sel)-X_ctrs(isuby,isubx)))+((Y(ind_sel)-Y_ctrs(isuby,isubx)).*(Y(ind_sel)-Y_ctrs(isuby,isubx)))));
    347     end
    348 end
    349 U_patch=sum(sum(U_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
    350 V_patch=sum(sum(V_smooth.*dist_ctre,3),2)./sum(sum(dist_ctre,3),2);
    351 
    352 %------------------------------------------------------------------------
    353 %fasshauer@iit.edu MATH 590 ? Chapter 19 32
    354 % X,Y initial coordiantes
    355 % XI vector, YI column vector for the grid of interpolation points
    356 function [U_smooth,U_tps,U_tps3]=tps_uvmat(X,Y,U,rho)
    357 %------------------------------------------------------------------------
    358 %rho smoothing parameter
    359 ep = 1;
    360 X=reshape(X,[],1);
    361 Y=reshape(Y,[],1);
    362 rhs = reshape(U,[],1);
    363 % if exist('FF','var')
    364 % test_false=isnan(rhs)|FF~=0;
    365 % else
    366 %     test_false=isnan(rhs);
    367 % end
    368 % X(test_false)=[];
    369 % Y(test_false)=[];
    370 % rhs(test_false)=[];
    371 %randn('state',3); rhs = rhs + 0.03*randn(size(rhs));
    372 rhs = [rhs; zeros(3,1)];
    373 dsites = [X Y];% coordinates of measurement sites
    374 ctrs = dsites;%radial base functions are located at the measurement sites
    375 DM_data = DistanceMatrix(dsites,ctrs);%2D matrix of distances between spline centres (=initial points) ctrs
    376 % if size(XI,1)==1 && size(YI,2)==1 % XI vector, YI column vector
    377 %      [XI,YI]=meshgrid(XI,YI);
    378 % end
    379 % [npy,npx]=size(XI);
    380 % epoints = [reshape(XI,[],1) reshape(YI,[],1)];
    381 IM_sites = tps(ep,DM_data);%values of thin plate at site points
    382 IM = IM_sites + rho*eye(size(IM_sites));%  rho=1/(2*omega) , omega given by fasshauer;
    383 PM=[ones(size(dsites,1),1) dsites];
    384 IM=[IM PM; [PM' zeros(3,3)]];
    385 %fprintf('Condition number estimate: %e\n',condest(IM))
    386 %DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
    387 %EM = tps(ep,DM_eval);%values of thin plate
    388 %PM = [ones(size(epoints,1),1) epoints];
    389 %EM = [EM PM];
    390 U_tps=(IM\rhs);
    391 PM = [ones(size(dsites,1),1) dsites];
    392 EM = [IM_sites PM];
    393 U_smooth=EM *U_tps;
    394 U_tps3=U_tps(end-2:end);
    395 U_tps=U_tps(1:end-3);
     366                [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),U(ind_sel(ind_ind_sel)),rho);
     367                [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),V(ind_sel(ind_ind_sel)),rho);
     368                U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub;
     369                V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub;
     370                nbpoints(isub)=numel(ind_ind_sel);
     371                X_tps(1:nbpoints(isub),isub)=X(ind_sel(ind_ind_sel));
     372                Y_tps(1:nbpoints(isub),isub)=Y(ind_sel(ind_ind_sel));
     373                U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub;
     374                V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub;
     375                nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1;
     376                display('good2')
     377                break
     378            end
     379        end
     380    end
     381end
     382ind_empty=find(test_empty);
     383%remove empty subdomains
     384if ~isempty(ind_empty)
     385    SubRangx(ind_empty,:)=[];
     386    SubRangy(ind_empty,:)=[];
     387    X_tps(:,ind_empty)=[];
     388    Y_tps(:,ind_empty)=[];
     389    U_tps(:,ind_empty)=[];
     390    V_tps(:,ind_empty)=[];
     391end
     392nb_select(nb_select==0)=1;%ones(size(find(nb_select==0)));
     393U_smooth=U_smooth./nb_select;
     394V_smooth=V_smooth./nb_select;
     395
     396
     397
     398
    396399
    397400
    … …  
    421424
    422425
    423 %   DM:     MxN matrix whose i,j position contains the Euclidean
    424 %              distance between the i-th data site and j-th center
    425   function DM = DistanceMatrix(dsites,ctrs)
    426   [M,s] = size(dsites); [N,s] = size(ctrs);
    427   DM = zeros(M,N);
    428   % Accumulate sum of squares of coordinate differences
    429   % The ndgrid command produces two MxN matrices:
    430   %   dr, consisting of N identical columns (each containing
    431   %       the d-th coordinate of the M data sites)
    432   %   cc, consisting of M identical rows (each containing
    433   %       the d-th coordinate of the N centers)
    434   for d=1:s
    435      [dr,cc] = ndgrid(dsites(:,d),ctrs(:,d));
    436      DM = DM + (dr-cc).^2;
    437   end
    438   DM = sqrt(DM);
    439 
    440 
    441   % rbf = tps(e,r)
    442 % Defines thin plate spline RBF
    443 function rbf = tps(e,r)
    444 rbf = zeros(size(r));
    445 nz = find(r~=0);   % to deal with singularity at origin
    446 rbf(nz) = (e*r(nz)).^2.*log(e*r(nz));
    447 
     426
     427
Note: See TracChangeset for help on using the changeset viewer.