Find the number of regions and calculate the area of each region

I used Gauss Random Field equation to create the bi-phase region, I want to find the number of regions separated by the equation and calculate the area of every single region of them. can anyone help? I attached the code and the image of the bi-phase field.

clear; format short e;
rng('shuffle');
%% Calculate
export = 1; % Whether or not to save data for ABAQUS
N = 1000; % number of waves to sum
mesh_size = 301; % coarseness of pattern
max_phi = 2*pi; % maximum phase shift
pct_stiff = 0.57; % fraction of material representing hard phase
lambda = 0.320; % characteristic wavelength (mm)
sample_size = 3.5 ; % characteristic length of volume (mm)
% Generate physical coordinates
x = linspace(-sample_size/2,sample_size/2,mesh_size);
y = linspace(-sample_size/2,sample_size/2,mesh_size);
% Create 3D coordinate matrices
[X, Y] = meshgrid(x, y);
phase = zeros(size(X));
for n=1:N
    % Generate random phase shift
    rand_phase = max_phi * rand(1,1);
    % Generate random vector
    rand_theta = 2 * pi * rand(1,1);
    rand_vec = [cos(rand_theta) sin(rand_theta)];
    % Dot direction vector with coordinate matrices
    vec = X*rand_vec(1) + Y*rand_vec(2);
    % Update the continuos phase field with new wave function
    phase = phase + sqrt(2/N)*cos(2*pi/lambda*vec + rand_phase);
end
%% Threshold
% Allocate space for new binary matrix
biphase = phase;
% Find average of all phase data
val_avg = mean(biphase,'all');
% Find standard deviation of all phase data
val_std = std(biphase,0,'all');
% Find the value at which the population splits into desired proportions
split = val_avg - sqrt(2) * erfinv((pct_stiff-0.5) * 2) * val_std;
% Apply threshold
biphase(phase>=split) = 1;
biphase(phase<split) = 0;
contourf(X,-Y,biphase)

Answer 1

[Hassaan]

clear; 
format short e;
rng('shuffle');

%% Parameters
export = 1; % Whether or not to save data for ABAQUS
N = 1000; % number of waves to sum
mesh_size = 301; % coarseness of pattern
max_phi = 2 * pi; % maximum phase shift
pct_stiff = 0.57; % fraction of material representing hard phase
lambda = 0.320; % characteristic wavelength (mm)
sample_size = 3.5; % characteristic length of volume (mm)

%% Generate physical coordinates
x = linspace(-sample_size / 2, sample_size / 2, mesh_size);
y = linspace(-sample_size / 2, sample_size / 2, mesh_size);
[X, Y] = meshgrid(x, y);
phase = zeros(size(X));

%% Generate Gaussian Random Field
for n = 1:N
    rand_phase = max_phi * rand(1,1);
    rand_theta = 2 * pi * rand(1,1);
    rand_vec = [cos(rand_theta) sin(rand_theta)];
    vec = X * rand_vec(1) + Y * rand_vec(2);
    phase = phase + sqrt(2/N) * cos(2 * pi / lambda * vec + rand_phase);
end

%% Threshold to create a bi-phase region
val_avg = mean(phase, 'all');
val_std = std(phase, 0, 'all');
split = val_avg - sqrt(2) * erfinv((pct_stiff - 0.5) * 2) * val_std;
biphase = double(phase >= split);

%% Find connected components
CC = bwconncomp(biphase);
numRegions = CC.NumObjects;
disp(['Number of Regions: ', num2str(numRegions)]);

%% Calculate the area of each region
pixelArea = (sample_size / mesh_size)^2; % Area represented by each pixel
areas = zeros(numRegions, 1);

for k = 1:numRegions
    numPixels = numel(CC.PixelIdxList{k});
    areas(k) = numPixels * pixelArea;
end

disp('Area of each region (mm^2):');
disp(areas);

%% Visualization
labelMatrix = labelmatrix(CC);
figure;
contourf(X, -Y, labelMatrix); % Using negative Y for correct orientation
colormap(jet(numRegions)); % Apply a colormap with a distinct color per region
colorbar;
title('Contour Plot of Regions');

[Huan]

Thank you very much for your answer. however, it seems it didn't find all the regions separated by the equation.