Matlab Process Images

hasDeepLearning = ~isempty(ver("nnet"));
hasComputerVision = ~isempty(ver("vision"));

imds = imageDatastore("data/images/", ...
    IncludeSubfolders=true, ...
    LabelSource="foldernames");  % Labels from folder names (e.g., cat/, dog/)

imds = imageDatastore("data/images/");

classNames = ["background", "tumor", "healthy"];
labelIDs = [0 1 2];
pxds = pixelLabelDatastore("data/labels/", classNames, labelIDs);

ds = combine(imds, pxds);

imds = imageDatastore("data/images/");

% Targets from a table or file (e.g., depth maps, angles, coordinates)
targets = readtable("data/targets.csv");  % columns: filename, value1, value2, ...
targetDs = arrayDatastore(targets{:, 2:end});

ds = combine(imds, targetDs);

inputImds = imageDatastore("data/noisy/");
targetImds = imageDatastore("data/clean/");
ds = combine(inputImds, targetImds);

targetSize = [256 256];  % [rows cols] — NOT [width height]

% For datastores — use transform
dsResized = transform(ds, @(data) resizeForNetwork(data, targetSize));

function out = resizeForNetwork(data, targetSize)
    out{1} = imresize(data{1}, targetSize);
    out{2} = imresize(data{2}, targetSize, "nearest");  % Labels: nearest!
end

function out = letterboxResize(img, targetSize)
    [h, w, ~] = size(img);
    scale = min(targetSize(1)/h, targetSize(2)/w);
    newSize = round([h w] * scale);
    resized = imresize(img, newSize);

    padValue = 114;  % common gray fill for detection networks
    out = ones(targetSize(1), targetSize(2), size(img,3), class(img)) * padValue;
    offset = floor((targetSize - newSize) / 2) + 1;
    out(offset(1):offset(1)+newSize(1)-1, offset(2):offset(2)+newSize(2)-1, :) = resized;
end

% ImageNet (ResNet, VGG, EfficientNet, etc.)
channelMeans = single([0.485, 0.456, 0.406]);
channelStds  = single([0.229, 0.224, 0.225]);

% COCO (detection/segmentation models pretrained on COCO)
channelMeans = single([0.471, 0.448, 0.408]);
channelStds  = single([0.234, 0.239, 0.242]);

% Apply: input must be single [0,1] first
img = im2single(img);  % uint8 [0,255] → single [0,1]
img = (img - reshape(channelMeans, 1, 1, 3)) ./ reshape(channelStds, 1, 1, 3);

% Step 1: Compute channel-wise statistics
imds = imageDatastore("train/", IncludeSubfolders=true);
pixelSum = single([0 0 0]);
pixelSqSum = single([0 0 0]);
nPixels = 0;

while hasdata(imds)
    img = single(read(imds));
    nPx = size(img, 1) * size(img, 2);
    pixelSum = pixelSum + sum(img, [1 2]);
    pixelSqSum = pixelSqSum + sum(img.^2, [1 2]);
    nPixels = nPixels + nPx;
end
reset(imds);

channelMeans = pixelSum / nPixels;
channelStds = sqrt(pixelSqSum / nPixels - channelMeans.^2);
save("normalization_params.mat", "channelMeans", "channelStds");

% Step 2: Apply in training pipeline
params = load("normalization_params.mat");
dsTrain = transform(dsTrain, @(data) normalizeForTraining(data, params));

function out = normalizeForTraining(data, params)
    img = single(data{1});
    img = (img - reshape(params.channelMeans, 1, 1, [])) ./ ...
          reshape(params.channelStds, 1, 1, []);
    out = {img, data{2}};  % {normalized image, label}
end

augmenter = imageDataAugmenter( ...
    RandXReflection=true, ...
    RandRotation=[-10 10], ...
    RandScale=[0.9 1.1]);

augDs = augmentedImageDatastore(targetSize, imds, ...
    DataAugmentation=augmenter);

ds = combine(imds, pxds);
dsAug = transform(ds, @augmentPair);

function out = augmentPair(data)
    img = data{1};
    lbl = data{2};
    % Apply same geometric transform to both
    if rand > 0.5
        img = fliplr(img);
        lbl = fliplr(lbl);
    end
    tform = randomAffine2d(Rotation=[-15 15], Scale=[0.8 1.2]);
    rout = affineOutputView(size(img), tform);
    img = imwarp(img, tform, OutputView=rout);
    lbl = imwarp(lbl, tform, OutputView=rout, Interpolation="nearest");
    out = {img, lbl};
end

layers = [
    imageInputLayer([256 256 3], Normalization="none")
    convolution2dLayer(3, 32, Padding="same")
    batchNormalizationLayer
    reluLayer
    maxPooling2dLayer(2, Stride=2)
    convolution2dLayer(3, 64, Padding="same")
    batchNormalizationLayer
    reluLayer
    globalAveragePooling2dLayer
    fullyConnectedLayer(numClasses)
    softmaxLayer
];
net = dlnetwork(layers);

net = dlnetwork;

% Main path
net = addLayers(net, [
    imageInputLayer([256 256 3], Name="input")
    convolution2dLayer(3, 32, Padding="same", Name="conv1")
    reluLayer(Name="relu1")
    ]);

% Branch path
net = addLayers(net, [
    convolution2dLayer(1, 32, Padding="same", Name="conv_branch")
    reluLayer(Name="relu_branch")
    ]);

% Merge
net = addLayers(net, concatenationLayer(3, 2, Name="concat"));

% Connect: branch takes same input, concat merges both paths
net = connectLayers(net, "input", "conv_branch");
net = connectLayers(net, "relu1", "concat/in1");
net = connectLayers(net, "relu_branch", "concat/in2");

[encoder, outputNames] = pretrainedEncoderNetwork("resnet18", 4);

inputSize = [256 256 1];  % [H W C]
numClasses = 3;

% U-Net
net = unet(inputSize, numClasses, EncoderDepth=4);

% 3-D U-Net (volumetric)
net = unet3d([64 64 64 1], numClasses, EncoderDepth=3);

% DeepLab v3+ with pretrained backbone
net = deeplabv3plus([256 256 3], numClasses, "resnet18");

% BiSeNet v2 — pretrained transfer learning (keeps pretrained weights)
net = bisenetv2(inputSize, numClasses);

% BiSeNet v2 — custom uninitialized network (no pretrained weights)
net = bisenetv2(inputSize, numClasses, NetworkType="custom");

net = unet([256 256 3], 2, EncoderDepth=4);  % Start with 2-class U-Net
% Replace final conv + softmax for single-channel regression output
net = replaceLayer(net, "encoderDecoderFinalConvLayer", ...
    convolution2dLayer(1, 1, Name="encoderDecoderFinalConvLayer"));
net = removeLayers(net, "FinalNetworkSoftmax-Layer");

[encoder, outputNames] = pretrainedEncoderNetwork("resnet18", 4);
net = encoderDecoderNetwork([256 256 3], encoder, outputNames, ...
    NumClasses=numClasses);

options = trainingOptions("adam", ...
    InitialLearnRate=1e-3, ...
    MaxEpochs=50, ...
    MiniBatchSize=8, ...
    Shuffle="every-epoch", ...
    ValidationData=valDs, ...
    ValidationFrequency=50, ...
    Plots="training-progress");

net = trainnet(dsTrain, net, "crossentropy", options);

% Single image
result = semanticseg(testImage, net);  % Returns categorical label map

% Batch via datastore
pxdsResults = semanticseg(testImds, net, WriteLocation="results/");

params = load("normalization_params.mat");
testImg = im2single(imread("test.png"));
testImg = (testImg - reshape(params.channelMeans, 1, 1, [])) ./ ...
          reshape(params.channelStds, 1, 1, []);
result = semanticseg(testImg, net);

bim = blockedImage(largeImage);

patchSize = [256 256];
overlap = [32 32];
result = apply(bim, @(block) semanticseg(block.Data, net), ...
    BlockSize=patchSize, ...
    BorderSize=overlap, ...
    PadPartialBlocks=true);

scores = minibatchpredict(net, dsTest);
YPred = scores2label(scores, classNames);

% ROC curve and AUC
rocObj = rocmetrics(YTest, scores, classNames);
plot(rocObj);

% Confusion matrix
confusionchart(YTest, YPred);

% Confusion matrix metrics
metrics = evaluateSemanticSegmentation(pxdsResults, pxdsGT);
disp(metrics.ClassMetrics);    % Per-class IoU, accuracy
disp(metrics.DataSetMetrics);  % Mean IoU, mean accuracy, weighted IoU

% Convert categorical to numeric labels
labelMap = uint8(result);

% Remove small regions (per-class)
for c = 1:numClasses
    classMask = labelMap == c;
    classMask = bwareaopen(classMask, minPixels);
    labelMap(labelMap == c & ~classMask) = 0;
end

% Fill holes in segmentation
for c = 1:numClasses
    classMask = labelMap == c;
    classMask = imfill(classMask, "holes");
    labelMap(classMask) = c;
end