Heatmap and surf plot of same data is looking different. Also is there a better way to 2D smooth

Hello, I am trying to estimate the background of this image. One way I would like to try is to break the image up into small regions and take the minium value in each region and then peform a surface fit to these. The aim would be to then subtract this from my main image (one thing I definetly need is smoothness, i.e no sharp jumps between regions)

My raw image showing the regions

The min in each region is obtained via blockproc (data attached)

blockMin =
183.00        185.00        184.00        195.00        198.00        185.00        203.00
194.00        204.00        188.00        192.00        217.00        190.00        212.00
202.00        213.00        204.00        191.00        208.00        210.00        209.00
204.00        205.00        225.00        204.00        202.00        210.00        211.00
225.00        240.00        219.00        209.00        210.00        227.00        210.00
225.00        213.00        217.00        241.00        223.00        226.00        223.00
219.00        222.00        232.00        233.00        217.00        238.00        231.00

Very strangely, the heatmap and surf plot of the same data set is showing a difference (I know I must rotate the surf plot, but its clear the peak regions are in different positions (shown by green circles). On the heatmap, the value of 241 is on the 3rd row, yet on the surf plot its on the 2nd row and on the wrong column.

This is my code (Im showing the pure surf plot above, not the smoothdata2 version)

            blockFwhm;
            figure   
            h=heatmap(blockFwhm, 'Colormap',cool,'FontSize',10);    % app.PanelScan,
            h.CellLabelFormat = '%.0f';

            figure
            ax1=subplot(1,2,1); ax2=subplot(1,2,2);
            surf(ax1,blockFwhm); colormap cool; colorbar
            
            % Try to smotth the data  *** Although this isn't very good  so
            % far
            B=smoothdata2(blockFwhm);
            s=surf(ax2,B); colormap cool;

[Jason]

I wonder if this is a better approach?

1: find the spots of as many as posible (x1,y1)

2: Dilate and put these values to zero

    binary = false(size(image));
    binary(sub2ind(size(image), round(y1), round(x1))) = true;
    
    %Dilate slightly
    SE=8;
    se = strel('disk',SE);
    mask=imdilate(binary,se,'same'); 
    figure
    ax1=subplot(1,3,1);  ax2=subplot(1,3,2);
    myImagesc(app,ax1,mask);    
 

3: Apply Mask To Image

    % Apply Mask to image
    I=image;
    I(mask) = 0;  
    myImagesc(app,ax2,I); 
    

But Im stuck now how to smmoth this image but ignore the zero values.

Answer 1

[Walter Roberson]

Notice that the heatmap shows the center of the blocks labeled 1 through 7, whereas the surf() shows 1 and 7 at the extreme edges of the plot.

heatmap() presents each value as-is in a different cell location.

surf() on the other hand treats the input values as vertices and interpolates face colors from the vertices. The number of faces is one less than the corresponding dimension.

You should be using imagesc() instead of surf()

[Jason]

Thankyou Walter for pointing this out, I never knew that.

Ideally, I would still like to 3D plot for visualisation as my ultimate goal is to fit some smooth surface to this and I'd prefer to see the goodness of the fit by 3D

[Walter Roberson]

The shortcut for viewing an array as a 3D surface, is to use warp

[Jason]

I couldnt get warp to work

[Jason]

I've nearly got scatter3 working - just can't get the colour mapping right. I want the colour to reflect the "Z" coordinate

  B=smoothdata2(blockFwhm);
  D=imgaussfilt(B,2);
  imagesc(ax2,D); colormap cool;
    
    
  % Try using Scatter3 to visualise data         
  [x,y] = meshgrid(1:ngrids,1:ngrids);   % ngrids is mumber of data points in x and y

  figure
  s=scatter3(x,y,D,'filled','MarkerEdgeColor','k');
  box on
  C = jet(numel(x));
  s.CData = C
  set(gca, 'colormap', C)
  colorbar()

But this gives:

[Walter Roberson]

pointsize = 30;
s=scatter3(x,y,D,pointsize,D,'filled','MarkerEdgeColor','k');

without setting s.CData

[Walter Roberson]

x = linspace(0, 1, 25);
y = linspace(2, 3, 25).';
z = x.^2 - cos(y);
warp(x, y, z, z)
colormap jet

●

[Jason]

Hi Walter, thanks for your suggestions. After a bit of searching I realised that x and y need to be 1D vectors. So I got this to work (realising my D data is actually an array)

D =
203.4923  204.3522  204.8332  205.0403  205.3193  205.8575  206.4547
206.1093  206.8875  207.3046  207.4248  207.5679  207.9480  208.4290
209.5499  210.2312  210.5939  210.6524  210.6793  210.9023  211.2564
213.3235  213.9204  214.2688  214.3278  214.3003  214.4067  214.6498
216.7701  217.3224  217.7213  217.8667  217.8735  217.9323  218.1024
219.3818  219.9351  220.4423  220.7446  220.8620  220.9405  221.0813
220.9900  221.5733  222.2104  222.6935  222.9567  223.0975  223.2430

ngrids=11;
[x,y] = meshgrid(1:ngrids,1:ngrids);
xy = [repelem(x(:), 1, 1) repmat(y(:), 1, 1)];    % This is the crucial step!
x=xy(:,1);                                        % Must use column vectors
y=xy(:,2);                                        % Must use column vectors

figure
s=scatter3(x,y,D(:),30,D(:),'filled','MarkerEdgeColor','k');
colormap turbo
grid on

I also wanted to try to fit and realised that fit and plot can do what I want

sf = fit([x, y],D(:),'poly33')
plot(sf,[x,y],D(:))
colormap turbo
grid on

So my final ask is to how do I find the fitted surface Z values at themarker positons so I can then subtract form my image. Remember my values D that I have "surface" fitted above are representative of the "regions"in my main image i.e. the cyan spots below.

Im guessing I need to somehow "stretch" the fit to math my image dimensions

I have tried this

[x,y] = ndgrid(1:11)
F = griddedInterpolant(x,y,D);
% Use a finer mesh to query the interpolant and improve the resolution.

% Assume that there are 10x10 pixels in each of my regions above. So need
% to make ngrid steps 1/10th I think
[xq,yq] = ndgrid(1:0.1:11);     
vq = F(xq,yq);
figure
surf(xq,yq,vq)

Which gives this:

But Im not sure this can be directly subtracted from my image, as there are 7 regions in my image and each region has 10 pixels ,in each dimension so I would expected this to extend to 70 in the x and y, allowing me to directly subtract this from my raw image and hence background correct it

[Jason]

I managed to get what I wanted - to "stretch" so to match to my raw image I just used imresize and used the number of pixels per ROI region as the parameter in the imresize.

This is my result - Im really pleased!

You can clearly see the min 5% intensity pixels in the raw image are towards the top left corner, but once the background correction is done its more uniform over the whole image

Shame I couldn't get the mask approach to work.