clearvars

sdirectory2 = 'C:\Users\Ellis.Newham\Desktop\academic work\cementum book\m104135\'; %insert folder containing straightened images here

tifFiles2 = dir([sdirectory2 '*.tif']);

NumberOfImages=length(tifFiles2);

for j=1:NumberOfImages
    clearvars imageTortuosity midLine
    filename=[sdirectory2 tifFiles2(j).name];
    img=imread(filename);
    %[m,n]=size(img);
    cem=img;
    [m,n]=size(cem);
    for i=1:n
        colVec=cem(:,i);
        [row,col]=find(colVec);
        clipOff(i)=min(row);
    end
    cem=cem(((max(clipOff))+5):end,:);
    [m,n]=size(cem);
    cem=cem(1:m,:);
    imshow(cem)
    %cem=cem(1:107,:);
    cem=double(cem);
    %Elis Newham
    %load image, call it cem
    %then normalize
    cemMean=mean(cem(:));
    fromAvr=cem-cemMean;
    stdev=std(cem(:));
    norm=fromAvr-stdev;
    cem=norm;
    cem1=cem;
    %call file "cem1", and change to double if uint8 format
    %image is now centred around mean - no "fromavr" formatting needed

    %height and depth matrices
    height=cem1;
    height(height<0.0001)=0;
    depth=cem1;
    depth(depth>-0.0001)=0;

    %Average height above mean
    heightVector=height(:);
    heightVector(heightVector<0.00001)=[];
    mean(heightVector);
    meanHeight=ans;
    meanHeightCem1=meanHeight;

    %root mean square height
    heightSqr=height.*height;
    heightSqrVector=heightSqr(:);
    meanSqrHeight=mean(heightSqrVector);
    rootMeanSqrHeight=sqrt(meanSqrHeight);
    rootMeanSqrHeightCem1=rootMeanSqrHeight;

    %maximum peak height
    peaks=extrema2(height);
    max(peaks);
    maxPeakHeight=ans;
    maxPeakHeightCem1=maxPeakHeight;

    %average maximum peak height
    sortedPeaks=sort(peaks,'descend');
    top10Peaks=sortedPeaks(1:10);
    mean(top10Peaks);
    avrMaximumPeakHeight=ans;
    avrMaximumPeakHeightcem1=avrMaximumPeakHeight;

    %volume of peaks
    peakVol=sum(peaks);

    %minimum valley depth
    depthSqr=depth.*depth;
    inverseDepth=sqrt(depthSqr);
    valleys=extrema2(inverseDepth);
    max(valleys);
    minValleyDepth=ans-(ans.*2);
    minValleyDepthCem1=minValleyDepth;


    %average minimum valley depth
    sortedValleys=sort(valleys,'descend');
    top10Valleys=sortedValleys(1:10);
    mean(top10Valleys);
    avrMinimumValleyDepth=ans-(ans.*2);
    avrMinimumValleyDepthCem1=avrMinimumValleyDepth;

    %average maximum height
    avrMaximumHeightCem1=avrMaximumPeakHeightcem1-avrMinimumValleyDepthCem1;

    %10 point average maximum height
    TenPointAvrMaximumHeightCem1=avrMaximumPeakHeightcem1+(-avrMinimumValleyDepthCem1);


    %maximum skew
    skew=skewness(cem1);
    max(skew);
    maximumSkew=ans;
    maximumSkewCem1=maximumSkew;


    %summit density
    [m,n]=size(peaks);
    [m2,n2]=size(cem1);
    cem1Area=m2.*n2;
    summitDensity=m/cem1Area;

    %maximum kurtosis
    kurtosisCem1=kurtosis(cem1);
    max(kurtosisCem1);
    maxkurtosis=ans;
    maxKurtosisCem1=maxkurtosis;


    %RMS gradient
    [Gmag,Gdir]=imgradient(cem1);
    GmagSqr=Gmag.*Gmag;
    GmagSqrVector=GmagSqr(:);
    mean(GmagSqrVector);
    RMSgrad=sqrt(ans);
    RMSgradCem1=RMSgrad;


    %non-dimensional (1 pixel equals 1) surface area
    dx=1;
    dy=1;
    Z=cem1;
    [m,n]=size(Z);
    areas = 0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ... 
    (dy*(Z(2:m,1:n-1) - Z(1:m-1,1:n-1))).^2) + ... 
    0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ... 
    (dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2);
    zMean = 0.25 * (Z(1:m-1,1:n-1) + Z(1:m-1,2:n) + Z(2:m,1:n-1) + Z(2:m,2:n));
    surfaceArea = sum(areas(:));
    surfaceAreaCem1=surfaceArea;

    %Developed interfacial area
    [m,n]=size(cem1);
    totalArea=m.*n;
    Sdr=surfaceArea/totalArea;


    %core depth calculations
    cemMin=min(cem(:));
    cemMin=cemMin+1;
    cem(cem<cemMin)=cemMin;
    h=cdfplot(cem(:));
    y=get(h,'YData');
    x=get(h,'XData');
    y=y(:);
    y2=y(1:2:end);

    y2(y2>0.9)=0;
    y2(y2<0.1)=0;
    [m,n]=size(y);
    figure(2),plot(y2)

    yDiff=diff(y2);
    [Bsort,Bidx]=getNLargeElements(yDiff,30);

    sortedIndex=sort(Bidx(:),'ascend');

    Index=sortedIndex;

    indexEdited=diff(Index);

    indexEdited=[1;indexEdited];

    indexEdited(indexEdited>3)=0;

    indexEdited(indexEdited>1)=1;

    indexIndice=indexEdited.*Index;

    indexIndice(indexIndice<1)=[];

    Indexindice=sort(indexIndice,'ascend');

    diffIndice=diff(Indexindice);

    diffIndice=[1;diffIndice];

    diffIndice(diffIndice>3)=0;

    diffIndice(diffIndice>1)=1;

    refinedIndice=Indexindice.*diffIndice;

    refinedIndice(refinedIndice<1)=[];

    refinedIndice=sort(refinedIndice,'descend');

    diffIndice2=-diff(refinedIndice);

    diffIndice2=[1;diffIndice2];

    diffIndice2(diffIndice2>3)=0;

    diffIndice2(diffIndice2>1)=1;

    refinedIndice2=diffIndice2.*refinedIndice;

    refinedIndice2(refinedIndice2<1)=[];

    Refinedindice2=sort(refinedIndice2,'descend');

    diffIndice3=-diff(Refinedindice2);

    diffIndice3=[1;diffIndice3];

    diffIndice3(diffIndice3>3)=0;

    diffIndice3(diffIndice3>1)=1;

    refinedIndice3=Refinedindice2.*diffIndice3;

    refinedIndice3(refinedIndice3<1)=[];

    diffIndice4=-diff(refinedIndice3);

    diffIndice4=[1;diffIndice4];

    diffIndice4(diffIndice4>3)=0;

    diffIndice4(diffIndice4>1)=1;

    refinedIndice4=refinedIndice3.*diffIndice4;

    refinedIndice4(refinedIndice4<1)=[];

    refinedIndice4(refinedIndice4<1)=[];

    diffIndice5=-diff(refinedIndice4);

    diffIndice5=[1;diffIndice5];

    diffIndice5(diffIndice5>3)=0;

    diffIndice5(diffIndice5>1)=1;

    refinedIndice5=refinedIndice4.*diffIndice5;

    refinedIndice5(refinedIndice5<1)=[];

    diffIndice6=-diff(refinedIndice5);

    diffIndice6=[1;diffIndice6];

    diffIndice6(diffIndice6>3)=0;

    diffIndice6(diffIndice6>1)=1;

    refinedIndice6=refinedIndice5.*diffIndice6;

    refinedIndice6(refinedIndice6<1)=[];

    diffIndice7=-diff(refinedIndice6);

    diffIndice7=[1;diffIndice7];

    diffIndice7(diffIndice7>3)=0;

    diffIndice7(diffIndice7>1)=1;

    refinedIndice7=refinedIndice6.*diffIndice7;

    refinedIndice7(refinedIndice7<1)=[];

    diffIndice8=-diff(refinedIndice7);

    diffIndice8=[1;diffIndice8];

    diffIndice8(diffIndice8>3)=0;

    diffIndice8(diffIndice8>1)=1;

    refinedIndice8=refinedIndice7.*diffIndice8;

    refinedIndice8(refinedIndice8<1)=[];

    diffIndice9=-diff(refinedIndice8);

    diffIndice9=[1;diffIndice9];

    diffIndice9(diffIndice9>3)=0;

    diffIndice9(diffIndice9>1)=1;

    refinedIndice9=refinedIndice8.*diffIndice9;

    refinedIndice9(refinedIndice9<1)=[];

    coreHeight=max(refinedIndice8);
    coreHeight=coreHeight+cemMin;

    coreDepth=min(refinedIndice8);
    coreDepth=coreDepth+cemMin;

    Core=coreHeight-coreDepth;

    T=table(coreHeight,coreDepth,Core)

    % Texture parameters following core depth calculations. core height and
    % depth must be known and inputted in as "coreHeight" and "coreDepth"
    % respectively. coreDepth must also be positive.

    %mean height of peaks above core height
    mean(cem1(:));
    cem1Mean=ans;
    cem1Mean=cem1Mean;
    Height=coreHeight-cem1Mean;
    Depth=cem1Mean-coreDepth;
    fromAvr=cem1-cem1Mean;
    aboveCore=fromAvr;
    aboveCore(aboveCore<0.000001)=0;
    peaksAboveCore=extrema2(aboveCore);
    peaksAboveCore(peaksAboveCore<(Height))=[];
    mean(peaksAboveCore);
    meanPeakHeightAboveCore=ans;
    meanPeakHeightAboveCoreCem1=meanPeakHeightAboveCore;


    %mean depth of valleys below core
    belowCore=fromAvr;
    belowCore(belowCore>-0.0000001)=0;
    belowCoreSqr=belowCore.*belowCore;
    rootbelowCore=sqrt(belowCoreSqr);
    valleysBelowCore=extrema2(rootbelowCore);
    n=Depth;
    valleysBelowCore(valleysBelowCore<(n))=[];
    mean(valleysBelowCore);
    meanValleysBelowCore=ans;
    meanValleysBelowCoreCem1=meanValleysBelowCore;


    % surface-bearing area ratio
    aboveCore=cem1-(cem1Mean+coreHeight);
    aboveCore(aboveCore<0.000001)=0;
    dx=1;
    dy=1;
    Z=aboveCore;
    [m,n]=size(Z);
    areas = 0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ... 
    (dy*(Z(2:m,1:n-1) - Z(1:m-1,1:n-1))).^2) + ... 
    0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ... 
    (dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2);
    zMean = 0.25 * (Z(1:m-1,1:n-1) + Z(1:m-1,2:n) + Z(2:m,1:n-1) + Z(2:m,2:n));
    surfaceArea = sum(areas(:));
    surfaceAreaAboveCoreheight=surfaceArea;
    [m,n]=size(cem1);
    totalArea=m.*n;
    SurfaceBearingAreaRatioCem1=surfaceAreaAboveCoreheight/totalArea;

    results2=table(meanPeakHeightAboveCoreCem1,meanValleysBelowCoreCem1,SurfaceBearingAreaRatioCem1); 


    %Autocorrelation length
    cemMean=mean(cem(:));
    fromAvr=cem-cemMean;
    Height=fromAvr;
    Height(Height<0.0001)=0;
    Depth=fromAvr;
    Depth(Depth>-0.0001)=0;
    heightSqr=Height.*Height;
    heightSqrVector=heightSqr(:);
    meanSqrHeight=mean(heightSqrVector);
    rootMeanSqrHeight=sqrt(meanSqrHeight);
    rootMeanSqrHeightCem1=rootMeanSqrHeight;
    b=rootMeanSqrHeightCem1;

    [m,n]=size(cem);

    for i=1:m
        a=bsxfun(@times,cem,cem(i,1:n));
        fromNorm=a-b;
        fromNorm(fromNorm<0)=0;
        avrFromNorm=mean(fromNorm);
        final=mean(avrFromNorm);
        horizVec(i,:)=final;

    end

    for i=1:n
        a2=bsxfun(@times,fromAvr,fromAvr(1:m,i));
        fromNorm2=a2-b;
        fromNorm2(fromNorm2<0)=0;
        avrFromNorm2=mean(fromNorm2);
        final2=mean(avrFromNorm2);
        vertVec(i,:)=final2;
    end

    figure(3),plot(ac2rc(horizVec))
    figure(4),plot(horizVec)
    AC=ac2rc(horizVec);
    [negVal negLoc]=min(AC);
    [posVal posLoc]=max(AC);
    negativeLength=(negLoc/m).*100;
    positiveLength=(posLoc/m).*100;
    ACTable=table(negativeLength, positiveLength);

    fullResults=table(cem1Mean,meanHeightCem1,rootMeanSqrHeightCem1,maxPeakHeightCem1,avrMaximumPeakHeightcem1,peakVol,minValleyDepthCem1,avrMinimumValleyDepthCem1,avrMaximumHeightCem1,TenPointAvrMaximumHeightCem1,maximumSkewCem1,summitDensity,maxKurtosisCem1,RMSgradCem1,surfaceAreaCem1,Sdr,coreHeight,coreDepth,Core,meanPeakHeightAboveCoreCem1,meanValleysBelowCoreCem1,SurfaceBearingAreaRatioCem1,negativeLength,positiveLength);
    
    fname = [sdirectory2 tifFiles2(j).name ' texture results.xlsx'];
    
    writetable(fullResults,fname);
end