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

3 views (last 30 days)
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.

Answers (2)

arushi
arushi on 21 Nov 2023
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.
  2 Comments
Giovanni de amici
Giovanni de amici on 25 Nov 2023
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
Bruno Luong on 25 Nov 2023
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.

Sign in to comment.


Bruno Luong
Bruno Luong on 25 Nov 2023
Edited: Bruno Luong on 25 Nov 2023
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

Categories

Find more on Interpolation in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!