function [blob, shapeID]= getimages(binaryImage, image)

[labeledImage, numberOfObjects] = bwlabel(binaryImage);
blobMeasurements = regionprops(labeledImage,...
    'Perimeter', 'FilledArea', 'BoundingBox');

blob = cell(1);
shapeID = length(numberOfObjects);
shapeID(1) = 0;

boundaries = bwboundaries(binaryImage);
hold on;
perimeters = [blobMeasurements.Perimeter];
filledAreas = [blobMeasurements.FilledArea];
% Calculate circularities:
circularities = perimeters .^2 ./ (4 * pi * filledAreas);

for blobNumber = 1 : numberOfObjects
    
    thisBoundary = boundaries{blobNumber};
    plot(thisBoundary(:,2), thisBoundary(:,1), 'r', 'LineWidth', 1);
    thisBlobsBoundingBox = blobMeasurements(blobNumber).BoundingBox;
    
    subBinaryImage = imcrop(binaryImage, thisBlobsBoundingBox);
    subImage = imcrop(image, thisBlobsBoundingBox);
    
    if circularities(blobNumber) < 1.2
        if iscirkle(subBinaryImage)
%             shape(blobNumber) = 'circle';
            shapeID(blobNumber) = 1;
            blob{blobNumber} = subImage;
        else
%             shape(blobNumber) = 'square < 1.2';
            shapeID(blobNumber) = 2;
            blob{blobNumber} = subImage;
        end
    elseif circularities(blobNumber) < 1.6
%         shape(blobNumber) = 'square < 1.6';
        shapeID(blobNumber) = 3;
        blob{blobNumber} = subImage;
    elseif circularities(blobNumber) > 1.6 && circularities(blobNumber) < 1.8
%         shape(blobNumber) = 'triangle';
        shapeID(blobNumber) = 4;
        blob{blobNumber} = subImage;
    else
%         shape(blobNumber) = 'something else';
        shapeID(blobNumber) = 5;
        blob{blobNumber} = subImage;
    end
%     disp(shape)
end