Using ridge detection to characterise grains in an SEM image

I am trying to extract data about grains from an SEM image of a substance I am working on. I would like to extract features such as average grain size, eccentricity, angle, surface roughness etc. The image region analyser gives me all the data I need, but I am struggling to segment the image in a way that makes the information available. This is the image I am looking at:

I tried using an edge detect (after some sharpening and contrast adjustment), but none of the edges are outputted as continuous, and some of the gaps between grains are quite large. Simple thresholding also did not work since some of the grain boundaries have the same darkness as the cuts and scratches on the individual grains, leading to erroneous results.

My current attempt first smooths and sharpens the image (some of the images can be a bit noisy), performs an iterative ridge filter, which is then thresholded using the maximum gradient in the colour-vs-intensity graph as as the threshold point. This is then multiplied by the original image to preserve large dark areas and thresholded again. The issue I have is that the output is pretty good, but not perfect. There are tiny gaps in the output that mean that when the region analyser is applied, one very large grain is detected, as they are all connected:

If I modify the parameters to close these, the code becomes too harsh and I end up with hundreds of tiny grains:

I have tried using morphological functions to close these gaps, but I get a similar problem. Does anyone know of any techniques to deal with this problem, or any algorithms that could do this better? Thanks in advance.

Here is my code, sorry if it is a bit unclear, I have been tacking things on and deleting things a lot.

   %Parameters. Play around with these until you get what you want
  wienersize=3; %size of wiener filter
  averagesize=3; %size of averaging filter
  ridgeiterations=1; %no. of ridge detection iterations. Increase if image has clean but faint lines. Decrease if image has lots of deep, wide pores.
  mingauss=1; %minimum sigma for the gaussian filter within the ridge detection. Increasing cuts off visibility of the small scale variations.
  maxgauss=3; %maximum sigma for the gaussian filter within the ridge detection. Increasing this increases the scales that the ridge filter can see, but if gaussstep is too low it will introduce some weird artifacts
  gaussstep=1; %step size of gaussian filter in ridge detection
  sizeavg=5; %size of smoothing region for pixel graph
  halfmaxormaxgrad='maxgrad'; %controls whether thresholding uses half of the maximum value, or the maximum value of gradient. Good thing to change if you're getting a lot of speckles. Options are 'halfmax' or 'maxgrad'
  smallest=0.001; %controls how small the smallest region is. Is a fraction of the largest area.
  linelength=3; %line length of the morphological closing of the image
  noangles=4; %no. of angles of lines for the morphological closing

%Read and smooth the image to make it easier to characterise
close all
grainbounds=rgb2gray(imread('Grain boundary test 4.png'));
img = wiener2(grainbounds,[wienersize wienersize]);
img = uint8(filter2(fspecial('average',averagesize),img));
img = imsharpen(img);

%Ridge finding loop: performs a series of gaussian filters to find ridges
%on different scales. Finds the eigenvalues of the hessian matrix for each
%pixel and outputs the smaller one (not sure why this works so well, should
%output the difference, but it does)
valleys=zeros(size(img));
ridges=zeros(size(img));
tempimg=img;
for j=1:ridgeiterations
    for i=mingauss:gaussstep:maxgauss
        tempimg2 = imgaussfilt(tempimg,i);
        [Hx, Hy] = gradient(double(imcomplement(tempimg2)));
        [Hxx, Hxy] = gradient(Hx);
        [Hyx, Hyy] = gradient(Hy);
        valleys=valleys+0.5*(Hxx+Hyy-sqrt(Hxx.^2+4.*Hxy.^2-2.*Hxx.*Hyy+Hyy.^2));
        ridges=ridges+0.5*(Hxx+Hyy+sqrt(Hxx.^2+4.*Hxy.^2-2.*Hxx.*Hyy+Hyy.^2));
    end
    tempimg=valleys;
    subplot(ridgeiterations,1,j)
    imshowpair(tempimg,img,'montage')
    title(['Ridge filter after iteration ',num2str(j)])
end

%Perform a threshold at the point where the gradient in the intensity vs
%pixel no. graph is brightest

%Rescale the image to 8 bytes
ridge_output = cast((tempimg - min(tempimg(:))) *255/ (max(tempimg(:)) - min(tempimg(:))),'uint8');

%Create a vector of pixel no.s for each intensity, and get rid of the outer
%extremes(usually large peaks)
distribution=zeros(256,1);
for i=0:255
    distribution(i+1)=sum(ridge_output(:)==i);
end
cutoff=10;
usable=distribution(cutoff:255-cutoff);
hold on

%Smooth the graph and plot
for i=sizeavg+1:size(usable,1)-sizeavg
    tempusable=usable;
    usable(i)=sum(tempusable(i-sizeavg:i+sizeavg))/(2*sizeavg+1);
end
plot(usable)
[maxvalue maxlocation]=max(usable);
line([maxlocation maxlocation], [0 maxvalue])
[halfmaxheight halfmaxheightlocation]=min(abs(usable(1:maxlocation)-maxvalue/2));
line([halfmaxheightlocation halfmaxheightlocation], [0 halfmaxheight])

%Calculate gradient of the pixel graph and find the maximum
gradusable=zeros(size(usable,1)-2*sizeavg-2,1);
for i=1:size(gradusable,1)
    gradusable(i)=(usable(i+sizeavg+2)-usable(i+sizeavg))/2;
end
[maxgrad maxgradlocation]=max(gradusable);
line([maxgradlocation+sizeavg+2 maxgradlocation+sizeavg+2], [0 usable(maxgradlocation+sizeavg+2)])

%Shift the location of the maximum to the correct place in the original
%vector
if halfmaxormaxgrad == 'halfmax'
    lower_bound=halfmaxheightlocation+sizeavg+cutoff+2;
elseif halfmaxormaxgrad == 'maxgrad'
    lower_bound=maxgradlocation+sizeavg+cutoff+2;
else print("Input Error: please set halfmaxormaxgrad to 'maxgrad' or 'halfmax' in the parameters section")
end

%Make all the pixels with a value higher than the max gradient white and
%all the pixel with a value lower black
threshold_output=imcomplement(imcomplement(ridge_output).*cast(ones(size(ridge_output))-(ridge_output >= lower_bound),'like',ridge_output));
threshold_output=threshold_output.*uint8(threshold_output == 255);
figure, imshowpair(threshold_output,grainbounds,'montage')
title('Mask after threshold')

%multiply the output by the smoothed image to preserve dark regions
masked=(cast(threshold_output,'double')./max(cast(threshold_output,'double'))).*(cast(img,'double')./max(cast(img,'double')));
figure, imshowpair(masked,grainbounds,'montage')
title('Original image multiplied by mask (allows deep pores to be conserved)')

%Perform the image segmentation
final_binary = masked>=graythresh(masked);
figure, imshowpair(grainbounds,final_binary)

grain_data=struct2array(regionprops((final_binary),'Area'));
areas=zeros(max(grain_data),1);
for i=1:max(grain_data)
    areas(i)=sum(grain_data(:)==i);
end
[ans boundlocation] = find(areas<=smallest*max(areas),1);
% areasfit=fit(transpose(linspace(1,boundvalue,boundvalue)),areas(1:boundvalue),'power1','StartPoint',[max(areas),log(areas(boundlocation))/log(max(areas)*boundlocation)])

%morphological operations

   final_binary=bwareaopen(final_binary,ans);
   figure, imshow(final_binary)
  final_binary=imcomplement(final_binary);
  for i=1:noangles
      lineel=strel('line',linelength,i*180/noangles);
      final_binary=imclose(final_binary,lineel);
  end
  final_binary=imcomplement(final_binary);

grain_data=struct2table(regionprops(final_binary,'Area','Centroid','Eccentricity','Orientation','MajorAxisLength','MinorAxisLength','Perimeter'));
areaavg=sum(grain_data.Area)./size(grain_data.Area,1);
areastd=sqrt(sum(grain_data.Area.^2)./size(grain_data.Area,1)-(sum(grain_data.Area)./size(grain_data.Area,1))^2);
areatot=sum(grain_data.Area);