Matlab Process Images

function fileType = triageImageFile(filename)
    % Determine file type and recommend read path
    %
    % Returns: "dicom", "raw", or "standard"

    % 1. Check for DICOM (files often have no extension)
    try
        info = dicominfo(filename);
        fprintf("DICOM file detected — Modality: %s, Patient: %s\n", ...
            info.Modality, info.PatientName.FamilyName);
        fileType = "dicom";
        return
    catch
        % Not DICOM — continue
    end

    % 2. Check for camera RAW by extension
    [~, ~, ext] = fileparts(filename);
    rawExts = [".nef",".cr2",".cr3",".crw",".arw",".dng",".orf",".rw2",".raf",".srw",".pef"];
    if any(strcmpi(ext, rawExts))
        info = rawinfo(filename);
        fprintf("Camera RAW file — Sensor: %s, CFA size: %s\n", ...
            info.CFASensorType, mat2str(info.CFAImageSize));
        fileType = "raw";
        return
    end

    % 3. Standard image format
    fileType = "standard";
end

files = dir(fullfile(folder, "*"));
files = files(~[files.isdir]);

isDicom = false(numel(files), 1);
for i = 1:numel(files)
    try
        dicominfo(fullfile(folder, files(i).name));
        isDicom(i) = true;
    catch
        % Not DICOM
    end
end

if any(isDicom)
    fprintf("%d of %d files are DICOM — use medicalVolume or dicomreadVolume\n", ...
        nnz(isDicom), numel(files));
end

% File-level metadata (no pixel loading)
info = imfinfo("image.tif");
fprintf("Format: %s\n", info.Format);
fprintf("Dimensions: %d x %d\n", info.Width, info.Height);
fprintf("Color type: %s\n", info.ColorType);   % "truecolor", "grayscale", "indexed"
fprintf("Bit depth: %d\n", info.BitDepth);
fprintf("File size: %.1f KB\n", info.FileSize/1024);

% STOP: Check dimensions before loading pixels
if info.Width * info.Height > 100e6
    % Large image (>100 MP) — see references/large-image-processing.md
    return
end

% Pixel-level characterization (requires loading)
img = imread("image.tif");
fprintf("Class: %s\n", class(img));
fprintf("Size: %s\n", mat2str(size(img)));
fprintf("Data range: [%g, %g]\n", min(img(:)), max(img(:)));

info = imfinfo(filename);
if strcmp(info.ColorType, "indexed")
    [X, map] = imread(filename);  % X = index matrix, map = Nx3 colormap
    fprintf("Palette size: %d colors\n", size(map, 1));
    % Convert to RGB for processing
    rgbImg = ind2rgb(X, map);     % Returns double [0,1]
end

img = imread(filename);
if isa(img, "uint16")
    maxVal = double(max(img(:)));
    sensorBits = [10 12 14 16];
    effectiveBits = sensorBits(find(maxVal <= 2.^sensorBits - 1, 1));
    fprintf("Stored as 16-bit, effective depth: %d-bit (max value: %d)\n", ...
        effectiveBits, maxVal);
end

imageshow(img, DisplayRangeMode="12-bit");   % Range [0, 4095]
imageshow(img, DisplayRangeMode="10-bit");   % Range [0, 1023]
imageshow(img, DisplayRangeMode="data-range"); % Range [min(img), max(img)]

% PNG: alpha is the third output of imread
info = imfinfo(filename);
hasAlpha = ~strcmp(info.Transparency, "none");
[img, ~, alpha] = imread(filename);  % alpha is non-empty if PNG has transparency

% TIFF: alpha is packed into the image array (M×N×4), third output is always empty
% SamplesPerPixel == 4 could be RGBA or CMYK — check ColorType to distinguish
info = imfinfo(filename);
if info.SamplesPerPixel == 4 && strcmp(info.ColorType, "truecolor")
    data = imread(filename);          % M×N×4
    img   = data(:,:,1:3);            % RGB
    alpha = data(:,:,4);              % Alpha
end

if size(img, 3) == 3
    if isequal(img(:,:,1), img(:,:,2)) && isequal(img(:,:,2), img(:,:,3))
        img = img(:,:,1);  % Collapse to true grayscale
        fprintf("Converted gray-as-RGB to grayscale (3x memory saved)\n");
    end
end

dataFolder = "path/to/Data";
imds = imageDatastore(dataFolder, IncludeSubfolders=true, ...
    LabelSource="foldernames");

%% 1. Dataset overview
fprintf("=== Dataset Overview ===\n");
fprintf("Total images: %d\n", numel(imds.Files));
fprintf("Classes:\n");
disp(countEachLabel(imds));

%% 2. File format distribution
exts = strings(numel(imds.Files), 1);
for i = 1:numel(imds.Files)
    [~, ~, exts(i)] = fileparts(imds.Files{i});
end
[uniqueExts, ~, idx] = unique(exts);
counts = accumarray(idx, 1);
fprintf("Formats: ");
for i = 1:numel(uniqueExts)
    fprintf("%s (%d)  ", uniqueExts(i), counts(i));
end
fprintf("\n");

%% 3. Sample characterization — read a few to understand properties
numSamples = min(10, numel(imds.Files));
sampleIdx = round(linspace(1, numel(imds.Files), numSamples));
fprintf("\n=== Sample Characterization ===\n");
for i = sampleIdx
    info = imfinfo(imds.Files{i});
    info = info(1);
    fprintf("  %s: %dx%d, %s, %d-bit\n", ...
        info.Format, info.Width, info.Height, info.ColorType, info.BitDepth);
end

%% 4. Pixel statistics from samples
fprintf("\n=== Pixel Statistics ===\n");
for i = sampleIdx(1:min(3,end))
    img = imread(imds.Files{i});
    fprintf("  %s — class: %s, range: [%g, %g]\n", ...
        imds.Files{i}, class(img), min(img(:)), max(img(:)));
end

%% 5. Visual overview — show sample grid
montage(imds, Indices=sampleIdx, ThumbnailSize=[128 128]);

% Quick corruption check
badFiles = {};
for i = 1:numel(imds.Files)
    try
        imread(imds.Files{i});
    catch
        badFiles{end+1} = imds.Files{i}; %#ok<SAGROW>
    end
end
if ~isempty(badFiles)
    fprintf("Found %d corrupt files\n", numel(badFiles));
end

gray = im2single(im2gray(img));

if ~isempty(ver("wavelet"))
    % Wavelet MAD estimator (more accurate, requires Wavelet Toolbox)
    [~, ~, ~, d] = dwt2(gray, "haar");         % Single-level 2-D DWT, diagonal detail
    noiseStd = median(abs(d(:))) / 0.6745;     % Robust MAD estimator
else
    % Laplacian-based estimator (Immerkaer 1996, IPT only)
    kern = [1 -2 1; -2 4 -2; 1 -2 1];
    filtered = imfilter(gray, kern, "symmetric");
    [H, W] = size(gray);
    noiseStd = sqrt(0.5*pi) / (6*W*H) * sum(abs(filtered(:)));
end
fprintf("Estimated noise std: %.4f\n", noiseStd);

gray = im2gray(img);
[counts, binLocs] = imhist(gray);

% Compute number of significant peaks (modes)
smoothCounts = smoothdata(double(counts), "gaussian", 7);  % 1-D smoothing
peaks = islocalmax(smoothCounts, MinProminence=max(smoothCounts)*0.05);
numModes = nnz(peaks);
fprintf("Histogram modes: %d\n", numModes);

gray = im2double(im2gray(img));
lap = imfilter(gray, fspecial("laplacian", 0), "replicate");
sharpness = var(lap(:));
fprintf("Sharpness (Laplacian variance): %.6f\n", sharpness);

gray = im2double(im2gray(img));

% Estimate background with large Gaussian
background = imgaussfilt(gray, 50);

% Measure variation in background
bgRange = max(background(:)) - min(background(:));
bgStd = std(background(:));
fprintf("Background range: %.3f, std: %.3f\n", bgRange, bgStd);

% Quick binary estimate
bw = imbinarize(im2gray(img));
bw = bwareaopen(bw, 20);  % Remove tiny specks
cc = bwconncomp(bw);

fprintf("Estimated objects: %d\n", cc.NumObjects);

if cc.NumObjects > 0
    props = regionprops("table", cc, "Area", "Eccentricity", "Solidity");
    fprintf("Area range: [%d, %d] pixels\n", min(props.Area), max(props.Area));
    fprintf("Median eccentricity: %.2f\n", median(props.Eccentricity));
    fprintf("Median solidity: %.2f\n", median(props.Solidity));

    % Detect touching objects (low solidity + large area suggests merged blobs)
    if any(props.Solidity < 0.7 & props.Area > median(props.Area)*3)
        fprintf("Warning: likely touching/overlapping objects — consider watershed\n");
    end
end

gray = im2gray(img);
[h, w, c] = size(img);

% Compute diagnostic features
edgeDensity = nnz(edge(gray, "Canny")) / numel(gray);
if isinteger(gray)
    uniqueColors = numel(unique(gray(:))) / double(intmax(class(gray)));
    entropy_val = entropy(gray);
else
    uniqueColors = numel(unique(gray(:))) / numel(gray);
    entropy_val = entropy(im2uint8(mat2gray(gray)));
end

fprintf("Edge density: %.3f\n", edgeDensity);
fprintf("Color utilization: %.1f%%\n", uniqueColors * 100);
fprintf("Entropy: %.2f\n", entropy_val);

fprintf("\n=== Summary ===\n");
fprintf("Format: %s, %dx%d, %s, %d-bit\n", format, width, height, colorType, bitDepth);
fprintf("Data class: %s, range: [%g, %g]\n", dataClass, minVal, maxVal);
fprintf("Sharpness: %.5f | BG uniformity: %.3f\n", sharpness, bgStd);
fprintf("Histogram modes: %d | Objects detected: %d\n", numModes, numObjects);

% Actionable recommendations
if sharpness < 0.001, fprintf("→ Image is blurred — consider deblurring\n"); end
if bgStd > 0.05,   fprintf("→ Uneven illumination — apply imflatfield\n"); end
if numModes == 2,   fprintf("→ Bimodal histogram — imbinarize should work well\n"); end