clearvars

sdirectory2 = '/Users/elisnewham/Desktop/Elis cementum work/8 bit straightened/k39/filter 2/'; %insert folder location for filtered images here
sdirectory3 = '/Users/elisnewham/Desktop/Elis cementum work/8 bit straightened/k39/filter 2/binaries/'; %insert folder location for binaries here
sdirectory4 = '/Users/elisnewham/Desktop/Elis cementum work/8 bit straightened/k39/'; %insert original foloder location for straightened images here
sdirectory5 = '/Users/elisnewham/Desktop/Elis cementum work/8 bit straightened/k39/increment models/'; %insert output folder location for increment models here

imgName = dir([sdirectory4]);
tifFiles2 = dir([sdirectory2 '*.tif']);
tifFiles3 = dir([sdirectory3 '*.tif']);
tifFiles4 = dir([sdirectory4 '*.tif']);
txtfiles5 = dir([sdirectory5 '*.tif']);
NumberFilterImages=length(tifFiles2);
NumberMaskImages=length(tifFiles3);
NumberOriginalImages=length(tifFiles4);
%tens=floor(NumberImages/10);
%Tens=tens*10;

for j=1:NumberOriginalImages

filterFilename=[sdirectory2 tifFiles2(j).name];
maskFileName=[sdirectory3 tifFiles3(j).name];
originalFileName=[sdirectory4 tifFiles4(j).name];


original=imread(originalFileName);
original=fliplr(original);
img=load(filterFilename);
mask=load(maskFileName);
mask(mask>0)=1;
imgDouble=im2double(img);
cem=immultiply(imgDouble,mask);
cem=fliplr(cem);
[m,n]=size(cem);
cemProject=cem;

firstIncrementImg=zeros(m,n);
secondIncrementImg=zeros(m,n);
thirdIncrementImg=zeros(m,n);
fourthIncrementImg=zeros(m,n);
fifthIncrementImg=zeros(m,n);
sixthIncrementImg=zeros(m,n);
seventhIncrementImg=zeros(m,n);
eitthIncrementImg=zeros(m,n);
ninthIncrementImg=zeros(m,n);
tenthIncrementImg=zeros(m,n);
%cemProject(cemProject<0.8)=0;
    for i=21:n
        colVec=cem(:,i);
        [row,col]=find(colVec);
        firstPixel=min(row);
        lastPixel=firstPixel+22;
        cemProject(firstPixel:lastPixel,i)=0;
        firstIncrementImg(firstPixel:lastPixel,i)=1;
    end

    firstIncrementBW=im2bw(firstIncrementImg);
    firstIncrementRealVals=immultiply(firstIncrementBW,cemProject);

    firstIncrementSaveProject=cemProject;

    colVec=cemProject(:,21);
    [row,col]=find(colVec);
    firstPixel=min(row);
    colVecDub=im2double(colVec);
    [val,loc,width]=findpeaks(colVecDub);
    firstPeak=loc(1);
    firstPeakVec(21)=firstPeak;
    if firstPeak<(firstPixel+3)
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
        firstPeakWidthVec(21)=firstPeakWidth;
        if firstPeakWidth<5
            lastPixel=firstPeak+firstPeakWidth;
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(21)=totalWidth;
            cemProject(firstPixel:lastPixel,21)=0;
            firstIncrementImg(firstPixel:lastPixel,21)=1;
        end
    end

    for i=22:n
        colVec=cemProject(:,i);
        [row,col]=find(colVec);
        firstPixel=min(row);
        colVecDub=im2double(colVec);
        [val,loc,width]=findpeaks(colVecDub);
sumVal=sum(val);
if sumVal>0
        firstPeak=loc(1);
        firstPeakVec(i)=firstPeak;
        if firstPeak<(firstPixel+3)
            firstPeakWidth=width(1);
            firstPeakWidth=round(firstPeakWidth);
            firstPeakWidthVec(i)=firstPeakWidth;
            if firstPeakWidth<5
                lastPixel=firstPeak+firstPeakWidth;
                totalWidth=lastPixel-firstPixel;
                totalWidthVec(i)=totalWidth;
                cemProject(firstPixel:lastPixel,i)=0;
                firstIncrementImg(firstPixel:lastPixel,i)=1;
            elseif firstPeakWidth>4
                lastPixel=firstPeak+4;
                totalWidth=lastPixel-firstPixel;
                totalWidthVec(i)=totalWidth;
                cemProject(firstPixel:lastPixel,i)=0;
                firstIncrementImg(firstPixel:lastPixel,i)=1;
            end
        end
    end
end

    firstIncrementBW=im2bw(firstIncrementImg);
    firstIncrementRealVals=immultiply(firstIncrementBW,cemProject);

    for i=21:n
        colVec=firstIncrementRealVals(:,i);
        colVecDub=im2double(colVec);
        [row,col]=find(colVecDub);
        [val,loc]=findpeaks(colVecDub);
        [m,n]=size(loc);
        peakCount=m;
        if peakCount>1
            secondPeak=loc(2);
            secondEnd=max(row);
            cemProject(secondPeak:secondEnd,i)=255;
            firstIncrementRealVals((secondPeak-1):secondEnd,i)=0;
        end
    end

    firstIncrementSaveProject=cemProject;

    %firstIncrementBW=im2bw(firstIncrementImg);
    %firstIncrementRealVals=immultiply(firstIncrementBW,cem);        

    %firstincrement now seperated
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    colVec=cemProject(:,21);
    [row,col]=find(colVec);
    firstPixel=min(row);
    firstPixelVec(21)=firstPixel;
    colVecDub=im2double(colVec);
    [val,loc,width]=findpeaks(colVecDub);
    firstPixelVec(24)=firstPixel;
    firstPeak=loc(1);
    firstPeakVec(21);
    firstPeakWidth=width(1);
    firstPeakWidth=round(firstPeakWidth);
    firstPeakWidthVec(21)=firstPeakWidth;
    lastPixel=firstPeak+firstPeakWidth;
    totalWidth=lastPixel-firstPixel;
    totalWidthVec(21)=totalWidth;
    cemProject(firstPixel:lastPixel,21)=0;
    secondIncrementImg(firstPixel:lastPixel,21)=1;

    [m,n]=size(cem);

    for i=22:n
colVec=cemProject(:,i);
sumColVec=sum(colVec);
if sumColVec>0
        [row,col]=find(colVec);
        firstPixel=min(row);
        firstPixelVec(i)=firstPixel;
        colVecDub=im2double(colVec);
        [val,loc,width]=findpeaks(colVecDub);
        firstPeak=loc(1);
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
widthDiff=firstPeak+firstPeakWidth;
        if widthDiff<m
        lastColVec=firstIncrementSaveProject(firstPeak:(widthDiff-1),i-1);
        sumLastColVec=sum(lastColVec);
        if sumLastColVec>1
        if (sum(lastColVec))<1
            cemProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            firstIncrementSaveProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            firstIncrementRealVals(firstPixel:(firstPeak+firstPeakWidth),i)=255;
            colVecDub(firstPixel:(firstPeak+firstPeakWidth))=0;
            %[val,loc,width]=findpeaks(colVecDub);
        end
sumColVec=sum(colVecDub);
if sumColVec>0
        [row,col]=find(colVecDub);
        firstPixel=min(row);
        [val,loc,width]=findpeaks(colVecDub);
sumVal=sum(val);
if sumVal>0
        firstPeak=loc(1);
        firstPeakVec(i)=firstPeak;
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
        firstPeakWidthVec(i)=firstPeakWidth;
        if firstPeakWidth<(firstPeakWidthVec(21))
            lastPixel=firstPeak+firstPeakWidth;
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            secondIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth>(firstPeakWidthVec(21))
            lastPixel=firstPeak+(firstPeakWidthVec(21));
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            secondIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth<((firstPeakWidthVec(21))+1)
            if firstPeakWidth>((firstPeakWidthVec(21))-1)
              lastPixel=firstPeak+(firstPeakWidthVec(21));
                totalWidth=lastPixel-firstPixel;
                totalWidthVec(i)=totalWidth;
                cemProject(firstPixel:lastPixel,i)=0;
                secondIncrementImg(firstPixel:lastPixel,i)=1;
            end
        end
    end
end
end
end
end
end

    secondIncrementBW=im2bw(secondIncrementImg);

[m,n]=size(secondIncrementBW);
[g,h]=size(cemProject);
difference=g-m;
if difference>0
addOn=zeros(difference,n);
secondIncrementBW=[addOn; secondIncrementBW];
end


    secondIncrementRealVals=immultiply(secondIncrementBW,cemProject); 

    secondIncrementSaveProject=cemProject;


    %second increment now seperated
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    colVec=cemProject(:,21);
    [row,col]=find(colVec);
    firstPixel=min(row);
    firstPixelVec(21)=firstPixel;
    colVecDub=im2double(colVec);
    [val,loc,width]=findpeaks(colVecDub);
    firstPixelVec(21)=firstPixel;
    firstPeak=loc(1);
    firstPeakVec(21);
    firstPeakWidth=width(1);
    firstPeakWidth=round(firstPeakWidth);
    firstPeakWidthVec(21)=firstPeakWidth;
    lastPixel=firstPeak+firstPeakWidth;
    totalWidth=lastPixel-firstPixel;
    totalWidthVec(21)=totalWidth;
    cemProject(firstPixel:lastPixel,21)=0;
    thirdIncrementImg(firstPixel:lastPixel,21)=1;

for i=22:n
colVec=cemProject(:,i);
sumColVec=sum(colVec);
if sumColVec>0
        [row,col]=find(colVec);
        firstPixel=min(row);
        firstPixelVec(i)=firstPixel;
        colVecDub=im2double(colVec);
        [val,loc,width]=findpeaks(colVecDub);
        firstPeak=loc(1);
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
widthDiff=firstPeak+firstPeakWidth;
        if widthDiff<m
        lastColVec=secondIncrementSaveProject(firstPeak:(widthDiff),i-1);
        sumLastColVec=sum(lastColVec);
        if sumLastColVec>1
        if (sum(lastColVec))<1
            cemProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            secondIncrementSaveProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            secondIncrementRealVals(firstPixel:(firstPeak+firstPeakWidth),i)=255;
            colVecDub(firstPixel:(firstPeak+firstPeakWidth))=0;
            %[val,loc,width]=findpeaks(colVecDub);
        end
sumColVec=sum(colVecDub);
if sumColVec>0
        [row,col]=find(colVecDub);
        firstPixel=min(row);
        [val,loc,width]=findpeaks(colVecDub);
sumVal=sum(val);
if sumVal>0
        firstPeak=loc(1);
        firstPeakVec(i)=firstPeak;
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
        firstPeakWidthVec(i)=firstPeakWidth;
        if firstPeakWidth<(firstPeakWidthVec(21))
            lastPixel=firstPeak+firstPeakWidth;
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            thirdIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth>(firstPeakWidthVec(21))
            lastPixel=firstPeak+(firstPeakWidthVec(21));
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            thirdIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth<((firstPeakWidthVec(21))+1)
            if firstPeakWidth>((firstPeakWidthVec(21))-1)
              lastPixel=firstPeak+(firstPeakWidthVec(21));
                totalWidth=lastPixel-firstPixel;
                totalWidthVec(i)=totalWidth;
                cemProject(firstPixel:lastPixel,i)=0;
                thirdIncrementImg(firstPixel:lastPixel,i)=1;
            end
        end
    end
end
end
end
end
end

thirdIncrementBW=im2bw(thirdIncrementImg);
[m,n]=size(thirdIncrementBW);
[g,h]=size(cemProject);
difference=g-m;
if difference>0
addOn=zeros(difference,n);
thirdIncrementBW=[addOn; thirdIncrementBW];
end

    thirdIncrementRealVals=immultiply(thirdIncrementBW,cemProject); 

    thirdIncrementSaveProject=cemProject;

    %third increment sorted
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    colVec=cemProject(:,21);
    [row,col]=find(colVec);
    firstPixel=min(row);
    firstPixelVec(21)=firstPixel;
    colVecDub=im2double(colVec);
    [val,loc,width]=findpeaks(colVecDub);
    firstPixelVec(21)=firstPixel;
    firstPeak=loc(1);
    firstPeakVec(21);
    firstPeakWidth=width(1);
    firstPeakWidth=round(firstPeakWidth);
    firstPeakWidthVec(21)=firstPeakWidth;
    lastPixel=firstPeak+firstPeakWidth;
    totalWidth=lastPixel-firstPixel;
    totalWidthVec(21)=totalWidth;
    cemProject(firstPixel:lastPixel,21)=0;
    fourthIncrementImg(firstPixel:lastPixel,21)=1;

    for i=22:n
        colVec=cemProject(:,i);
sumColVec=sum(colVec);
if sumColVec>0
        [row,col]=find(colVec);
        firstPixel=min(row);
        firstPixelVec(i)=firstPixel;
        colVecDub=im2double(colVec);
        [val,loc,width]=findpeaks(colVecDub);
        firstPeak=loc(1);
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
widthDiff=firstPeak+firstPeakWidth;
        if widthDiff<m
        lastColVec=thirdIncrementSaveProject(firstPeak:(widthDiff),i-1);
        sumLastColVec=sum(lastColVec);
        if sumLastColVec>1
        if (sum(lastColVec))<1
            cemProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            thirdIncrementSaveProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            thirdIncrementRealVals(firstPixel:(firstPeak+firstPeakWidth),i)=255;
            colVecDub(firstPixel:(firstPeak+firstPeakWidth))=0;
            %[val,loc,width]=findpeaks(colVecDub);
        end
sumColVec=sum(colVec);
if sumColVec>0
        [row,col]=find(colVecDub);
        firstPixel=min(row);
        [val,loc,width]=findpeaks(colVecDub);
sumVal=sum(val);
if sumVal>0
        firstPeak=loc(1);
        firstPeakVec(i)=firstPeak;
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
        firstPeakWidthVec(i)=firstPeakWidth;
        if firstPeakWidth<(firstPeakWidthVec(21))
            lastPixel=firstPeak+firstPeakWidth;
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            fourthIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth>(firstPeakWidthVec(21))
            lastPixel=firstPeak+(firstPeakWidthVec(21));
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            fourthIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth<((firstPeakWidthVec(21))+1)
            if firstPeakWidth>((firstPeakWidthVec(21))-1)
              lastPixel=firstPeak+(firstPeakWidthVec(21));
                totalWidth=lastPixel-firstPixel;
                totalWidthVec(i)=totalWidth;
                cemProject(firstPixel:lastPixel,i)=0;
                fourthIncrementImg(firstPixel:lastPixel,i)=1;
            end
        end
    end
end
end
end
end
end

fourthIncrementBW=im2bw(fourthIncrementImg);
[m,n]=size(fourthIncrementBW);
[g,h]=size(cemProject);
difference=g-m;
if difference>0
addOn=zeros(difference,n);
fourthIncrementBW=[addOn; fourthIncrementBW];
end
    fourthIncrementRealVals=immultiply(fourthIncrementBW,cemProject);

    fourthIncrementSaveProject=cemProject;

    for i=21:n
        colVec=fourthIncrementRealVals(:,i);
        [row,col]=find(colVec);
        firstPixel=min(row);
        lastPixel=max(row);
        if lastPixel>firstPixel
            midPixel=firstPixel+((lastPixel-firstPixel)/2);
            midLine(i)=midPixel;
        elseif lastPixel<(firstPixel+1)
            midPixel=firstPixel;
            midLine(i)=midPixel;
        end
    end

    %fourth increment sorted
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [m,n]=size(cem);

    colVec=cemProject(:,21);
    [row,col]=find(colVec);
    firstPixel=min(row);
    firstPixelVec(21)=firstPixel;
    colVecDub=im2double(colVec);
    [val,loc,width]=findpeaks(colVecDub);
    firstPixelVec(21)=firstPixel;
    firstPeak=loc(1);
    firstPeakVec(21);
    firstPeakWidth=width(1);
    firstPeakWidth=round(firstPeakWidth);
    firstPeakWidthVec(21)=firstPeakWidth;
    lastPixel=firstPeak+firstPeakWidth;
    totalWidth=lastPixel-firstPixel;
    totalWidthVec(21)=totalWidth;
    cemProject(firstPixel:lastPixel,21)=0;
    fifthIncrementImg(firstPixel:lastPixel,21)=1;

    for i=22:n
        colVec=cemProject(:,i);
sumColVec=sum(colVec);
if sumColVec>0
        [row,col]=find(colVec);
        firstPixel=min(row);
        firstPixelVec(i)=firstPixel;
        colVecDub=im2double(colVec);
        [val,loc,width]=findpeaks(colVecDub);
        firstPeak=loc(1);
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
widthDiff=firstPeak+firstPeakWidth;
        if widthDiff<m
        lastColVec=fourthIncrementSaveProject(firstPeak:(widthDiff),i-1);
        sumLastColVec=sum(lastColVec);
        if sumLastColVec>1
        if (sum(lastColVec))<1
            cemProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            fourthIncrementSaveProject(firstPixel:(firstPeak+firstPeakWidth),i)=0;
            fourthIncrementRealVals(firstPixel:(firstPeak+firstPeakWidth),i)=255;
            colVecDub(firstPixel:(firstPeak+firstPeakWidth))=0;
            %[val,loc,width]=findpeaks(colVecDub);
        end
sumColVec=sum(colVecDub);
if sumColVec>0
        [row,col]=find(colVecDub);
        firstPixel=min(row);
        [val,loc,width]=findpeaks(colVecDub);
valSum=sum(val);
if sumVal>0
        firstPeak=loc(1);
        firstPeakVec(i)=firstPeak;
        firstPeakWidth=width(1);
        firstPeakWidth=round(firstPeakWidth);
        firstPeakWidthVec(i)=firstPeakWidth;
        if firstPeakWidth<(firstPeakWidthVec(21))
            lastPixel=firstPeak+firstPeakWidth;
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            fifthIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth>(firstPeakWidthVec(21))
            lastPixel=firstPeak+(firstPeakWidthVec(21));
            totalWidth=lastPixel-firstPixel;
            totalWidthVec(i)=totalWidth;
            cemProject(firstPixel:lastPixel,i)=0;
            fifthIncrementImg(firstPixel:lastPixel,i)=1;
        elseif firstPeakWidth<((firstPeakWidthVec(21))+1)
            if firstPeakWidth>((firstPeakWidthVec(21))-1)
              lastPixel=firstPeak+(firstPeakWidthVec(21));
                totalWidth=lastPixel-firstPixel;
                totalWidthVec(i)=totalWidth;
                cemProject(firstPixel:lastPixel,i)=0;
                fifthIncrementImg(firstPixel:lastPixel,i)=1;
            end
        end
    end
end
end
end
end
end

    fifthIncrementBW=im2bw(fifthIncrementImg);
[m,n]=size(fifthIncrementBW);
[g,h]=size(cemProject);
difference=g-m;
if difference>0
addOn=zeros(difference,n);
fifthIncrementBW=[addOn; fifthIncrementBW];
end
    fifthIncrementRealVals=immultiply(fifthIncrementBW,cemProject);

    fifthIncrementSaveProject=cemProject;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clearvars midLine gaussLine

[m,n]=size(cem);
firstIncrementTextureImage=zeros(m,n);
for i=21:n
    colVec=firstIncrementBW(:,i);
    [row,col]=find(colVec);
    firstPixel=min(row);
    lastPixel=max(row);
    if lastPixel>firstPixel
        midPixel=firstPixel+((lastPixel-firstPixel)/2);
        midLine(i)=midPixel;
    elseif lastPixel<(firstPixel+1)
        midPixel=firstPixel;
        midLine(i)=midPixel;
    end
end

[m,n]=size(cem);
[g,h]=size(midLine);
midLineRefined=midLine(21:h);
midLineRefined=midLineRefined(2:end);
[m,n]=size(midLineRefined);
cutOffLength=floor(n/5);
gaussLine=zeros(1,h);
startpoint=(cutOffLength*(1-1))+1;
endPoint=(cutOffLength*1);
cutOfMidline=midLineRefined(startpoint:endPoint);
[m,n]=size(cutOfMidline);
x=[1:n];
x=x';
cutOfMidline=cutOfMidline';
f=fit(x,cutOfMidline,'poly4');
y=feval(f,x);
gaussLine(startpoint:endPoint)=y;
nextStartPixel=gaussLine(endPoint);
for i=2:5
    startpoint=(cutOffLength*(i-1))+1;
    endPoint=(cutOffLength*i);
    cutOfMidline=midLineRefined(startpoint:endPoint);
    [m,n]=size(cutOfMidline);
    x=[1:n];
    x=x';
    cutOfMidline=cutOfMidline';
    f=fit(x,cutOfMidline,'poly4');
    y=feval(f,x);
    gaussStartPoint=y(1);
    if gaussStartPoint>nextStartPixel
        difference=gaussStartPoint-nextStartPixel;
        y=y-difference;
    elseif gaussStartPoint<nextStartPixel
        difference=nextStartPixel-gaussStartPoint;
        y=y+difference;
    end
    gaussLine(startpoint:endPoint)=y;
    nextStartPixel=gaussLine(endPoint);
end
minVal=min(gaussLine);
if minVal<0
PositiveMinVal=(round(sqrt(minVal.*minVal)))+1;
gaussLine=gaussLine+PositiveMinVal;
end
[q,r]=size(gaussLine);
[m,n]=size(cem);
yStartPoint=(n-r)-1;
yEndPoint=r;
gaussLineImg=zeros(m,n);
[row,col]=find(gaussLine);
endPoint=max(col);
for i=1:endPoint
midlinePoint=yStartPoint+i;
midLineVal=round(gaussLine(i));
if midLineVal<1
midLineVal=round(gaussLine(i-1));
if midLineVal<1
    midLineVal=round(gaussLine(i-2));
end
end
if midLineVal>5
    gaussLineImg((midLineVal-5):(midLineVal+5),i)=1;
elseif midLineVal<6
    gaussLineImg((midLineVal-(midLineVal-1)):(midLineVal+(midLineVal-1)),i)=1;
end
end

gaussLineBW=im2bw(gaussLineImg);
firstIncrementGaussBW=gaussLineBW;
%finalFirstincrementTexture=immultiply(gaussLineBW,original);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=21:n
    colVec=secondIncrementBW(:,i);
    [row,col]=find(colVec);
    firstPixel=min(row);
    lastPixel=max(row);
    if lastPixel>firstPixel
        midPixel=firstPixel+((lastPixel-firstPixel)/2);
        midLine(i)=midPixel;
    elseif lastPixel<(firstPixel+1)
        midPixel=firstPixel;
        midLine(i)=midPixel;
    end
end
[m,n]=size(cem);
[g,h]=size(midLine);
midLineRefined=midLine(21:h);
midLineRefined=midLineRefined(2:end);
[m,n]=size(midLineRefined);
cutOffLength=floor(n/5);
gaussLine=zeros(1,h);
startpoint=(cutOffLength*(1-1))+1;
endPoint=(cutOffLength*1);
cutOfMidline=midLineRefined(startpoint:endPoint);
[m,n]=size(cutOfMidline);
x=[1:n];
x=x';
cutOfMidline=cutOfMidline';
f=fit(x,cutOfMidline,'poly4');
y=feval(f,x);
gaussLine(startpoint:endPoint)=y;
nextStartPixel=gaussLine(endPoint);
for i=2:5
    startpoint=(cutOffLength*(i-1))+1;
    endPoint=(cutOffLength*i);
    cutOfMidline=midLineRefined(startpoint:endPoint);
    [m,n]=size(cutOfMidline);
    x=[1:n];
    x=x';
    cutOfMidline=cutOfMidline';
    f=fit(x,cutOfMidline,'poly4');
    y=feval(f,x);
    gaussStartPoint=y(1);
    if gaussStartPoint>nextStartPixel
        difference=gaussStartPoint-nextStartPixel;
        y=y-difference;
    elseif gaussStartPoint<nextStartPixel
        difference=nextStartPixel-gaussStartPoint;
        y=y+difference;
    end
    gaussLine(startpoint:endPoint)=y;
    nextStartPixel=gaussLine(endPoint);
end
minVal=min(gaussLine);
if minVal<0
PositiveMinVal=(round(sqrt(minVal.*minVal)))+1;
gaussLine=gaussLine+PositiveMinVal;
end


[q,r]=size(gaussLine);
[m,n]=size(cem);
yStartPoint=(n-r)-1;
yEndPoint=r;
gaussLineImg=zeros(m,n);
[row,col]=find(gaussLine);
endPoint=max(col);
for i=1:endPoint
midlinePoint=yStartPoint+i;
midLineVal=round(gaussLine(i));
if midLineVal<1
midLineVal=round(gaussLine(i-1));
if midLineVal<1
    midLineVal=round(gaussLine(i-2));
end
end
if midLineVal>5
    gaussLineImg((midLineVal-5):(midLineVal+5),i)=1;
elseif midLineVal<6
    gaussLineImg((midLineVal-(midLineVal-1)):(midLineVal+(midLineVal-1)),i)=1;
end
end

gaussLineBW=im2bw(gaussLineImg);
secondIncrementGaussBW=gaussLineBW;
%finalSecondincrementTexture=immultiply(gaussLineBW,original);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=21:n
    colVec=thirdIncrementBW(:,i);
    [row,col]=find(colVec);
    firstPixel=min(row);
    lastPixel=max(row);
    if lastPixel>firstPixel
        midPixel=firstPixel+((lastPixel-firstPixel)/2);
        midLine(i)=midPixel;
    elseif lastPixel<(firstPixel+1)
        midPixel=firstPixel;
        midLine(i)=midPixel;
    end
end
[m,n]=size(cem);
[g,h]=size(midLine);
midLineRefined=midLine(21:h);
midLineRefined=midLineRefined(2:end);
[m,n]=size(midLineRefined);
cutOffLength=floor(n/5);
gaussLine=zeros(1,h);
startpoint=(cutOffLength*(1-1))+1;
endPoint=(cutOffLength*1);
cutOfMidline=midLineRefined(startpoint:endPoint);
[m,n]=size(cutOfMidline);
x=[1:n];
x=x';
cutOfMidline=cutOfMidline';
f=fit(x,cutOfMidline,'poly4');
y=feval(f,x);
gaussLine(startpoint:endPoint)=y;
nextStartPixel=gaussLine(endPoint);
for i=2:5
    startpoint=(cutOffLength*(i-1))+1;
    endPoint=(cutOffLength*i);
    cutOfMidline=midLineRefined(startpoint:endPoint);
    [m,n]=size(cutOfMidline);
    x=[1:n];
    x=x';
    cutOfMidline=cutOfMidline';
    f=fit(x,cutOfMidline,'poly4');
    y=feval(f,x);
    gaussStartPoint=y(1);
    if gaussStartPoint>nextStartPixel
        difference=gaussStartPoint-nextStartPixel;
        y=y-difference;
    elseif gaussStartPoint<nextStartPixel
        difference=nextStartPixel-gaussStartPoint;
        y=y+difference;
    end
    gaussLine(startpoint:endPoint)=y;
    nextStartPixel=gaussLine(endPoint);
end
minVal=min(gaussLine);
if minVal<0
PositiveMinVal=(round(sqrt(minVal.*minVal)))+1;
gaussLine=gaussLine+PositiveMinVal;
end


[q,r]=size(gaussLine);
[m,n]=size(cem);
yStartPoint=(n-r)-1;
yEndPoint=r;
gaussLineImg=zeros(m,n);
[row,col]=find(gaussLine);
endPoint=max(col);
for i=1:endPoint
midlinePoint=yStartPoint+i;
midLineVal=round(gaussLine(i));
if midLineVal<1
midLineVal=round(gaussLine(i-1));
if midLineVal<1
    midLineVal=round(gaussLine(i-2));
end
end
if midLineVal>5
    gaussLineImg((midLineVal-5):(midLineVal+5),i)=1;
elseif midLineVal<6
    gaussLineImg((midLineVal-(midLineVal-1)):(midLineVal+(midLineVal-1)),i)=1;
end
end

gaussLineBW=im2bw(gaussLineImg);
thirdIncrementGaussBW=gaussLineBW;
%finalThirdincrementTexture=immultiply(gaussLineBW,original);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=21:n
    colVec=fourthIncrementBW(:,i);
    [row,col]=find(colVec);
    firstPixel=min(row);
    lastPixel=max(row);
    if lastPixel>firstPixel
        midPixel=firstPixel+((lastPixel-firstPixel)/2);
        midLine(i)=midPixel;
    elseif lastPixel<(firstPixel+1)
        midPixel=firstPixel;
        midLine(i)=midPixel;
    end
end
[m,n]=size(cem);
[g,h]=size(midLine);
midLineRefined=midLine(21:h);
midLineRefined=midLineRefined(2:end);
[m,n]=size(midLineRefined);
cutOffLength=floor(n/5);
gaussLine=zeros(1,h);
startpoint=(cutOffLength*(1-1))+1;
endPoint=(cutOffLength*1);
cutOfMidline=midLineRefined(startpoint:endPoint);
[m,n]=size(cutOfMidline);
x=[1:n];
x=x';
cutOfMidline=cutOfMidline';
f=fit(x,cutOfMidline,'poly4');
y=feval(f,x);
gaussLine(startpoint:endPoint)=y;
nextStartPixel=gaussLine(endPoint);
for i=2:5
    startpoint=(cutOffLength*(i-1))+1;
    endPoint=(cutOffLength*i);
    cutOfMidline=midLineRefined(startpoint:endPoint);
    [m,n]=size(cutOfMidline);
    x=[1:n];
    x=x';
    cutOfMidline=cutOfMidline';
    f=fit(x,cutOfMidline,'poly4');
    y=feval(f,x);
    gaussStartPoint=y(1);
    if gaussStartPoint>nextStartPixel
        difference=gaussStartPoint-nextStartPixel;
        y=y-difference;
    elseif gaussStartPoint<nextStartPixel
        difference=nextStartPixel-gaussStartPoint;
        y=y+difference;
    end
    gaussLine(startpoint:endPoint)=y;
    nextStartPixel=gaussLine(endPoint);
end
minVal=min(gaussLine);
if minVal<0
PositiveMinVal=(round(sqrt(minVal.*minVal)))+1;
gaussLine=gaussLine+PositiveMinVal;
end


[q,r]=size(gaussLine);
[m,n]=size(cem);
yStartPoint=(n-r)-1;
yEndPoint=r;
gaussLineImg=zeros(m,n);
[row,col]=find(gaussLine);
endPoint=max(col);
for i=1:endPoint
midlinePoint=yStartPoint+i;
midLineVal=round(gaussLine(i));
if midLineVal<1
midLineVal=round(gaussLine(i-1));
if midLineVal<1
    midLineVal=round(gaussLine(i-2));
end
end
if midLineVal>5
    gaussLineImg((midLineVal-5):(midLineVal+5),i)=1;
elseif midLineVal<6
    gaussLineImg((midLineVal-(midLineVal-1)):(midLineVal+(midLineVal-1)),i)=1;
end
end

gaussLineBW=im2bw(gaussLineImg);
fourthIncrementGaussBW=gaussLineBW;
%finalFourthincrementTexture=immultiply(gaussLineBW,original);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[m,n]=size(cem);
[g,h]=size(midLine);
midLineRefined=midLine(21:h);
midLineRefined=midLineRefined(2:end);
[m,n]=size(midLineRefined);
    x=[1:n];
    x=x';
    midLineRefined=midLineRefined';
    f=fit(x,midLineRefined,'poly4');
    y=feval(f,x);
gaussLine=y;
[q,r]=size(gaussLine);
[m,n]=size(cem);
yStartPoint=(n-q)-1;
[m,n]=size(cem);
[g,h]=size(midLine);
midLineRefined=midLine(21:h);
midLineRefined=midLineRefined(2:end);
[m,n]=size(midLineRefined);
cutOffLength=floor(n/5);
gaussLine=zeros(1,h);
startpoint=(cutOffLength*(1-1))+1;
endPoint=(cutOffLength*1);
cutOfMidline=midLineRefined(startpoint:endPoint);
[m,n]=size(cutOfMidline);
x=[1:n];
x=x';
cutOfMidline=cutOfMidline';
f=fit(x,cutOfMidline,'poly4');
y=feval(f,x);
gaussLine(startpoint:endPoint)=y;
nextStartPixel=gaussLine(endPoint);
for i=2:5
    startpoint=(cutOffLength*(i-1))+1;
    endPoint=(cutOffLength*i);
    cutOfMidline=midLineRefined(startpoint:endPoint);
    [m,n]=size(cutOfMidline);
    x=[1:n];
    x=x';
    cutOfMidline=cutOfMidline';
    f=fit(x,cutOfMidline,'poly4');
    y=feval(f,x);
    gaussStartPoint=y(1);
    if gaussStartPoint>nextStartPixel
        difference=gaussStartPoint-nextStartPixel;
        y=y-difference;
    elseif gaussStartPoint<nextStartPixel
        difference=nextStartPixel-gaussStartPoint;
        y=y+difference;
    end
    gaussLine(startpoint:endPoint)=y;
    nextStartPixel=gaussLine(endPoint);
end
minVal=min(gaussLine);
if minVal<0
PositiveMinVal=(round(sqrt(minVal.*minVal)))+1;
gaussLine=gaussLine+PositiveMinVal;
end


[q,r]=size(gaussLine);
[m,n]=size(cem);
yStartPoint=(n-r)-1;
yEndPoint=r;
gaussLineImg=zeros(m,n);
[row,col]=find(gaussLine);
endPoint=max(col);
for i=1:endPoint
midlinePoint=yStartPoint+i;
midLineVal=round(gaussLine(i));
if midLineVal<1
midLineVal=round(gaussLine(i-1));
if midLineVal<1
    midLineVal=round(gaussLine(i-2));
end
end
if midLineVal>5
    gaussLineImg((midLineVal-5):(midLineVal+5),i)=1;
elseif midLineVal<6
    gaussLineImg((midLineVal-(midLineVal-1)):(midLineVal+(midLineVal-1)),i)=1;
end
end

gaussLineBW=im2bw(gaussLineImg);
fifthIncrementGaussBW=gaussLineBW;
%finalFifthincrementTexture=immultiply(gaussLineBW,original);
firstIncrementFileName=[sdirectory5 tifFiles2(j).name ' first increment.txt'];
secondIncrementFileName=[sdirectory5 tifFiles2(j).name ' second increment.txt'];
thirdIncrementFileName=[sdirectory5 tifFiles2(j).name ' third increment.txt'];
fourthIncrementFileName=[sdirectory5 tifFiles2(j).name ' fourth increment.txt'];
fifthIncrementFileName=[sdirectory5 tifFiles2(j).name ' fifth increment.txt'];

dlmwrite(firstIncrementFileName,firstIncrementGaussBW);
dlmwrite(secondIncrementFileName,secondIncrementGaussBW);
dlmwrite(thirdIncrementFileName,thirdIncrementGaussBW);
dlmwrite(fourthIncrementFileName,fourthIncrementGaussBW);
dlmwrite(fifthIncrementFileName,fifthIncrementGaussBW);

end