I don't think that what's requested is necessarily the only way to get around the problem, but consider the following image:
%% fit the object to an orthogonal grid
% an object in an image
A = imread('businesscard.jpg');
imshow(A)
You could choose to just correct the image for the perspective distortion. At that point, you don't need to deal with foreshortening anymore. You can easily overlay a regular grid or subdivide the object.
% these are the coordinates of the box corners
% you can get these using getpts() or impixelinfo() or datatips
boxm = [1115 955; 2890 1167; 2898 2221; 821 1918]; % clockwise from NW corner
% assert that this is where they're supposed to be in the output geometry
szo = [384 640]; % business cards are roughly 1:5/3
boxf = [1 1; szo(2) 1; szo(2) szo(1); 1 szo(1)]; % same order as boxm
% transform the image
TF = fitgeotrans(boxm,boxf,'projective');
outview = imref2d(szo);
B = imwarp(A,TF,'fillvalues',255,'outputview',outview);
% now the object can be subdivided easily
imshow(B); hold on
% for example, show a grid on the image
% i'm choosing to scale N to suit the object aspect ratio
% so that the object subdivisions are roughly square
N = [17 10]; % number of points per edge [x y]
x = linspace(1,szo(2),N(1)); % just use output image geometry
y = linspace(1,szo(1),N(2));
[X Y] = meshgrid(x,y);
plot(X,Y,'m.','markersize',20)
Instead, you could try to work in the original image space as requested:
%% fit an orthogonal grid to the object
% an object in an image
A = imread('businesscard.jpg');
% a uniform grid of points
% i'm choosing to scale N to suit the object aspect ratio
% so that the object subdivisions are roughly square
N = [17 10]; % number of points per edge [x y]
x = linspace(0,1,N(1)); % using unit-scale
y = linspace(0,1,N(2));
[X Y] = meshgrid(x,y);
% these are the coordinates of the grid corners
boxm = [0 0; 1 0; 1 1; 0 1]; % same order as boxf
% assert that this is where they're supposed to be in the image geometry
% these are the same points used for boxm in the prior example
boxf = [1115 955; 2890 1167; 2898 2221; 821 1918]; % clockwise from NW corner
% transform the grid
TF = fitgeotrans(boxm,boxf,'projective');
[Xout,Yout] = transformPointsForward(TF,X(:),Y(:));
% show the transformed grid on the image
imshow(A); hold on
plot(Xout,Yout,'m.','markersize',40)
