interpolate a 2D map (array) of density samples and preserve the total ?

Good day;

I have a 2D map of data (think: signal strenght per unit sold angle). I want to resample to a finer grid (say from 1x1 deg to 0.1x0.1 deg) while preserving the total strenght.

I tried with interp2, but it does not work. as an example:

datamap = rand(5,7)*10  ;
Xcoarse = [0:6]  ;
Ycoarse = [0:4] ;
[X, Y] = meshgrid (Xcoarse, Ycoarse)  ;
sampleratio = 10  ;
Xfine = linspace (min(Xcoarse), max(Xcoarse), numel(Xcoarse)*sampleratio-(sampleratio-1))  ;
Yfine = linspace (min(Ycoarse), max(Ycoarse), numel(Ycoarse)*sampleratio-(sampleratio-1))  ;
[finemeshX, finemeshY] = meshgrid(Xfine, Yfine)   ;
datamapfine = interp2 (X, Y, datamap, finemeshX, finemeshY, 'linear')  ;

When I integrate (i.e out = sum(sum(data)) * dx * dy) 'datamap' and 'datamapfine' (accounting for the different surface units), 'datamapfine' has picked up another ~50% of signal strenght. The same effect persists for any type of interpolation: "cubic', 'spline', etc...

For a series of reason too long to explain here, taking the integral of the two maps and re-normalization of the resampled map after-the-fact is not an acceptable solution (it would be too easy if it were ! ) to the problem.

Can anyone offer a solution? thanks.

Answer 1

[arushi]

Hi Giovanni,

I understand you want to preserve the total signal strength directly during the resampling process without relying on post-resampling normalization. To achieve this, you can use a different approach that directly adjusts the interpolation weights to ensure the preservation of the total signal strength.

One approach to accomplish this is by modifying the interpolation weights based on the ratio of the areas covered by the original and resampled grid cells. This can be done by scaling the interpolation weights to account for the change in area.

Here's an example of how you can achieve this in MATLAB:

% Given data 

datamap = rand(5, 7) * 10; 
Xcoarse = 0:6; 
Ycoarse = 0:4; 
sampleratio = 10; 

% Create finer grid 
[X, Y] = meshgrid(Xcoarse, Ycoarse); 
[Xq, Yq] = meshgrid(linspace(min(Xcoarse), max(Xcoarse), numel(Xcoarse) * sampleratio), ... 

                    linspace(min(Ycoarse), max(Ycoarse), numel(Ycoarse) * sampleratio)); 

 
% Calculate the areas of the original and resampled grid cells 
coarse_area = (Xcoarse(2) - Xcoarse(1)) * (Ycoarse(2) - Ycoarse(1)); 
fine_area = (Xq(1,2) - Xq(1,1)) * (Yq(2,1) - Yq(1,1)); 

% Calculate the interpolation weights adjustment factor 
weight_adjustment = fine_area / coarse_area; 

% Perform the interpolation with adjusted weights 
datamapfine = interp2(X, Y, datamap, Xq, Yq, 'linear') * weight_adjustment; 

In this approach, the interpolation weights are adjusted by a factor that accounts for the change in area when resampling to a finer grid and the signal strength is reduced to around 30%. By directly modifying the interpolation weights, this method aims to preserve the total signal strength during the resampling process without the need for post-resampling normalization.

Please note that this method assumes a regular grid and may need further adjustments if the grid is irregular or if additional considerations are required for the specific nature of your data.

Hope this helps.

[Giovanni de amici]

good day Arushi;

thanks for the suggestion.

unfortunately your solution is equivalent at conflating the dx and dy integrands into one single term (dx*dy, or weight_adjustment in your code) and that does nothing to solve my problem. I have become convinced that while both interp1 and interp2 do a decent job at resampling a smooth-changing function (e.g. a sinusoidal), interp2 is unable to handle a sequence of numbers with a non-continuous derivative.

as an example I post a code snippet that contains two distributions. while the resampling returns in both case an integral that is within 3% of the original, the shape of the resampled random-number distribution is nothing like its original one

datamap = zeros( 5, 7); 
Ycoarse = linspace(0.5, 4.5, 5)   ;
Xcoarse = linspace (0.5, 6.5, 7)  ;
Yphase = Xcoarse/5.5 * pi  ;
Xphase = Ycoarse/7.5 * pi  ;  
for i=1:5
    for j = 1:7
        datamap(i,j) = sin(Yphase(j))*sin(Xphase(i));
    end
end
datamap = rand( 5, 7) * 10;  % comment this out for case 1
sampleratio = 10; 
% Create finer grid 
[X, Y] = meshgrid(Xcoarse, Ycoarse); 
Yfine = linspace (0, 4.9, 5*sampleratio)  ;
Xfine = linspace (0, 6.9, 7*sampleratio)  ;
[Xq, Yq] = meshgrid( Xfine,  Yfine )   ;
% Calculate the areas of the original and resampled grid cells 
coarse_area = (Xcoarse(2) - Xcoarse(1)) * (Ycoarse(2) - Ycoarse(1)); 
fine_area = (Xq(1,2) - Xq(1,1)) * (Yq(2,1) - Yq(1,1)); 
% Calculate the interpolation weights adjustment factor 
weight_adjustment = fine_area / coarse_area; 
% Perform the interpolation with adjusted weights 
datamapfine = interp2(X, Y, datamap, Xq, Yq, 'spline') * weight_adjustment;

[Bruno Luong]

@Giovanni de amici

Your code does not only do interpolation but extrapolation beyond the original sample point; by spline function.

MATLAB spline (with not a knot bc) is quite unstable and you should not any consistent result with that extrapolation.

If you use sum as integration, the underlined function is constant by pixel centedred at your data grid. You should use nearest interpolation to be coherent with your integrationscheme.

You use the moste unstable extrapolation NOT compatible with your integration scheme.

So either you use nearest extrapolation, or with my answer do not use extrapolation but with linear interpolation.

You want both smooth (C^2) interpolation/extrapolation and at the same time preserve integration. To me you want too much and both requirements can't have a solution.

Answer 2

[Bruno Luong]

Use trapz for the integration (not sum), it should match (up to round off) with linear interpolation and interger subsamping

format long
m = 5;
n = 7;
dx = 1;
dy = 1;
Xcoarse = (0.5:n-0.5)*dx;
Ycoarse = (0.5:m-0.5)*dy;;
data = rand(m,n);

Icoarse = trapz(trapz(data,2),1)*(dx*dy)

Icoarse = 
  11.238805786450900

nx = 8;
ny = 8;

Xfine = linspace(min(Xcoarse),max(Xcoarse),(length(Xcoarse)-1)*nx+1);
Yfine = linspace(min(Ycoarse),max(Ycoarse),(length(Ycoarse)-1)*ny+1);
datafine = interp2(Xcoarse, Ycoarse', data, Xfine, Yfine', 'linear');

dxfine = mean(diff(Xfine)); % dx / nx;
dyfine = mean(diff(Yfine)); % dy / ny;
Ifine = trapz(trapz(datafine,2),1)*(dxfine*dyfine)

Ifine = 
  11.238805786450902