clc; clear, close all;
images = {'book-cropped.jpg'};
for nImage = 1 : length (images)

% Load image
[pathStr, name, ext] = fileparts (images{nImage});
img = imread (images{nImage});

% Define parameters for local Otsu
winSize = 32;
stepSize = 16;
heightCrop = floor (size(img,1) / winSize) * winSize;
widthCrop = floor (size (img,2) / winSize) * winSize;
img = img (1 : heightCrop, 1: widthCrop, :);
imgGray = rgb2gray (img);
imgDouble = double (img);
imgGrayDouble = double (imgGray);

% Adaptive binarization
localThresh = 60;
[imgBW, localUsed, localThresh] = localAdaptiveOtsu ( imgGray, winSize, stepSize, localThresh );
imgBW = 1 - imgBW;
figure; clf;
set (gcf, 'Position', [150 150 300 200]);
subplot (2, 3, 1); imshow (imgBW); title ('Binarized Image');
subplot (2, 3, 2); imshow (localThresh, [0 255]); title ('Local Thresholds');
subplot (2, 3, 3); imshow (localUsed) ; title ('Local versus Global Thresholded Tiles');
imwrite (imgBw, [name '-binarize.jpg']);

%Keep eccentric blobs
imgBwLabel = bwlabel(imgBW);
shapeProps = regionprops (imgBWLabel, 'MajorAxisLength', 'MinorAxisLength');
for nRegion = 1 : length (shapeProps)
idx = find (imgBWLabel == nRegion);
eccentricity = shapeProps (nRegion) . MajorAxisLength / shapeProps (nRegion) .MinorAxisLength;
if eccentricity < 5
imgBW(idx) = 0;
end
end %nRegion
subplot (2, 3, 4); imshow (imgBW) ; title ('Eccentric Blobs');
imwrite (imgBW, [name '-binarize-eccentricity.jpg']);
% Find major axis length histogram
imgBWLabel = bwlabel (imgBW);
shapeProps = regionprops (imgBWLabel, 'MajorAxisLength');
majorLengths = zeros (1, length (shapeProps));
for nRegion = 1 : length (shapeProps)
majorLengths (nRegion) = shapeProps (nRegion) . MajorAxisLength;
end % nRegion
histBins = linspace (0, max (majorLengths));
majorHist = hist (majorLengths, histBins);
[maxHist, maxHistIdx] = max (majorHist);
maxMajorLength = histBins (maxHistIdx);

% Keep blobs with major axis length close to histogram-peak length
for nRegion = 1 : length (shapeProps)
	idx = find (imgBWLabel == nRegion);
	if abs (shapeProps (nRegion) .MajorAxisLength - maxMajorLength) > 0.4 * maxMajorLength
	imgBW (idx) = 0;
end
end % nRegion

subplot (2, 3, 5); imshow (imgBW); title ('Filtered by Major Axis Length');
imwrite (imgBW, [name '-binarize-length.jpg']);
% Find angle histogram
imgBWLabel = bwlabel (imgBW);
shapeProps = regionprops (imgBWLabel, 'Orientation');
orientations = zeros (1, length(shapeProps));
for nRegion = 1 : length (shapeProps)
	orientations (nRegion) = shapeProps (nRegion) .Orientation;
end
histBins = linspace (-90,90);
orientationHist = hist (orientation, histBins);
[maxHist, maxHistIdx] = max (orientationHist);
maxAngle = histBins (maxHistIdx);

% Keep blobs with orientation close to histogram- peak orientation
for nRegion = 1 : length (shapeProps)
	idx = find (imgBWLabel == nRegion);
	if abs (shapeProps (nRegion) .Orienatation - maxAngle) > 15
end
end

subplot (2, 3, 6); imshow (imgBW); title ('Filtered by Orientation');
imwrite (imgBW,[name '-binarize-angle.jpg']);

end