Meshgrid and interp2 to make a mask for planefitting

For my bachelor assignment I do research on piezoelectric materials (materials that are able to mechanically deform under influence of an electric applied field). When comparing planefits for different regions the actuation of the piezoelecric material in Z-direcction could be obtained. The deformation has been measured in the Z-direction for a default region X=[1:1000] Y = [1:1000], so for each different measurement(I can vary the voltage) I get a .dat file containing nearly 1E6 Z points. The data set is given by the white light interferometer in three columns X,Y,Z. In this particular measurement the x-values run from 0:943 and y-values 0:999, so X,Y and Z have lengths 944000. Now I would like to do a planefitting on a specific region of that 1E6 points, making a meshgrid and interpolation because of missing values(NaN's). The mesghrid on the region [X,Y] = (400:600,50:450). An error returned:

"Error using griddedInterpolant
The grid vectors do not define a grid of points that match the given values."

Furthermore, If one makes a meshgrid, the dataset is scaled from X,Y to Xmesh,Ymesh. so also the region on which one tries do do the planefitting is scaled. How to deal with this? Thanks in advance for your help.

M=load('WLI_m_17oct_PILL002EL3_0V_2.xyz'); %
%ref = load('ref0V_P3_EL2.xyz');
breedtex = M(:,1);
dieptey = M(:,2);
hoogtez = M(:,3);

x = 0:breedtex(end); % definieer de breedte van het meetbereik
y = 0:dieptey(end); % definieer de diepte van het meetbereik

Mvierkant=zeros(1000,1000);                 % maak een lege nieuwe matrix aan
for i = 1 : length(hoogtez)
       Mvierkant(breedtex(i)+1,dieptey(i)+1) = hoogtez(i);
end

sc =0.5                                                       %schaalfactor =sc
[X,Y] = meshgrid(400:600,50:450); 
Mvierkantmesh = Mvierkant(400:600,50:450);
[Xmesh,Ymesh] = meshgrid(400:sc:600,50:sc:450);              
Mres = interp2(X,Y,Mvierkantmesh,Xmesh,Ymesh,'linear');
Mresi = Mres(:);

Ax = zeros(201,401);
Ay = zeros(201,401);
for maski = 1:201                                           % x = [400,600]
    for maskj = 1:401                                       % y = [50,450]
        Ax(maski,maskj) = Xmesh(round((1/sc)*944*(50+maskj)+(1/sc)*maskj)); 
        Ay(maski,maskj) = Ymesh(round((1/sc)*944*(50+maskj)+(1/sc)*maskj));         
        Az(maski,maskj) = Mres(round((1/sc)*944*(50+maskj)+(1/sc)*maskj));   
    end                                                       
end
Ax = Ax(:);
Ay = Ay(:);
Az = Az(:);
sf = fit([Ax, Ay],Az,'poly11');

figure(2); 
surf(Xmesh,Ymesh,Mres)
hold on
surf(Mvierkant)