I'm trying to implement Otsu binarization technique on document images such as the one shown:

Could someone please tell me how to implement the code in MATLAB?

Taken from Otsu's method on Wikipedia

I = imread('cameraman.tif');

Step 1. Compute histogram and probabilities of each intensity level.

nbins = 256; % Number of bins
counts = imhist(I,nbins); % Each intensity increments the histogram from 0 to 255
p = counts / sum(counts); % Probabilities

Step 2. Set up initial omega_i(0) and mu_i(0)

omega1 = 0;
omega2 = 1;
mu1 = 0;
mu2 = mean(I(:));

Step 3. Step through all possible thresholds from 0 to maximum intensity (255)

Step 3.1 Update omega_i and mu_i

Step 3.2 Compute sigma_b_squared

for t = 1:nbins
    omega1(t) = sum(p(1:t));
    omega2(t) = sum(p(t+1:end));
    mu1(t) = sum(p(1:t).*(1:t)');
    mu2(t) = sum(p(t+1:end).*(t+1:nbins)');
end

sigma_b_squared_wiki = omega1 .* omega2 .* (mu2-mu1).^2; % Eq. (14)
sigma_b_squared_otsu = (mu1(end) .* omega1-mu1) .^2 ./(omega1 .* (1-omega1)); % Eq. (18)

Step 4 Desired threshold corresponds to the location of maximum of sigma_b_squared

[~,thres_level_wiki] = max(sigma_b_squared_wiki);
[~,thres_level_otsu] = max(sigma_b_squared_otsu);

There are some differences between the wiki-version eq. (14) in Otsu and the eq. (18), and I don't why. But the thres_level_otsu correspond to the MATLAB's implementation graythresh(I)