Implementing 2D Otsu Thresholding

I'm trying to binarize some very large (8548x10434) 2P images of whole mouse brains. I've been using Otsu (graythresh) with some degree of success but my images do not give bimodal histograms (image attached).

I think 2D Otsu will be an improvement but I'm unsure if it's the best method. I'm trying to use a 2D Otsu method I've come across here:

 %hists is a 256x256 2D-histogram of grayscale value and neighborhood average grayscale value pair.
%total is the number of pairs in the given image.
%threshold is the threshold obtained.
function threshold = Otsu2D(hists, total)
maximum = 0.0;
threshold = 0;
w_0=0;mu_ti=0;mu_tj=0;mu_i=0;mu_j=0;mu_is=0;mu_js=0;
for i =1:256
    for j=1:256
        P(i,j)=hists(i,j)/total;
    end
end
for i=1:256
    for j=1:256
        mu_ti=mu_ti+i*P(i,j);
        mu_tj=mu_tj+j*P(i,j);
    end
end
for ii=1:256
    for jj=1:256
        w_0=w_0+P(ii,jj);
        mu_is=mu_is+ii*P(ii,jj);
        mu_i(ii,jj)=mu_is;
        mu_js=mu_js+jj*P(ii,jj);
        mu_j(ii,jj)=mu_js;
        tr = ((mu_i(ii,jj)-(w_0*mu_ti))^2 + (mu_j(ii,jj)-(w_0*mu_tj))^2)/(w_0*(1-w_0));

        if ( tr >= maximum )
            threshold = ii;
            maximum = tr;
        end
    end
end

Im using the GLCM to compute the 2D histogram as per an example posted here previously by Image Analyst. The relevant code is:

%Now get the GLCM with graycomatrix().
kernel = [0, D; -D, D; -D, 0; -D, -D];
glcm = graycomatrix(grayImage, 'NumLevels', integerValue,...
  'GrayLimits', [],'offset',kernel);
numberOfDirections = size(glcm, 3);
for k = 1 : numberOfDirections
  glcm(:,:,k) % Display in the command window.
end
% Combine all directions into one.
% Sum along the third dimension.
glcm = sum(glcm, 3);
% Display the 2D glcm as a matrix.
subplot(1, 2, 2);
rgbImage = ind2rgb(glcm, hot(integerValue));
imshow(rgbImage);
% imshow(glcm, []);
axis on;
title('Image of the

 GLCM', 'FontSize', fontSize);

The problem is my threshold output from Otsu2D is always 0.