I had the same problem with surface DEM's. The best solution I found in Matlab was using the scatteredInterpolant class, it is inbuilt in Matlab. It allows Natural neighbour interpolation (that is a class of weighted distance interpolation as suggested in previous comments). The usage is like this:
% IF X and Y are the coords of each grid cell (you can use meshgrid to get
invalidValue=NaN; % or -99999 or whatever is in your case
gapIdxs= DEM==invalidValue; %Indexes of the void cells
interpolant = scatteredInterpolant(X(~gapIdxs),Y(~gapIdxs),DEM(~gapIdxs),'natural');
% interpolant is now a function, and we can evaluate it at any point you
% want, in this case we evaluate it in the void cells
% Now we patch the DEM
That produce a result way better than cubic interpolation in my experience. If that is not enough, before you start doing manual picking of features you can explore Kriging interpolation.
Finally if you what to force directionality, you can rotate your data (using matlab imrotate) so that your main direction is vertical and then interpolate column by column (using matlab interp1).