Memory saving method to interpolate a large scattered dataset

2016 12 8
1 个回答
回答已采纳
更新时间：2019 6 2
10 次查看（30 天）

0 个投票

Hello,
I have a quite large dataset of about 57 million uniformly gridded density samples in 3D space (four column vectors x, y, z and d of length 5.7e7). For a reference frame transformation, I have to apply a rotation and translation on x, y and z which kind of destroys the gridded structure. In order to fix this, I wanted to run an interpolation on a new uniform grid. I attached an illustration with a 2D simplification of the problem. The original gridded dataset given in black is first transformed and then a new gridded dataset given in red should be sampled from the transformed original dataset.
Because the transformed dataset does not comply with the meshgrid format, I cannot use interp3 and have go for a scattered data interpolation approach like griddata or scatteredInterpolant. The problem is that my memory (16 GB) as well as my swap (another 16 GB) fill up quite quickly. So I think the approach to estimate an interpolant for all data points is not the best way to go here.
Does anyone have a better approach for this problem?

请先登录，再回答此问题。

 采纳的回答

Jot We
Jot We 2016-12-13

1 个投票

Just to close this question: I finally found a solution by dividing my dataset into slices and using inverse distance weighting in combination with a k-d tree to generate an interpolation. The computation time and memory usage are quite reasonable.

2 个评论
显示 无 隐藏 无

puni
puni 2019-5-31
Jot, can you please provide an example?
Jot We
Jot We 2019-6-2
Hi puni,
This is a quite old question, but I found the following code snippet in my archive. Maybe it helps.
% Set parameters
powerParameter = 4;
% Initialize coordinates and densities
[globalCoordinatesX, globalCoordinatesY, globalCoordinatesZ] = meshgrid( ...
0:cubeLength:((dimensions(1) - 1) * cubeLength), ...
0:cubeLength:((dimensions(2) - 1) * cubeLength), ...
0:-cubeLength:-((dimensions(3) - 1) * cubeLength) ...
);
globalCoordinates = [globalCoordinatesX(:), globalCoordinatesY(:), globalCoordinatesZ(:)];
clear globalCoordinatesX globalCoordinatesY globalCoordinatesZ;
globalCoordinates = sortrows(globalCoordinates, [-3, 2, 1]);
globalDensityMap = densityMap(:);
clear densityMap;
% Transform and interpolate density map for the thigh
thigh.length = norm(joints.HJ_R - joints.KJ_R);
thigh.landmarks = thighContour.landmarks;
thigh.joints = thighContour.joints;
thigh.hipPlane = thighContour.hipPlane;
thigh.kneePlane = thighContour.kneePlane;
xAxis = cross(landmarks.LFC_R - joints.HJ_R, landmarks.MFC_R - joints.HJ_R);
xAxis = xAxis / norm(xAxis);
yAxis = joints.HJ_R - joints.KJ_R;
yAxis = yAxis / norm(yAxis);
zAxis = cross(xAxis, yAxis);
zAxis = zAxis / norm(zAxis);
translationVector = joints.HJ_R;
rotationMatrix = [xAxis, yAxis, zAxis]';
transformedCoordinates = rotationMatrix * bsxfun(@minus, globalCoordinates', translationVector);
sortedData = sortrows([transformedCoordinates', globalDensityMap], [2, 1, 3, 4]);
transformedCoordinates = sortedData(:, 1:3);
transformedDensityMap = sortedData(:, 4);
clear sortedData;
thigh.limitsX = [floor(min(min(thighContour.contour(:, :, 1))) / cubeLength) * cubeLength, ceil(max(max(thighContour.contour(:, :, 1))) / cubeLength) * cubeLength];
thigh.limitsY = [floor(min(min(thighContour.contour(:, :, 2))) / cubeLength) * cubeLength, ceil(max(max(thighContour.contour(:, :, 2))) / cubeLength) * cubeLength];
thigh.limitsZ = [floor(min(min(thighContour.contour(:, :, 3))) / cubeLength) * cubeLength, ceil(max(max(thighContour.contour(:, :, 3))) / cubeLength) * cubeLength];
[thighCoordinatesX, thighCoordinatesY, thighCoordinatesZ] = meshgrid( ...
thigh.limitsX(1):cubeLength:thigh.limitsX(2), ...
thigh.limitsY(2):-cubeLength:thigh.limitsY(1), ...
thigh.limitsZ(1):cubeLength:thigh.limitsZ(2) ...
);
thighCoordinates = [thighCoordinatesX(:), thighCoordinatesY(:), thighCoordinatesZ(:)];
clear thighCoordinatesX thighCoordinatesY thighCoordinatesZ;
thighCoordinates = sortrows(thighCoordinates, [2, 1, 3]);
thighDensityMap = zeros(size(thighCoordinates, 1), 1);
valuesY = thigh.limitsY(2):-cubeLength:thigh.limitsY(1);
pool = parpool('local');
for valueIndex = 1:length(valuesY);
% Get value range on y-axis
startValueY = valuesY(valueIndex);
endValueY = startValueY + cubeLength;
% Find start and end indeces for sorted transformed coordinates in an
% extended value range
startIndex1 = find(transformedCoordinates(:, 2) >= (startValueY - (sqrt(3) * cubeLength)), 1, 'first');
endIndex1 = find(transformedCoordinates(:, 2) < (endValueY + (sqrt(3) * cubeLength)), 1, 'last');
% Create search tree for coordinates within the extended value range
kdTree = KDTreeSearcher(transformedCoordinates(startIndex1:endIndex1, :));
% Find start and end indeces for sorted thigh coordinates
startIndex2 = find(thighCoordinates(:, 2) >= startValueY, 1, 'first');
endIndex2 = find(thighCoordinates(:, 2) < endValueY, 1, 'last');
% Run interpolation based on inverse distance weighting with a given
% power parameter
parfor coordinateIndex = startIndex2:endIndex2
[indices, distances] = knnsearch(kdTree, thighCoordinates(coordinateIndex, :), 'K', 8);
if any(distances > sqrt(3) * cubeLength)
thighDensityMap(coordinateIndex) = 0;
elseif any(distances == 0)
thighDensityMap(coordinateIndex) = transformedDensityMap(startIndex1 + indices(distances == 0) - 1);
else
weights = 1 ./ distances.^powerParameter;
thighDensityMap(coordinateIndex) = sum(weights .* transformedDensityMap(startIndex1 + indices - 1)') / sum(weights);
end
end
% Print status
fprintf('STATUS: %.1f %% of thigh interpolation done.\n', 100 * valueIndex / length(valuesY));
end
delete(pool);
% Reshape density map for the thigh
fprintf('STATUS: Reshaping thigh density map.\n');
thigh.densityMap = zeros(thigh.limitsX(2) - thigh.limitsX(1) + 1, thigh.limitsY(2) - thigh.limitsY(1) + 1, thigh.limitsZ(2) - thigh.limitsZ(1) + 1);
valuesX = thigh.limitsX(1):cubeLength:thigh.limitsX(2);
valuesY = thigh.limitsY(2):-cubeLength:thigh.limitsY(1);
for indexY = 1:size(thigh.densityMap, 2)
indicesY = (thighCoordinates(:, 2) == valuesY(indexY));
for indexX = 1:size(thigh.densityMap, 1)
indicesYX = indicesY & (thighCoordinates(:, 1) == valuesX(indexX));
thigh.densityMap(indexX, indexY, :) = thighDensityMap(indicesYX);
end
end
clear transformedCoordinates transformedDensityMap thighCoordinates thighDensityMap thighContour indicesY indicesYX valuesX valuesY;

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Interpolation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by