Memory saving method to interpolate a large scattered dataset

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?

Answer 1

[Jot We]

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.

[puni]

Jot, can you please provide an example?

[Jot We]

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;