function Localcurvature
%Program to compute curvature patterns for characteristic leaf surface (principal curvatures, mean curvature, Gaussian curvature). 
%Next program to run: Display_Characteristicleaf_heights_in2D (optional) or
%go directly to Display_local_curvature_data.
%lines 119- 190: this part was slightly adapted from code originally
%written by Dumais and Kwiatkowska (Plant Journal 2002, 31:229-241)

days=7:19;
%days=input('Enter DAS numbers contained in this data set ');

%load(['S:\Camille Girard-Bock\Shape paper_Summer2013\From Anne-Gaelle\MeanleavesDAS7_19cSep24v4.mat'])
%load(['S:\Camille Girard-Bock\Shape paper_Summer2013\From Anne-Gaelle\Sep 25 - new mean leaf data\MeanleavesDAS7_19cSep24v5testborder.mat']);
%load('S:\Camille Girard-Bock\Shape paper_Summer2013\Lauren - mean outlines\meanoutlines3D_DAS7-19.mat');
Path1=uigetdir(cd,'Select Reference data folder');%AG July 2013: for saving mean leaf data
load([Path1 '\meanoutlines_DAS7-19.mat']);
Pathname=uigetdir(cd,'Select main project folder in  which all shape results and figures are stored');%AG July 2013: for saving mean leaf data
load([Pathname '\3D surface and figures\Mean3Dsurface.mat']); 

resdir='Local curvature data and figures';
if ~exist([Pathname '\' resdir],'dir')
    mkdir(Pathname,resdir);
end
Pathname=[Pathname '\' resdir];
%PathName=uigetdir(pwd,'Select folder in  which to save local curvature data for characteristic leaf shape');


%load('C:\Users\Administrator\Dropbox\Lauren_shapepaper\SEP29_2013_GoodDataforMS\TESTMEANLEAF_SEP27_noextrap_borderclassok_V9\MeanleavesDAS7_19cSep27v9noextra.mat');
%load('C:\Users\Administrator\Dropbox\Lauren_shapepaper\SEP29_2013_GoodDataforMS\meanoutlines3D_DAS7-19.mat');

    Vtemp_all={}; Curv1_all={}; Curv2_all={}; curv1_all={}; curv2_all={}; gcurv_all={}; mcurv_all={}; alpha1_all={}; alpha2_all={}; XIg_all={}; YIg_all={}; ZIg_all={};
for ii=days(1):days(end)
    %DAS=ii+6;
    XI=leaf_XIg{ii};
    YI=leaf_YIg{ii};
    ZI=Mode_leaf_ZIg{ii};
    sizeXI=size(XI);
    mX=sizeXI(1,1);
    nX=sizeXI(1,2);
    
    sizeYI=size(YI);
    mY=sizeYI(1,1);
    nY=sizeYI(1,2);
    
    sizeZI=size(ZI);
    mZ=sizeZI(1,1);
    nZ=sizeZI(1,2);
    
    XI2=reshape(XI,mX*nX,1);
    YI2=reshape(YI,mY*nY,1);
    ZI2=reshape(ZI,mZ*nZ,1);
    indNan=isnan(XI2);
    XI2(indNan)=[];
    YI2(indNan)=[];
    ZI2(indNan)=[];
    XYZ=[];
    XYZ=[XI2,YI2,ZI2];
    %XYZG= true scale data realigned in 3D
    draw_mesh=1;
    if draw_mesh==1
        step=(abs(floor(min(avx{ii})))+ceil(max(avx{ii})))/40;
        ix=floor(min(avx{ii})):step:ceil(max(avx{ii}));
        iy=floor(min(avy{ii})):step:ceil(max(avy{ii}));
        [XIg,YIg] = meshgrid(ix,iy);
        
        %AG - added Sep 24 2013
        pts=[XYZ(:,1),XYZ(:,2),XYZ(:,3)];
        pts=unique(pts,'rows');
        indx=find(~isnan(pts(:,3)));
        pts=pts(indx,:);
        %F=scatteredInterpolant(pts(:,1),pts(:,2),pts(:,3),'linear','linear');
        F=scatteredInterpolant(pts(:,1),pts(:,2),pts(:,3),'linear','none');
        %F=TriScatteredInterp(pts(:,1),pts(:,2),pts(:,3));
        ZIg=F(XIg,YIg);
        
        [IN ON] = inpolygon(XIg,YIg,avx{ii},avy{ii});
        INON=IN+ON;
        ZIg(INON==0)=NaN;
        
        
%          F=TriScatteredInterp(XYZ(:,1),XYZ(:,2),XYZ(:,3));
%          ZIg=F(XIg,YIg);
    end
    axis equal
    %%%%%CALC LOCAL CURVATURES AND THEIR ANGLE VALUES%%%%%
    n=5;
    Vtemp={};
    mapc=jet(10);
    a=1;
    for i=1:size(XIg,1) % go through grid rows
        for j=1:size(XIg,2); % go through grid columns
            a=a+1;
            if a==11
                a=1;
            end;
            %only keep points in polygon
%             [IN ON] = inpolygon(XIg(i,j),YIg(i,j),avx{ii},avy{ii});
%             if IN+ON==0
%                 ZIg(i,j)=NaN;
%             end;
            if ~isnan(ZIg(i,j))
                %Vtemp{i,j}=[XIg(i,j),YIg(i,j),ZIg(i,j)];
                Vtemp{i,j}=[];
                %2 left
                indleft=[j-n:j]; indleft=indleft(indleft>0);
                % 2right
                indright=[j+1:j+n]; indright=indright(indright<size(XIg,2));
                indlr=cat(2,indleft,indright);
                % 2top
                indbot=[i-n:i]; indbot=indbot(indbot>0);
                %2 down
                indtop=[i+1:i+n]; indtop=indtop(indtop<size(XIg,1));
                indbt=cat(2,indbot,indtop);
                for k=indlr
                    for l=indbt
                        Vtemp{i,j}=cat(1,Vtemp{i,j},[XIg(l,k),YIg(l,k),ZIg(l,k)]);
                    end
                end
                V=Vtemp{i,j};
                %AG April 2013: remove neighbours with NaN values
                
                %% This part was slightly adapted from code originally written by Dumais and Kwiatkowska (Plant Journal 2002, 31:229-241 )
                indd=find(isnan(V(:,3)));
                V(indd,:)=[];
                Vt=ones(size(V));
                s3=size(V,1);
                Vxy=[ones(s3,1), V(:,1), V(:,2)];
                vz=V(:,3);
                r=inv(Vxy'*Vxy)*Vxy'*vz;
                kt=-[r(2),r(3),-1]/norm([r(2),r(3),-1]);
                it=[1, 1, -(kt(1)+kt(2))/kt(3)]; it=it/norm(it);
                jt=cross(kt,it);
                normal=kt;
                CT=[it;jt;kt];
                c=mean(V,1);
                c=(CT*c')';
                for l=1:size(V,1),
                    Vt(l,:)=(CT*V(l,:)')';
                    Vt(l,:)=Vt(l,:)-c;
                end
                X=ones(size(V,1), 6);
                for l=1:size(V,1)
                    X(l,1)=Vt(l,1)^2;
                    X(l,2)=Vt(l,1)*Vt(l,2);
                    X(l,3)=Vt(l,2)^2;
                    X(l,4)=Vt(l,1);
                    X(l,5)=Vt(l,2);
                end
                b=inv(X'*X)*X'*Vt(1:size(V,1),3);
                xu=[1 0 b(4)];
                xv=[0 1 b(5)];
                xuu=[0 0 2*b(1)];
                xuv=[0 0 b(2)];
                xvv=[0 0 2*b(3)];
                E=dot(xu,xu);
                F=dot(xu,xv);
                G=dot(xv,xv);
                e=det([xuu; xu; xv])/sqrt(E*G-F^2);
                f=det([xuv; xu; xv])/sqrt(E*G-F^2);
                g=det([xvv; xu; xv])/sqrt(E*G-F^2);
                A=F*g-G*f;
                B=E*g-G*e;
                C=E*f-F*e;
                %get roots of A x^2 - B x +C
                delta=(B^2)- (4*A*C);
                X1=(-B-sqrt(delta))/(2*A);
                X2=(-B+sqrt(delta))/(2*A);
                k1=(e+f*X1)/(E+F*X1);
                k2=(e+f*X2)/(E+F*X2);
                K1=[1, X1, 0];
                K2=[1, X2, 0];
                
                if abs(k1)<abs(k2)
                    k1new=k2;
                    k2new=k1;
                    K1new=K2;
                    K2new=K1;
                    K1=K1new;
                    K2=K2new;
                    k1=k1new;
                    k2=k2new;
                end

                K1=CT'*(K1)'; K1=k1*K1'/norm(K1);
                K2=CT'*(K2)'; K2=k2*K2'/norm(K2);
                Curv1{i,j}=K1;
                Curv2{i,j}=K2;
                curv1(i,j)=k1;
                curv2(i,j)=k2;
                gcurv(i,j)=k1*k2;
                mcurv(i,j)=(k1+k2)/2;
                
                %% End of part adapted from Dumais and Kwiatkowska 2002.
                %angles for K1 & K2:
                alpha1(i,j)=atan2(K1(2),K1(1))*180/pi;
                alpha2(i,j)=atan2(K2(2),K2(1))*180/pi;
                if alpha1(i,j)<-90
                    alpha1(i,j)=alpha1(i,j)+180;
                end
                if alpha1(i,j)>90
                    alpha1(i,j)=alpha1(i,j)-180;
                end
                if alpha2(i,j)<-90
                    alpha2(i,j)=alpha2(i,j)+180;
                end
                if alpha2(i,j)>90
                    alpha2(i,j)=alpha2(i,j)-180;
                end
                
            else Vtemp{i,j}=[NaN NaN NaN];
                Curv1{i,j}=[NaN NaN NaN];
                Curv2{i,j}=NaN;
                curv1(i,j)=NaN;
                curv2(i,j)=NaN;
                gcurv(i,j)=NaN;
                mcurv(i,j)=NaN;
                alpha1(i,j)=NaN;
                alpha2(i,j)=NaN;
            end
        end
    end

    %this above calc cuvartures
    %record the data for all days:
    Vtemp_all{ii}=Vtemp; %Curv1_all{ii}=Curv1; Curv2_all{ii}=Curv2; 
    curv1_all{ii}=curv1; curv2_all{ii}=curv2; gcurv_all{ii}=gcurv; 
    mcurv_all{ii}=mcurv; alpha1_all{ii}=alpha1; alpha2_all{ii}=alpha2;
    XIg_all{ii}=XIg; YIg_all{ii}=YIg; ZIg_all{ii}=ZIg;
    clear Vtemp Curv1 Curv2 curv1 curv2 gcurv mcurv alpha1 alpha2 XIg YIg ZIg;
end

cd(Pathname)
save([Pathname '\Characteristicleaves_localcurvaturedata.mat'], 'Vtemp_all','curv1_all','curv2_all','gcurv_all','mcurv_all','alpha1_all','alpha2_all','XIg_all','YIg_all','ZIg_all');

disp(['Curvature data have been saved to ' Pathname '\Characteristicleaves_localcurvaturedata.mat']);
% % % we register the curvature data then can display the information needed with another program




