Sharpening a 2D histogram and finding peaks

Can you help me?

I have a bunch of 2D histograms that I'd like to sharpen so that I can identify some features. Here is an example image:

Each vertical line in the image is a scan (row) in my data and I really want to operate on a per-row basis. The sharpening needs to be primarily in the vertical (value) direction.

I'm trying to identify the center (i.e., median, mode) of the upper and the lower peak for each scan. Any intermediate peaks (for example in scans 3-10) are of lesser interest. Finding the center of the the lower peak is the primary challenge.

I have attached a data file containing the data I used to generate the plot. Here is the plotting code:

function PlotHist(plotData, plotTitle)
    figure;
    surf(plotData); shading("flat");    
    xlim([0, 4096]); xticks(0 : 512 : 4096); xlabel("Laser Value");
    ylim([0, size(plotData, 1)]); ylabel("Scan");
    cb = colorbar; cb.Label.String = '# points';
    view([90, -90])
    title(plotTitle, Interpreter = 'none');
end

I've reviewed and tried a number of the Matlab answers on deblurring and sharpening, but I'm struggling to focus my results in the vertical direction.

Answer 1

[Steven Lord]

You're not plotting a histogram using histogram or histogram2, but did you create the data that you passed into surf using one of those two tools (or the histcounts* equivalents)?

You might want to see if you can use the islocalmax2 and islocalmin2 functions to locate your peaks and valleys and use that information to compute the center locations. The FlatSelection name-value argument may even give you what you're looking for directly from islocalmax2 or islocalmin2.

Alternately, if you want to do effectively 1-dimensional local extrema finding, try using islocalmax and islocalmin with the dim input argument.

[Jim McIntyre]

Yes, I created the histogram using histcounts and I've already used findpeaks with some options to tune the results. I also did some convolutional smoothing on the data in an attempt to refine the peaks. These are all great functions.

The problem is that I have a wide variety of kinds of data. Some of the histograms have two three, or even four well-defined peaks (so that's where tuning findpeaks comes in). Some others have either no upper or lower peak, or a kind of smeared appearance in the middle. I don't want to get into posting all of the various cases, but in this case, which is relatively common, Findpeaks doesn't really have enough to go on at the bottom. So I'm trying to sharpen the peaks to make it work better.

Given the kinds of data I'm working with, local maxes seem to be less useful than findpeaks. I'll take a look at FlatSelection, though.

Answer 2

[dpb]

load HistSmoothed
surf(HistSmooth); shading("flat");
xlim([-inf inf]), ylim(xlim)
xlabel('X'), ylabel('Y'), zlabel('Z')

figure
surf(HistSmooth); shading("flat");
xlabel('X'), ylabel('Y'), xlabel('Z')
HMax=max(HistSmooth(:))

HMax = 89.4286

zlim([HMax/3 inf])
view([90 -90])
xlim([-inf inf]), ylim(xlim)

You can then try

xlim([2000 3800])
box on

or how about

hAx=gca;
hAx.XDir='reverse';

If you try to manipulate the scale in the displayed plane, it's confusing as the dickens...instead, you would need to "zoom" in the z-axis direction. The above displays just the top 2/3rds which expands the scale by 1/3rd over full scale.

However, it's not at all clear to me as to what it is you're wanting to do that makes this the chosen view although it does show that there are some paths across the 2D histogram that do have a significantly fewer number of counts in those channels.

Hmmm....maybe

figure
subplot(2,1,1)
surf(HistSmooth); shading("flat");
xlabel('X'), ylabel('Y'), zlabel('Z')
HMax=max(HistSmooth(:))

HMax = 89.4286

zlim([HMax/3 inf])
view([90 -90])
xlim([2000 3800])
box on
%hAx=gca;
%hAx.XDir='reverse';

subplot(2,1,2)
H1D=sum(HistSmooth,2);        % marginal distribution with X
plot(H1D)
xlabel('Y'), ylabel('Z')
hAx=gca;
hAx.YDir='reverse';

?

[Jim McIntyre]

Thank you. I've looked at that. Functionally it's pretty similar to seeking a local max as discussed in the post above. The actual areas that I'm trying to get information on are these.

What I'm hoping to be able to do is to sharpen all of the peaks, essentially concentrating the energy around the local medians so that findpeaks or islocalmax can do a better job.

I've graphed this stuff six ways from Friday, including both of the options you suggested, and while I can subjectively "see" the peaks there's a bunch of hidden energy that gets in the way of objectively identifying these robustly. I will literally need to do this programmatically for millions of scans, and this is one of the simpler cases.

What I'm really looking for is advice on how to tune a sharpening algorithm to concentrate the vertical distribution of energy around each local peak.

FWIW, I've also looked at edge detection algoriths, and those convince me that I need to "focus" the vertical direction.

[dpb]

load HistSmoothed
surf(HistSmooth); shading("flat");
xlabel('X'), ylabel('Y'), zlabel('Z')
ylim([15 50]); xlim([2500 3000]), zlim([0 60])
clim(zlim)

So it's that kinda' thing you're trying to do something with, I gather?

Answer 3

[Star Strider]

'I'm trying to identify the center (i.e., median, mode) of the upper and the lower peak for each scan.'

I am not absolutely certain what you want. I selected a series of rows that appear to correspond to one of the circled areas, and then did some simpleindexing and statistics. I plotted the data I presented to ghe functions, and then commented it out here. Uncomment those lines to plot the data. I added the standard deviation.

Try something like this --

load('HistSmoothed.mat')
% whos('HistSmoothed.mat')

rowv = (1:5)+20;

xv = 0:size(HistSmooth,1)-1;
yv = 0:size(HistSmooth,2)-1;
[Xm,Ym] = ndgrid(xv,yv);

figure
plot3(Xm(rowv,:).', Ym(rowv,:).', HistSmooth(rowv,:).')
grid
axis('padded')

idxrngv = -200:200;
for k1 = rowv %1:size(HistSmooth,1)
    k1
    [pks,locs] = findpeaks(HistSmooth(k1,:), MinPeakProminence=1.5)
    for k2 = 1:numel(locs)
        idxrng = idxrngv + locs(k2);
        muest(k1,k2) = mean(HistSmooth(k1,idxrng));
        stdest(k1,k2) = std(HistSmooth(k1,idxrng));
        modest(k1,k2) = mode(HistSmooth(k1,idxrng));
        
        % figure
        % plot(idxrng, HistSmooth(k1,idxrng))
        % % hold on
        % % plot(HistSmooth(k1,idxrng))
        % % hold off
        % grid
        % xlabel("Row Index Range")
        % ylabel("Histogram Value")
        % title("Column = "+k1+", Peak = "+k2)
    end
end

k1 = 21

pks = 1×2

3.6667 89.4286

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2574 3320

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

k1 = 22

pks = 1×2

3.9524 87.7143

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2600 3343

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k1 = 23

pks = 1×2

3.4286 86.1905

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2631 3366

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k1 = 24

pks = 1×2

3.4762 86.1429

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2664 3386

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k1 = 25

pks = 1×2

3.3810 85.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2690 3407

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

format shortG
disp(muest)

            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
      0.91414       8.9736
      0.87329       8.9538
      0.89514       9.0049
      0.89835        9.046
      0.88089       9.0064

disp(stdest)

            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
       1.0375       20.855
        1.082       20.655
       1.0434       20.672
       1.0189       20.596
       1.0427        20.28

disp(modest)

     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0

format short

.

EDIT -- (27 Nov 2025 at 16:54)

I'm not certain what you intend by 'smoothing'. One option is to filter the individual peaks, as this illustrates. Choose the cutoff frequency ('Fco') for the lowpass filter to give the result you want.

Alternatively, try something like this --

rowv = (1:5)+20;
idxrngv = -200:200;
Fco = 1/50;
Fs = 1/mean(diff(idxrngv));
for k1 = rowv %1:size(HistSmooth,1)
    k1
    [pks,locs] = findpeaks(HistSmooth(k1,:), MinPeakProminence=1.5)
    for k2 = 1:numel(locs)
        idxrng = idxrngv + locs(k2);
        HistSmoothFilt(k1,:) = lowpass(HistSmooth(k1,:), Fco, Fs);
        muest(k1,k2) = mean(HistSmoothFilt(k1,idxrng));
        stdest(k1,k2) = std(HistSmoothFilt(k1,idxrng));
        modest(k1,k2) = mode(HistSmoothFilt(k1,idxrng));
        
        % figure
        % plot(idxrng, HistSmoothFilt(k1,idxrng))
        % % hold on
        % % plot(HistSmooth(k1,idxrng))
        % % hold off
        % grid
        % xlabel("Row Index Range")
        % ylabel("Filtered Histogram Value")
        % title("Column = "+k1+", Peak = "+k2)
    end
end

k1 = 21

pks = 1×2

3.6667 89.4286

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2574 3320

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k1 = 22

pks = 1×2

3.9524 87.7143

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2600 3343

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k1 = 23

pks = 1×2

3.4286 86.1905

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2631 3366

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k1 = 24

pks = 1×2

3.4762 86.1429

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2664 3386

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k1 = 25

pks = 1×2

3.3810 85.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

2690 3407

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

format shortG
disp(muest)

            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
      0.91406       8.9738
      0.87326       8.9541
      0.89511       9.0051
       0.8982       9.0461
        0.881       9.0064

disp(stdest)

            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
            0            0
       1.0353       20.844
       1.0796       20.645
       1.0415       20.662
       1.0172       20.585
       1.0412        20.27

disp(modest)

     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0

format short

.

[dpb]

I was also confused by the subject of wanting to "sharpen" peaks wherein the circled areas seem to be areas of smaller magnitudes that globally don't look much like single peaks. The only idea I had would be to decimate by a factor in those areas combining counts and then recasting the total into the center of (say) 3x3 or 5x5 areas with the outer cells then being zeros. That would certainly create peaks where there aren't specifically identified ones now. Does that constitute "sharpening"? blockproc might be useful...

Answer 4

[Jim McIntyre]

Thank you for all of the comments. I've been away from the office for a few days, but I will definitely be looking at these.