Automatic Thresholding Using Second Derivative: The threshold is set as the x value of the maximum peak to the right of the max peak of second derivative of ksdensity of insulin pixels. Results: https://youtu.be/sptYuctKS5U

automatic_threshold_insulin.m

figure;
for n=1:size(img_names,1)
    clf
    % Load image
    img = imread([imgs_path img_names{n}]);
    insulin = imgaussfilt(double(img(:,:,2)),10);
    cyto = imgaussfilt(double(img(:,:,1)),1);
    nuc = imgaussfilt(double(img(:,:,3)),1);
    rgb = zeros(size(nuc,1),size(nuc,2),3);
    rgb(:,:,1) = cyto;
    rgb(:,:,2) = insulin;
    rgb(:,:,3) = nuc;
    
    % Display original color image
    subplot(2,2,1);
    imshow(uint8(rgb),[])
    title('Original Image')

    % Plot ksdensity of insulin values
    subplot(2,2,2);
    [f,xi] = ksdensity(insulin(:),'bandwidth',1);
    plot(xi,f)
    title('ksdensity of insulin pixels')
    spacing = xi(2)-xi(1); % Store spacing between sliding window jumps of ksdensity

    % Display second derivative of the ksdensity curve
    subplot(2,2,4);
    dx=gradient(f,spacing);
    dxx=gradient(dx,spacing);
    plot(xi,dxx)
    title('second derivative of the ksdensity curve')

    % Maximum peak to the right of the max peak of second derivative of ksdensity of insulin pixels
    [pks, locs] = findpeaks(dxx);
    locs = ((locs-1)*spacing)+min(xi); % Set correct range of x values. The derivative necessitates the -1.
    [max_peak,max_peak_loc]=max(f);
    max_peak_loc = ((max_peak_loc-1)*spacing)+min(xi);
    pks = pks(locs>max_peak_loc);
    locs = locs(locs>max_peak_loc);
    [max_peak,max_peak_loc]=max(pks);
    thresh = locs(max_peak_loc);

    % Plot red line to visualize threshold on second derivative plot
    hold on
    line([thresh thresh], [min(dxx) max(dxx)],'Color','red')

    % Plot red line to visualize threshold on ksdensity plot
    subplot(2,2,2);
    hold on
    line([thresh thresh], [min(f) max(f)],'Color','red')

    % Display thresholded insulin image
    subplot(2,2,3);
    imshow(insulin>thresh,[])
    title(['loop = ' num2str(n) ', thresh = ' num2str(thresh)])
    pause
end