Inflection points on grayscale image histogram

I have a grayscale image, then plotted the histogram using imhist. Now, I need to find the first inflection point of the histogram. That point halfway up, when the bell-shape curve starts to change (see red line in attached jpg). I calculated first and second derivative of the histogram, but I'm stuck on how to get that point. My code looks something like this:

I = imread('dog3.jpg');
%figure; imshow(I)
%Calculating and plotting histogram of image
y = imhist(I);
f1 = diff(y);
f2 = diff(f1);
figure; imhist(I)
figure; plot(y)
figure; plot(f2)

Answer 1

[Star Strider]

Use the gradient (link) function to calculate the numerical derivative, not the diff function.

See if the approach in my Answer to How to draw tangent line at infleciton point? (link) will work for you.

[Veronica Morales]

haha I am actually trying to use your approach, but I don't have an equation. I perform the gradient or diff on the y = imhist(I) variable. I tried the interp1 steps and got NaN.

[Star Strider]

An equation is not necessary for my code.. The gradient function calculates the numerical derivative.

If you got NaN values from the interp1 calls, tell interp1 to extrapolate.

For example:

yi = interp1(x, y, xi, 'linear','extrap');

That should eliminate the NaN results.

I don’t have your data, so I can’t write specific code for it.

[Veronica Morales]

Thank you so much, yeah eliminated the NaN! Here's my image. (very generic)

[Star Strider]

As always, my pleasure!

[Image Analyst]

Why are you thresholding that image? What do you want to find? And why do you think the inflection point is the right threshold as compared to anything else, like the full-width-half-max, the bottom corner (like from a triangle threshold which would be a better one to detect the background dark noise), the Otsu threshold, or simply doing a manual/interactive threshold like with my visual thresholding app: Thresholding applet? What's so special about that gray level?

[Star Strider]

@ Veronica Morales —

The findpeaks (link) function can calculate FWHM.

See the documentation section on Peak Prominences (link) for an illustration and discussion.

[Veronica Morales]

Thanks again. I will check this out and give it a try. I'm sure I will have more questions.

[Star Strider]

As always, my pleasure.

[Veronica Morales]

Hi again...I'm still trying to solve this. Attached is the plot I get now using 'spline' instead of 'linear' interpolation, althougt they look exactly the same. I am missing something... :(Here's my code:

clear all

clc

% Script to find the dark level of a histogram on a greyscale image

I = imread('dog3.tif');

%figure; imshow(I)

%Calculating and plotting histogram of image

y = imhist(I);

f1 = gradient(y); %f1 = diff(y);

f2 = gradient(f1); %f2 = diff(f1);

figure; plot(y)

%====================================================

% https://www.mathworks.com/matlabcentral/answers/295156-how-to-find-the-inflection-point-of-a-curve-in-matlab

% [d,s,r] = xlsread('cloudy snow 30ppmGE.xlsx');

% I = -d(:,1); % Current

% E = d(:,2); % Potential

% t = E(E<=0);

% y = I(E<=0);

%[b,S,mu] = polyfit(x, y, 6);

% fy = polyval(b,t,S,mu);

% y = fy;

%

d1y = gradient(y); % Numerical Derivative

d2y = gradient(d1y); % Numerical Second Derivative

index = find(d1y == 1); % I'm trying to find where f'=1, but maybe I'm wrong here

x_infl = interp1(d1y, max(d1y), 'spline', 'extrap'); % Why do I get 0 with linear?

y_infl = interp1(y, x_infl, 'spline', 'extrap'); % Find ‘y’ At Maximum Of First Derivative

slope = interp1(d1y, x_infl, 'spline', 'extrap'); % Slope Defined Here As Maximum Of First Derivative

intcpt = y_infl - slope*x_infl; % Calculate Intercept

tngt = slope + intcpt; % Calculate Tangent Line

figure(1)

plot(y)

hold on

%plot(fy)

%plot(d1y, '-.m', d2y, '--c') % Plot Derivatives (Optional)

plot(tngt, '-r', 'LineWidth',1) % Plot Tangent Line

plot(x_infl, y_infl, 'bp') % Plot Maximum Slope

hold off

grid

legend('y(t)', 'y(t) Fit', 'dx/dy', 'd^2x/dy^2', 'Tangent', 'Location','E')

axis([xlim min(min(y),intcpt) ceil(max(y))])

%====================================================

[Star Strider]

I cannot run your code since I do not have the files.

A linear interpolation and a spline interpolation may not look (or be) significantly different if the points are closely spaced and the curve is relatively linear in the region over which you are interpolating it.

I doubt that you are missing anything.

[Veronica Morales]

this is the image I'm using. See attached jpg

[Star Strider]

I am still doing my best to figure out what you are doing in your code with respect to the inflection points.

Meanwhile, one way to find FWHM is to use the midcross (link) function. (It was introduced in R2012a.)

specifically:

mcy = midcross(y);

produces:

mcy =
    23.4983
    38.7014

So:

FWHM = diff(mcy)
FWHM =
    15.2031

Those values are in terms of the ‘x’ variable being defined as:

    x = 1:numel(y);

so essentially the indices of ‘y’.

Answer 2

[Veronica Morales]

Whoa! Thanks for asking. Well my intension is automatization of a process. I need to select that value as the dark level, when the histogram start changing. Don't have more info of why it has to be that point. AS of today we are doing it manually, and want to get away from that, because I have to do the same with hundreds of images...kinda' tedious heh?! I have no clue about all those other methods you mentioned (well getting the FWHM I know). So based on what I was doing, selecting that point visually, I thought it was similar as selecting an inflection point.

[Veronica Morales]

hum...you are right (obviously) maybe what I need is the FWHM...now I need to figure out how to do it.