%% =========================
% [3.1] Load and display image
% =========================
imOriginal = imread("Section3.png");
figure;
imshow(imOriginal);
title("Original MRI Image");
% If image is RGB, convert to grayscale (usually this file is already grayscale)
if ndims(imOriginal) == 3
imOriginal = rgb2gray(imOriginal);
end
%% =========================
% [3.1] Remove salt-and-pepper noise using median filter
% =========================
imMed = medfilt2(imOriginal); % 2D median filter
figure;
imshowpair(imOriginal, imMed, "montage", "scaling", "none");
title("Original (left) vs Median Filtered (right)");
%% =========================
% [3.1] Compare with averaging filter (why median works better for S&P noise)
% =========================
avgFilt = (1/9) * ones(3,3);
% filter2 returns double, so cast back for display if needed
imAvg = filter2(avgFilt, imOriginal);
imAvg = uint8(imAvg);
figure;
imshowpair(imOriginal, imAvg, "montage", "scaling", "none");
title("Original (left) vs Averaging Filter (right)");
%% =========================
% [3.1] Histogram of median-filtered image
% Lab says use median-filtered output for intensity analysis
% =========================
figure;
imhist(imMed);
title("Histogram of Median-Filtered Image");
%% =========================
% [3.1] Create logical mask for NON-background pixels
% Black background pixel value = 0 (for uint8 grayscale)
% =========================
logicArr = (imMed ~= 0); % 1 = non-background, 0 = black background
%% =========================
% [3.1] Contrast adjustment using ONLY non-background pixels
% Steps:
% 1) extract non-background pixels
% 2) compute intensity limits from those pixels
% 3) use imadjust on the whole image with those limits
% =========================
nonBgPixels = imMed(logicArr);
% Simple version (directly from image values)
inLow = double(min(nonBgPixels)) / 255;
inHigh = double(max(nonBgPixels)) / 255;
% Optional better version (more robust to outliers):
% pr = prctile(nonBgPixels, [1 99]);
% inLow = double(pr(1))/255;
% inHigh = double(pr(2))/255;
% Contrast-adjust image
imAdj = imadjust(imMed, [inLow inHigh], []);
%% Histogram after adjustment (should use more of the intensity range)
figure;
imhist(imAdj);
title("Histogram After Contrast Adjustment");
%% =========================
% [3.1] Final comparison: original vs median-filtered + contrast-adjusted
% =========================
figure;
imshowpair(imOriginal, imAdj, "montage", "scaling", "none");
title("Original (left) vs Median Filtered + Contrast Adjusted (right)");
