%% Sandbox Project
% 1) Organize Abnormal Image Data based on following classifications
    % Classfications
  % CIRC  Well-defined/circumscribed masses
  % SPIC  Spiculated masses
  % MISC  Other, ill-defined masses
  % ARCH  Architectural distortion
  % ASYM  Asymmetry
  % NORM  Normal
    % Background Tissue
  % F = Fatty
  % G = Fatty Glandular
  % D = Dense
    % Location and Size of Abnormality
  % x = xcoordinates
  % y = ycoordinates
  % r = radius
% 2) view how different gray level thresholds affect
%    mass detection
% 3) Test Images for which threshold value
%    yields correct location of abnormality
% 4) Export all data into the file: Mammogram Data.txt
% 5) Analyze the data in MammogramStats.m
close all
clear
clc
%% Getting Image Info 
% Opening the file "Mammogram Info With Extra.txt
% Extracting image information and reformatting to array form
fileID = fopen('Mammogram Info With Extra.txt');
C = textscan(fileID,'%3s %3d %s %4s %s %3d %3d %2d');
TextInfo = string([string(C{2}) C{3} C{4} C{5}]);
% order of columns: image number, tissue type, mass description, benign/malignant
TextLocations = double([C{2} C{6} C{7} C{8}]);
% order of columns: image number, xcoord, ycoord, radisu
fclose(fileID);
%% Importing Abnormal Images from Mammogram Database
% Accessing all .pgm files in the Mammogram Database
Dir = 'pwd';
GetDir = dir('/MATLAB Drive/Mammogram Database/*.pgm');
% Extracting xcoord, ycoord, r, and image num arrays.
% Finding the number of images in database
imageNum = nonzeros(TextLocations(:,1));
x = nonzeros(TextLocations(:,2));
y = nonzeros(TextLocations(:,3));
r = nonzeros(TextLocations(:,4));
numImages = length(imageNum);
%% Testing Abnormal Images
% initialize variable i, the size of numImages
% starting for loop
for i = 1:1:numImages
    close all % to close all images once the image is processed
    % Only abnormal mammograms are processed by the database
    % By finding finding consecutive imagenum from the imageNum array
    currentImage = imageNum(i);
    file = GetDir(currentImage).name;
    X = imread(file);
    figure;
axis on
%% Filtering Image for Better Contrast
% adapthisteq is used to improve contrast to 
% better visualize bright spots
figure;
enhanced_contrast = adapthisteq(X);
title('Enhanced Contrast Image')
imshow(enhanced_contrast);

%% Thresholding
% Varies thresholding of the image until
% the area and location of abnormality matches
% actual values
axis on
plotNumber = 0;
%% Using Varying Thresholds to find Masses
% Make another function to test classifier
% Use this file to find values
% figure out a way to organize the data
% need location, size, and threshold in one array
% initializing for loop counter "n"
n = 0;
TF = 0;
for k = .750:.025:.875
binaryimage = imbinarize(enhanced_contrast,k);
%% Using the Strel Function to remove noise
    SE = strel('disk',4);
    SE.Neighborhood
    % The imopen function removes any spots smaller
    % than a radius of 4 pixels. 
    % The function is applied to the negative of the filtered
    % image to remove black spots.
    open = imopen(~(binaryimage),SE);
    % The new image is shown on figure 2
    % as the original cropped, binarized, and filtered 
    % image for comparison
%% Using BW Label to find Abnormality
    [Labeled, number_blobs] = bwlabel(~open,8);
    % Use regionprops to find the centroid and boundary of the 
    % abnormality
    props = regionprops(~open);
    % find the areas of the regions in props
    areas = [props.Area];
    [sortedAreas,indicies] = sort(areas,'descend');
    % finding largest blob
    % This is most likely the abnormality
    % Includes coordinates of largest area (x,y), k (threshold), Area
    Identified = ismember(Labeled, indicies(sortedAreas>500));
    Spots = Identified>0;
%% Plotting Images of Varying Thresholds With Identified Masses
    figure(4);
    axis xy
    plotNumber = plotNumber +1;
    subplot(2,3,plotNumber)
    imshow(Spots)
    title('Identified Masses')
    hold on
    % use region props to plot the location 
    % of each identified mass
    % plot the largest identified mass
    n = n+1;
    Data = zeros();
    limit = 500;
    for idx = 1:length(areas)
        if areas(idx)>limit
        xcoordinate = props(idx).Centroid(1);
        ycoordinate = props(idx).Centroid(2);
        Data(idx,1) = xcoordinate;
        Data(idx,2) = ycoordinate;
        plot(xcoordinate,ycoordinate,'bx','DisplayName','Masses Identified by Classifier')
        end 
    end
%% Plotting actual region of abnormality
center = [x(i),1024-y(i)];
th = linspace(0,2*pi);
xcoord = center(1) + r(i)*cos(th);
ycoord = center(2) + r(i)*sin(th);
hold on
plot(xcoord, ycoord,'r','DisplayName','Diagnosed Mass')
%% Testing Found Areas to actual value
% Determine best threshold
% Each x and y coordinate is tested to determine
% if the values are within the radius of the known mass size
% the coordinates are saved if this is true
  [rows,columns] = size(Data);
  imagenumber = imageNum(i);
    LBx = x(i)-1.5*r(i);
    UBx = x(i)+1.5*r(i);
    LBy = y(i)-1.5*r(i);
    UBy = y(i)+1.5*r(i);
  if Data ==0
      fid2 = fopen('Mammogram Data.txt', 'a');
            fprintf(fid2,'%3d 0000 %s \n', imagenumber, TextInfo(i,2));
            fclose(fid2);
  else
    DataX = Data(:,1);
    DataY = Data(:,2);
    xc = DataX(DataX(:,1)~=0);
    yc = DataY(DataY(:,1)~=0);
    correctx = xc >= LBx & xc <= UBx;
    correcty = yc >= LBy & yc <= UBy;
    correct = correctx==1 & correcty==1;
    Thresholds = [];
    for s = 1:length(xc)
        % finds all thresholds for which the mass is recognized by the classifier
        % the values are printed in a separate excel file and text file
        if correct(s) == 1
           Thresholds(s) = k;
        end
    end
  end
    sortThresh= sort(Thresholds,'descend');
        if isempty(sortThresh) == 0
            min = sortThresh(1);
            fid = fopen('Mammogram Results.txt', 'a');
            fprintf(fid,'%.3f Threshold is Adequate for %s Classified Tissue in Image %3d \n', min, TextInfo(i,2),imagenumber);
            fclose(fid);
            fid2 = fopen('Mammogram Data.txt','a');
            fprintf(fid2,'%03d %.3f %s \n', imagenumber, min, TextInfo(i,2));
            fclose(fid2); 
        end
end

end
        
