main.m

%% Main Metrics Code
% Comment out these lines if using the GUI
% clear all
% clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 1: Converts the  to binary
%this section loads in
%loads the image and mask, Creates the brain and mask matrix, Information about files

% These lines are for testing without the GUI
%    imfile = ('Dicom.nii.gz');  
%    imbrain = (niftiread(imfile)); 
%    brain_info = niftiinfo(imfile);
% 
%    maskfile = ('Mask.nii.gz');  
%    immask = (niftiread(maskfile));    
%    mask_info = niftiinfo(maskfile);

% These lines are for using the GUI
     brain_info = niftiinfo(imbrainf);
     mask_info = niftiinfo(immaskf);

%gets info about sizes of image    

    imagesize = size(imbrain);
    r = imagesize(1,3);
    m = imagesize(1,1);   
    n = imagesize(1,2);

%Initalize variables to count # of slices that dont have aneurysms and counts slices with the aneurysms in it
    non_ia_slices = 0;
    num_ia_slices = 0;
     
for i = 1:r
%makes the image double so that it can be
    mask_slice = -double(immask(:,:,i));  
    brain(:,:,i) = double(imbrain(:,:,i));
   
    %goes through slices to make them 0's and 1's
         for l = 1:m
            for j = 1:n
                if mask_slice(l,j) ==  32768 %Check other images for correct thresh
                   binary_mask_i(l,j) = 0;
                else
                   binary_mask_i(l,j) = 1;
                end
            end
         end
         
      binary_mask(:,:,i) = binary_mask_i;
      nnz__mask = nnz(binary_mask_i);%number of non-zeros in the slice
%if there are no non-zero pixels then there is no aneurysm array that keeps track ofs what slices the aneurysm is in
        if nnz__mask == 0
          non_ia_slices =  non_ia_slices +1;
        else
          num_ia_slices = num_ia_slices +1;
          SLICE_NUMBER(num_ia_slices) = i;
        end
end
% Pass in the array of slices that contain the aneurysms
    slices = length(SLICE_NUMBER); %can edit this number to have slices before and after anuerysm
    SLICE_NUMBER_MIN = min(SLICE_NUMBER)-1; %can edit this number to have slices before and after anuerysm
    firstANSlice = min(SLICE_NUMBER);
%the code below allows the aneurysm slices to be extracted and put into one mXnXi 3D image
    slice = 0; %code debugger
    post_size_x= 0; %allows the size of an array to increase
    count = 1; %code debugger to count how many aneurysm pixels in one slice
    number_of_total_aneurysm_pixels = 1; %ocunts total number of pixels that exist
    newx = 0;
    newy = 0;
    per_slice_count = 0;
   
for i=1:slices
    ia_slice = SLICE_NUMBER_MIN+i;%chooses one slice in our 3D set is analyzed and extracts i
    binary_mask_with_aneurysm(:,:,i) = binary_mask(:,:,ia_slice); %saves all the slices of the mask with the aneurysm in i
    brain_with_aneurysm(:,:,i) = brain(:,:,ia_slice);
    [x,y]= find(binary_mask_with_aneurysm(:,:,i));

%if there is a value in x then the if statement will run to save all of the
%informations of the locations of the aneurysm in various arrays
            %counts how many anuerysm pixels exist on the image            
            %adds #pixels to the previous number to find the exact number
            %of pixels in all of the slices that hold the aneurysm
       
         if x >= 1
            size_x = length(x);
            size_y = length(y);
            newx = size_x + newx;
            newy = size_y + newy;
           
                for a = 1:size_x
                    locations_x(count) = x(a,1); %x locations of aneurysms in each slice
                    locations_y(count) = y(a,1); %y locations of aneurysms in each slice
                    pixel_value(count) = brain_with_aneurysm(x(a,1),y(a,1)); %pixel values of all the aneurysm locations
                    count = count+1;
                    per_slice_count = per_slice_count+1;  
                end
               
            slice = 1+slice;
            pixel_per_slice(slice) = per_slice_count; %number of pixels that make up the aneurysm per slice
            per_slice_count = 0;  
           
         end
   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%% Section 2: Thresholding the brain image
%for loop to make all of the slices binary that are before
    threshold = median(pixel_value);
%the mode, max, and min calculations below are used to extract the loaction
%of the anuerysm in hopes of finding the blood vessel attached to it
    mode_x = mode(locations_x);
    mode_y = mode(locations_y);
    max_x = max(locations_x);
    max_y = max(locations_y);
    min_x = min(locations_x);
    min_y = min(locations_y);
%the for loop below thresholds the brain image so that it becomes binary
for i = 1:slices
    braini = i+SLICE_NUMBER_MIN;
    for a = 1:m   %make universal by taking the size of the image eventually
        for b = 1:n
            if brain_with_aneurysm(a,b,i) >= threshold %unsure of what the trheshold should actually be. this needs to be more looked at
                brain_with_aneurysm_binary(a,b,i) = 1;
            else
                brain_with_aneurysm_binary(a,b,i) = 0;
            end
        end
    end  
     
end

number = 0; %code debugger
BVC = 0;
%this for loop allows the blood vessel to be extracted
for i = 1:slices
       brain = brain_with_aneurysm_binary(:,:,i);
%finds and labels all of the connected components with a specific pixel value, unless there are no
%connected components in which this function will leave the pixel value as zero      
       label = bwlabel(brain);
% % % % %        s = regionprops(label, 'Centroid');
% % % % %        for k = 1:numel(s)
% % % % %             c = s(k).Centroid;
% % % % %             text(c(1),c(2), sprintf('%d', k),...
% % % % %             'HorizontalAlignment', 'center',...
% % % % %             'VerticalAlignment', 'middle');
% % % % %         end
% % % % %         hold off
% % % % %         number = fprintf('Number of Grains of Rice = %d\n',k);
            
       val = label(min_x, min_y); %determines the pixel values conneccted to the aneurysm
%if the anuerysm is not on the blood vessel then the pixel value should have been 0, ...
%so if this occurs theis while loop will keep searching until it finds a pixel value that is not zero 
if BVC == 0 
        while  val == 0
            number = number + 1;
            xplustest = mode_x + number;
            yplustest = mode_y + number; 
            
            if xplustest > m || yplustest> n
                number = 0;
                
                while val == 0;
                    number = number +1;
                    xminustest = max_x - number;
                    yminustest = max_y - number;
                    val = label(xminustest, yminustest);
                end 
                
            else               
                val = label(xplustest, yplustest);
            end 
        
        end
        display(val);
else
    while val == 0
        number = number + 1;
        val = label(floor(s.Centroid(1,2)+number), floor(s.Centroid(1,1)+number));
    end
             
end 
%on the brain_with_aneurysm image it keeps only the components connected to the anuerysm
        for a = 1:m
               for b = 1:n
                    if val == label(a,b)
                       blood_vessel(a,b,i) = 1;
                   else
                       blood_vessel(a,b,i) = 0;
                   end
               end
        end
        s = regionprops(blood_vessel(:,:,1), 'Centroid');
        BVC = 1;
       % imshow(blood_vessel(:,:,i)); % commented out for Audi Presentation
        number = 0;     
end

%the below for loop overlays the aneurysm mask and the blood vessel and
%turns it into grayscale so it can all be saved as one mxnxi matrix
for s = 1: slices
    v = SLICE_NUMBER_MIN+s;
    mask  = binary_mask_with_aneurysm(:,:,s);
    blood_V = blood_vessel(:,:,s);
    over_layed = imfuse(mask, blood_V);
    overlayed_images(:,:,s) = rgb2gray(over_layed);
    BV_AN_Brain = imfuse(overlayed_images(:,:,s), imbrain(:,:,v)); % put this in the GUI
    four_d(:,:,:,s) = BV_AN_Brain;%outputs a 4D image such that we can use to try and put it through a 3D viewer and have it look mildly normal
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Uses function "Metrics" to analyze the ellipse surrounding the blood vessel
%This for loop is to analyze all of the pixel values in one slice to find
%the values that make up the blood vessel versus the values that make up
%the the aneurysm
for s = 1:slices
    for a= 1:m
        for b = 1:n  
               if overlayed_images(a,b,s) ~= 0
                   value = overlayed_images(a,b,s);
               end
        end
    end
    mode_bood_vessel_pixel_value = mode(value);
%The if statement below says:
%if the pixel value is above zero, the pixels will be analzyed from there
    for a= 1:m
            for b = 1:n  
                    if overlayed_images(a,b,s) == mode_bood_vessel_pixel_value
                          just_BV(a,b,s) = 1;   %record value at a,b,s in overlayed images
                    else
                          just_BV(a,b,s) = 0;
                    end
            end
    end
   
        for a= 1:m
            for b = 1:n  
                    if overlayed_images(a,b,s) ~= mode_bood_vessel_pixel_value && overlayed_images(a,b,s) ~= 0
                          just_AN(a,b,s) = 1;   %record value at a,b,s in overlayed images
                    else
                          just_AN(a,b,s) = 0;
                    end
            end
        end
end
PreviousLengthAN = 0;
PreviousLengthBV = 0;
for s = 1:slices

    %Blood Vessel calculations
        [MajorAxisBV, MinorAxisBV, OrientationBV, xRotatedBV, yRotatedBV, CenterBV] = metrics(just_BV(:,:,s));
    %calculation to break up the major axis into x and y components to find the dimensions of the major axis in mm
        verticalBV   = MajorAxisBV* (sin(OrientationBV)) * brain_info.PixelDimensions(1,1);
        horizontalBV = MajorAxisBV * (cos(OrientationBV)) * brain_info.PixelDimensions(1,2);
        subtractedBV = MajorAxisBV - MinorAxisBV;
        LengthBV(s) = sqrt(verticalBV^2 + horizontalBV^2); %pythagorean theorem  
       
%         fprintf("\\nBV: Slice Number:%f \nLength in mm:%f \\n",s, LengthBV, subtractedBV);
%         figure
%         imshow(just_BV(:,:,s))
%         hold on
%         plot(xRotatedBV, yRotatedBV, 'LineWidth', 2);
%         plot(CenterBV(1,1), CenterBV(1,2),'*');
%         hold off

    %Aneurysm Cacluations
        [MajorAxisAN, MinorAxisAN, OrientationAN, xRotatedAN, yRotatedAN, CenterAN] = metrics(just_AN(:,:,s));
    %Calculation to break up the major axis into x and y components to find the dimensions of the major axis in mm
        verticalAN   = MajorAxisAN* (sin(OrientationAN)) * mask_info.PixelDimensions(1,1);
        horizontalAN = MajorAxisAN * (cos(OrientationAN)) * mask_info.PixelDimensions(1,2);
        subtractedAN = MajorAxisAN - MinorAxisAN;
        LengthAN(s) = sqrt(verticalAN^2 + horizontalAN^2);%pythagorean theorem
       
%         fprintf("\\nAN: Slice Number:%f \nLength in mm:%f \\n ",s, LengthAN, subtractedAN);
%         %dislay the images
%         figure
%         imshow(just_AN(:,:,s))
%         hold on
%         plot(xRotatedAN, yRotatedAN, 'LineWidth', 2);
%         plot(CenterAN(1,1), CenterAN(1,2),'*');
%         hold off

end
%Largest Diameter:
BVLargestLength = max(LengthBV);
    largestBVs = find(LengthBV == BVLargestLength);
    largestBVs = largestBVs(1,1);
ANLargestLength = max(LengthAN);
    largestANs = find(LengthAN == ANLargestLength);
    largestANs = largestANs(1,1);
   
% Size Ratio
SizeRatio = ANLargestLength/BVLargestLength;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Putting Metrics in GUI
GUIBVLargestLength = num2str(BVLargestLength);
GUIANLargestLength = num2str(ANLargestLength);
GUIsizeRatio = num2str(SizeRatio);

if largestBVs == largestANs
    v = s+largestBVs;
    BV_AN_Brain = imfuse(overlayed_images(:,:,largestBVs), imbrain(:,:,v));
else
    v = firstANSlice-1 + largestANs;
    BV_AN_Brain = imfuse(overlayed_images(:,:,largestANs), imbrain(:,:,v));
end

%cropping the picture:
beginningy = floor(CenterAN(1,1))- 10;
endy = floor(CenterAN(1,1)) + 10;
beginningx = floor(CenterAN(1,2)) - 10;
endx = floor( CenterAN(1,2)) + 10;

BV_AN_Zoomed = BV_AN_Brain(beginningx:endx, beginningy:endy, :);