function WheatSpikelet_v1(D1,D2)
switch nargin%change D1 & D2 according to the data
  case 0
    D1 = 80;%D1, distance between spikelet (for single grain)
    D2 = 60;%D2, distance between spikelet (for grouped grains)
  case 1
    D2 = D1;
  case 2
  otherwise
    error('2 inputs are accepted.')
end
%% calculate grain traits
imageRes = 0.077;%image resolution, mm
% Load binary images to a 3d stack:ThreeD
listing = dir('*.tif');
im = imread(listing(1).name);
[a, b] = size(im);
ThreeD = zeros(a,b,length(listing));
sumThreeD = zeros(length(listing),1);
for t = 1:length(listing)
    temp = imread(listing(t).name);
    sumThreeD(t) = sum(sum(logical(temp)));
    ThreeD(:,:,t) = (temp);
end
%find the top and bottom of the wheat ear.
DDD = find(sumThreeD);
earBottom = max(DDD);
earTop = min(DDD);
earLength = (earBottom-earTop)*imageRes;%length of ear (in voxel).

% label the grains
CC = bwconncomp(ThreeD,6);
%calculate grain trais
stats = regionprops3(CC,"Volum",'SurfaceArea','PrincipalAxisLength','Orientation','Centroid');
grainVolume = stats.Volume.*imageRes^3;
grainnumber = size(stats,1);%number of grains
totalVolume = sum(grainVolume);
averagegrainVolume = totalVolume/grainnumber;
grainSurfaceArea = stats.SurfaceArea.*imageRes^2;
grainAxis = stats.PrincipalAxisLength.*imageRes;
grainCentroid = stats.Centroid.*imageRes;
graindata = [grainVolume,grainSurfaceArea, grainAxis,grainCentroid];
wheateardata = [earLength grainnumber totalVolume averagegrainVolume];
save('wheateardata_mm.txt','wheateardata','-ascii');
save('graindata_mm.txt','graindata','-ascii');
%-------------------------------------------------------
%% find spikelet
grainlocation = stats.Centroid;%get the 3d location of all the grains
grains = [1:grainnumber;1:grainnumber]';%creat a matrix to stroe the output
distmat = pdist(grainlocation);%calculate distance of each grain pairs
Z = squareform(distmat);%show the distance in 2D matrix
Z(Z<=0) = inf;% remove the 0 values, avoid problem in finding nearest neighbors
L = tril(Z); %lower left triangle of the distance matrix
% find nearest neighbor
for i = 1:grainnumber
    temp1 = find(Z(:,i) == min(Z(:,i)));
    if (Z(i,temp1)<D1)
        grains(i,2) = grains(temp1,1);
    end
end
% group to spikelets: find cluster
spikelet = zeros(grainnumber,7);
for i = 1:grainnumber
    if(~ismember(grains(i,1),spikelet))
        t = 1;
        spikelet(i,t) = grains(i,1);
        temp2 = grains(i,2);
        while(~ismember(temp2,spikelet(i,:)))
            t = t + 1;
            spikelet(i,t) = temp2;
            temp2 = grains(temp2,2);
        end
    end
end
% mearge spikelets with common grains
mergedspikelet1 = spikelet;
for i = 1:grainnumber
    [a,~] = find(mergedspikelet1 == i);
    gg = unique(mergedspikelet1(a,:));
    mm = size(gg(:),1)-1;% -1 is for the 0 value
    mergedspikelet1(a(1),1:mm) = gg(2:end);
    mergedspikelet1(a(2:end),:) = 0;
end
mergedspikelet1(mergedspikelet1(:,1) == 0,:) = [];
% merger mis-separated spikelets
mergedspikelet2 = mergedspikelet1;
for i = 1:(size(mergedspikelet2,1)-1)
    nog1 = sum(mergedspikelet2(i,:)>0);%
    nog2 = sum(mergedspikelet2(i+1,:)>0);
    if ((nog1+nog2 <= 5) && (nog1 > 0) && (nog2 > 0))
        grainID = [mergedspikelet2(i,1:nog1),mergedspikelet2(i+1,1:nog2)];
        grainDis = Z(mergedspikelet2(i,1:nog1),mergedspikelet2(i+1,1:nog2));
        if min(grainDis(:))<D2
            mergedspikelet2(i,1:nog1+nog2) = grainID;
            mergedspikelet2(i+1,:) = 0;
        end
    end
end
mergedspikelet2(mergedspikelet2(:,1) == 0,:) = [];
% save data to txt file
fileID = fopen('mergedspikelet_v1.txt','w');
fprintf(fileID, [repmat('%d ', 1, 7) '\n'], mergedspikelet2');
fclose(fileID);