Make calculation on line profile (Intensity) through object on image

Hello, I have an image that I want to approx the size of the central blob.

For this I just create a line profile fo the intensity and just calculate the width corresponding to e.g. 20% of the max-min value (i.e." Y20")

(Data attached)

My function to calc this is below. However it doesn't always perform the caluculation on the correct peak, as in this case. It should be on the more central peak, so the calc I want is the length of the green line and not the blue line.

       function [xx,xx2,dx] = calcWidth2(app,scol, AX,frac,mx,mn) 

               % --------Calcs the 13.5% width (1/e^2), but allows for any 
               %Outputs x coordinate of 13.5 (or what ever fract) is used
               %for the width (dx)
              
               % scol=data   (vector)
               % AX = axis for plots
               % frac=  i.e. 0.135 for 13.5% width (1/e^2)            
                     
               % mx = max(max(scol(:)));
               % mn = max(min(scol(:)));
               
               cross=find(scol > val ) % Find values above the Y20 line
               
               % ***********************************************************************
               % Perhaps here some logic is required that if there is more
               % than 1 peak, use the central one?
               %************************************************************************
               
               %Use interp to get better (finer) resolution
               n=numel(cross);  %get number of point above, then multiply by 20 for finer res
               xdataFine=(linspace(cross(1)-1,cross(end)+1,20*n))';  
               xq=xdataFine; x=1:length(scol); v=scol;
               vq2 = interp1(x,v,xq,'linear');  %spline     %rloess   sgolay
               cross=find(vq2 > val );             
               xx=xdataFine(cross(1));
               xx2=xdataFine(cross(end));  
               dx=xx2-xx;
                       
               plot(AX,xdataFine,vq2,'k:','LineWidth',1); 
               
               xline(AX,xx,'Color','r','LineStyle','--');
               xline(AX,xx2,'Color','r','LineStyle','--');   
               %text(AX,0.02,0.94,['X',num2str(frac*100),'% = ',num2str(dx,'%.1f'),' pxls'],'Units','normalized','FontSize',11,'Color','r','HorizontalAlignment','left');  
               hold(AX,'off');
       end
       

I did try this, but it doesn't work if the unwanted peak is at the beginning

       cross=find(scol > val ) 
       % Account for multiple peaks   
       D=diff(cross);
       idx=find(D>5);
       %size(cross)
       cross(idx:end)=[];

Thanks for any help

Answer 1

[Mathieu NOE]

hello Jason

hope you are doing well

see my suggestion below

hope it helps !

% load data
y = readmatrix('LineProfile');
n=numel(y);
x=(0:n-1);

% define threshold line
threshold = 0.2*max(y);

[max_val,locs] = findpeaks(y,'MinPeakHeight',threshold); % find the peak(s) and x axis location(s)

% if multiple peaks are selected , select the one closest to mid x position
if numel(locs)>1
    dist = abs(locs-round(n/2));
    [m,ind] = min(dist);
    locs = locs(ind);
    max_val = max_val(ind);
end

% find "threshold" crossing values and indexes (NB linear interpolation on x values for better accuracy)
[x0_pos1,s0_pos1,x0_neg1,s0_neg1]= crossing_V7(y,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

% if multiple values of x0_pos1 are returned , select the one closest to mid x position
if numel(x0_pos1)>1
    dist = abs(x0_pos1-round(n/2));
    [m,ind] = min(dist);
    x0_pos1 = x0_pos1(ind);
    s0_pos1 = s0_pos1(ind);
end

% if multiple values of x0_neg1 are returned , select the one closest to mid x position
if numel(x0_neg1)>1
    dist = abs(x0_neg1-round(n/2));
    [m,ind] = min(dist);
    x0_neg1 = x0_neg1(ind);
    s0_neg1 = s0_neg1(ind);
end


figure(1)
plot(x,y,'b',x,threshold*ones(size(x)),'--r',x(locs),max_val,'dr',x0_pos1,s0_pos1,'db',x0_neg1,s0_neg1,'dg','linewidth',2,'markersize',12);grid on
legend('signal','threshold','signal peak value','signal positive slope crossing points','signal negative slope crossing points'  );

% x difference between crossing points
peak_width = - x0_pos1 + x0_neg1;

% add text to the plot
for ck = 1:numel(peak_width)
    x_text = x0_pos1(ck);
    y_text = 1.05*max_val(ck);
    text(x_text,y_text,['Peak Value : ' num2str(max_val(ck))] ); 
    y_text = 1.1*max_val(ck);
    text(x_text,y_text,['Peak Width : ' num2str(peak_width(ck))] ); 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%


% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if ~isempty(ind)

    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)

                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end

    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);

    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end

[Mathieu NOE]

Just to test the code in another case with threshold set at 10% , the line crosses multiple times the data , but the results are correct

% define threshold line
threshold = 0.1*max(y);

results :

[Jason]

Thankyou, but i dont have findpeaks since no signal processing toolbox

[Mathieu NOE]

ok, no problemo

use this code instead

% [max_val,locs] = findpeaks(y,'MinPeakHeight',threshold); % find the peak and its x axis location
[locs,max_val] = peakseek(y,10,threshold); % find the peak and its x axis location

the function peakseek is attached (or if you prefer you can find it at the FEX page : PeakSeek - File Exchange - MATLAB Central )

[Jason]

Thankyou Mathieu. I've managed to utilise your fantastic code, however, when my data is like this, its not picking up the correct green marker positon

(I'm ust consider the red curve for now, but have also just realised there will be a problem with the data in the green curve as on the left there is no crossing point - perhaps in this scenario, its just put as x=0)

[Jason]

Perhaps this addiotional line is required:

if numel(x0_pos1)>1
    x0_pos1(x0_pos1<round(n/2))=[]     %  THIS Additional Line?
    dist = abs(x0_pos1-round(n/2));
    [m,ind] = min(dist);
    x0_pos1 = x0_pos1(ind);
    s0_pos1 = s0_pos1(ind);
end

[Jason]

Its also struggling to find the postive crossing point with this dataset

[Jason]

Actually, once I've got the main peak x location (x(locs)), this seems to work well

   %Upper Crossing point
   cross=find(c1 < threshold ); 
   xx=x(cross)';
   xx(xx<round(x(locs)))=[];
   Upper=xx(1);
   xline(ax5,Upper);

   %Lower Crossing Point
   cross=find(c1 < threshold ); 
   xx=x(cross)';
   xx(xx>round(x(locs)))=[];
   Lower=xx(end);
   xline(ax5,Lower);

[Image Analyst]

Getting a line profile, then doing signal analysis on it to find the bottoms of the peaks seems like a lot more work than simply using regionprops ad I suggested. Plus the line profile method only gives the "size" along the sampling direction, so it would not be as accurate as something that looked at the whole blob, like regionprops.

[Mathieu NOE]

hello again @Jason

try this modified code :

% load data
y = readmatrix('LineProfile3.txt');
n=numel(y);
x=(0:n-1);

% define threshold line
threshold = 0.2*max(y);

% [max_val,locs] = findpeaks(y,'MinPeakHeight',threshold); % find the peak and its x axis location
[locs,max_val] = peakseek(y,10,threshold); % find the peak and its x axis location

% if multiple peaks are selected , use the one closest to mid x range
if numel(locs)>1
    dist = abs(locs-round(n/2));
    [m,ind] = min(dist);
    locs = locs(ind);
    max_val = max_val(ind);
end

% find "threshold" crossing values and indexes (NB linear interpolation on x values for better accuracy)
[x0_pos1,s0_pos1,x0_neg1,s0_neg1]= crossing_V7(y,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

% if multiple values of x0_pos1 are returned , use the one closest to mid x range
% but also eliminate values of x0_pos1 that are after the peak (non valid)
dist = (x0_pos1-round(n/2));
ii = find(dist<0);
[~,ind] = max(dist(ii));
x0_pos1 = x0_pos1(ii(ind));
s0_pos1 = s0_pos1(ii(ind));

% if multiple values of x0_neg1 are returned , use the one closest to mid x range
% but also eliminate values of x0_neg1 that are before the peak (non valid)
dist = (x0_neg1-round(n/2));
ii = find(dist>0);
[~,ind] = min(dist(ii));
x0_neg1 = x0_neg1(ii(ind));
s0_neg1 = s0_neg1(ii(ind));

figure(1)
plot(x,y,'b',x,threshold*ones(size(x)),'--r',x(locs),max_val,'dr',x0_pos1,s0_pos1,'db',x0_neg1,s0_neg1,'dg','linewidth',2,'markersize',12);grid on
legend('signal','threshold','signal peak value','signal positive slope crossing points','signal negative slope crossing points' ,'FontSize',14 );

% x difference between crossing points
peak_width = - x0_pos1 + x0_neg1;

% add text to the plot
for ck = 1:numel(peak_width)
    x_text = x0_pos1(ck);
    y_text = 1.05*max_val(ck);
    text(x_text,y_text,['Peak Value : ' num2str(max_val(ck))] ); 
    y_text = 1.1*max_val(ck);
    text(x_text,y_text,['Peak Width : ' num2str(peak_width(ck))] ); 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%


% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if ~isempty(ind)

    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)

                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end

    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);

    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end

[Jason]

Thankyou

[Mathieu NOE]

as always, my pleasure !

Answer 2

[Image Analyst]

To get the size (area), equivalent circular diameter, mean, max, and min intensities of the blob, why not simply threshold and call regionprops. Something like

mask = grayImage > someThreshold;
props = regionprops(mask, 'Area', 'Perimeter', 'MeanIntensity', 'MaxIntensity', 'MinIntensity', 'EquivDiameter');

See my Image Segmentation Tutorial in my File Exchange:

Image Analyst's File Exchange

It's a generic, general purpose demo of how to threshold an image to find blobs, and then measure things about the blobs, and extract certain blobs based on their areas or diameters.