| 1 |
function [ow]=gcmfaces_quadmap(px,py,ox,oy); |
| 2 |
%[ow]=gcmfaces_quadmap(px,py,x,y); |
| 3 |
%object: compute bilinear interpolation coefficients for x(i,:),y(i,:) |
| 4 |
% in px(i,:),py(i,:) by remapping x(i,:),y(i,:) along with the |
| 5 |
% px(i,:),py(i,:) quadrilateral to the 0-1,0-1 square. |
| 6 |
%inputs: px,py are Mx4 matrices where each line specifies one quad |
| 7 |
%(optional) ox,oy are MxP position matrices |
| 8 |
%outputs: pw are the MxPx4 bilinear interpolation weights |
| 9 |
|
| 10 |
doDisplay=0; |
| 11 |
if nargin==0; |
| 12 |
%the following test case is based upon https://www.particleincell.com/2012/quad-interpolation/ |
| 13 |
%(see end of routine for alternatives) |
| 14 |
px = [-1, 8, 13, -4]; |
| 15 |
py = [-1, 3, 11, 8]; |
| 16 |
ox=0; oy=6; |
| 17 |
doDisplay=1; |
| 18 |
end; |
| 19 |
|
| 20 |
if nargin==2; ox=[]; oy=[]; end; |
| 21 |
|
| 22 |
%solve linear problem for a,b vectors (knowing px,py) |
| 23 |
% logical (l,m) to physical (x,y) mapping is then |
| 24 |
% x=a(1)+a(2)*l+a(3)*m+a(2)*l*m; |
| 25 |
% y=b(1)+b(2)*l+b(3)*m+b(2)*l*m; |
| 26 |
% A=[1 0 0 0;1 1 0 0;1 1 1 1;1 0 1 0]; AI = inv(A); |
| 27 |
% AI=[1 0 0 0;-1 1 0 0;-1 0 0 1; 1 -1 1 -1]; |
| 28 |
% a = AI*px'; |
| 29 |
% b = AI*py'; |
| 30 |
tmp1=px(:,1); |
| 31 |
tmp2=-px(:,1)+px(:,2); |
| 32 |
tmp3=-px(:,1)+px(:,4); |
| 33 |
tmp4=px(:,1)-px(:,2)+px(:,3)-px(:,4); |
| 34 |
a=[tmp1 tmp2 tmp3 tmp4]; |
| 35 |
% |
| 36 |
tmp1=py(:,1); |
| 37 |
tmp2=-py(:,1)+py(:,2); |
| 38 |
tmp3=-py(:,1)+py(:,4); |
| 39 |
tmp4=py(:,1)-py(:,2)+py(:,3)-py(:,4); |
| 40 |
b=[tmp1 tmp2 tmp3 tmp4]; |
| 41 |
|
| 42 |
%chose between the two mapping solutions dep. on sum of interior angles |
| 43 |
[angsum]=gcmfaces_polygonangle(px,py); |
| 44 |
sgn=NaN*px(:,1); |
| 45 |
ii=find(abs(angsum-360)<1e-3); sgn(ii)=1; |
| 46 |
ii=find(abs(angsum+360)<1e-3); sgn(ii)=-1; |
| 47 |
ii=find(isnan(angsum)); |
| 48 |
if length(ii)>0; |
| 49 |
warning('edge point was found'); |
| 50 |
keyboard; |
| 51 |
end; |
| 52 |
|
| 53 |
%solve non-linear problem for pl,pm (knowing px,py,a,b) |
| 54 |
% physical (x,y) to logical (l,m) mapping |
| 55 |
% |
| 56 |
% quadratic equation coeffs, aa*mm^2+bb*m+cc=0 |
| 57 |
if ~isempty(ox); |
| 58 |
x=[px ox]; y=[py oy]; |
| 59 |
else; |
| 60 |
x=px; y=py; |
| 61 |
end; |
| 62 |
a=reshape(a,[size(a,1) 1 size(a,2)]); a=repmat(a,[1 size(x,2) 1]); |
| 63 |
b=reshape(b,[size(b,1) 1 size(b,2)]); b=repmat(b,[1 size(x,2) 1]); |
| 64 |
sgn=repmat(sgn,[1 size(x,2)]); |
| 65 |
% |
| 66 |
aa = a(:,:,4).*b(:,:,3) - a(:,:,3).*b(:,:,4); |
| 67 |
bb = a(:,:,4).*b(:,:,1) -a(:,:,1).*b(:,:,4) + a(:,:,2).*b(:,:,3) ... |
| 68 |
- a(:,:,3).*b(:,:,2) + x.*b(:,:,4) - y.*a(:,:,4); |
| 69 |
cc = a(:,:,2).*b(:,:,1) -a(:,:,1).*b(:,:,2) + x.*b(:,:,2) - y.*a(:,:,2); |
| 70 |
%compute m = (-b+sqrt(b^2-4ac))/(2a) |
| 71 |
det = sqrt(bb.*bb - 4.*aa.*cc); |
| 72 |
pm = (-bb+sgn.*det)./(2.*aa); |
| 73 |
%compute l by substitution in equation system |
| 74 |
pl = (x-a(:,:,1)-a(:,:,3).*pm)./(a(:,:,2)+a(:,:,4).*pm); |
| 75 |
|
| 76 |
ow=[]; |
| 77 |
if ~isempty(ox); |
| 78 |
tmp1=(1-pl(:,5:end)).*(1-pm(:,5:end)); |
| 79 |
tmp2=pl(:,5:end).*(1-pm(:,5:end)); |
| 80 |
tmp3=pl(:,5:end).*pm(:,5:end); |
| 81 |
tmp4=(1-pl(:,5:end)).*pm(:,5:end); |
| 82 |
ow=cat(3,tmp1,tmp2,tmp3,tmp4); |
| 83 |
end; |
| 84 |
|
| 85 |
if doDisplay; |
| 86 |
cols='brgm'; |
| 87 |
% |
| 88 |
figureL; |
| 89 |
%plot original quad and obs |
| 90 |
subplot(1,2,1); |
| 91 |
plot([px px(1)],[py py(1)],'k-','LineWidth',2); hold on; |
| 92 |
for pp=1:4; |
| 93 |
plot(px(pp),py(pp),[cols(pp) '.'],'MarkerSize',64); |
| 94 |
end; |
| 95 |
aa=axis; aa(1)=aa(1)-1; aa(2)=aa(2)+1; aa(3)=aa(3)-1; aa(4)=aa(4)+1; axis(aa); |
| 96 |
grid on; plot(ox,oy,'k.','MarkerSize',64); |
| 97 |
% |
| 98 |
subplot(1,2,2); |
| 99 |
plot([pl(1:4) pl(1)],[pm(1:4) pm(1)],'k-','LineWidth',2); hold on; |
| 100 |
for pp=1:4; plot(pl(pp),pm(pp),[cols(pp) '.'],'MarkerSize',64); end; |
| 101 |
aa=axis; aa(1)=aa(1)-1; aa(2)=aa(2)+1; aa(3)=aa(3)-1; aa(4)=aa(4)+1; axis(aa); |
| 102 |
grid on; plot(pl(5:end),pm(5:end),'k.','MarkerSize',64); |
| 103 |
end; |
| 104 |
|
| 105 |
%swap px and py: |
| 106 |
% ppx=px; ppy=py; |
| 107 |
% px=ppy; py=ppx; |
| 108 |
% |
| 109 |
%shift quad corners: |
| 110 |
% nn=0; |
| 111 |
% px=circshift(px,[0 nn]); |
| 112 |
% py=circshift(py,[0 nn]); |
| 113 |
% |
| 114 |
%flip quad corners (see sum of interior angles): |
| 115 |
% ppx=flipdim(px,2); ppy=flipdim(py,2); |
| 116 |
% px=ppx; py=ppy; |
| 117 |
% |
| 118 |
%test cases for x,y locations: |
| 119 |
% ox=0; oy=6;%inside point : sum(ang2)=+-360 |
| 120 |
% ox=(px(1)+px(2))/2; oy=(py(1)+py(2))/2;%edge limit case : sum(ang2)=+-180 |
| 121 |
% ox=(px(1)+px(4))/2; oy=(py(1)+py(4))/2;%edge limit case : sum(ang2)=-180 |
| 122 |
% ox=px(1); oy=py(1);%corner limit case : sum(ang2)=NaN |
| 123 |
% ox=0; oy=12;%outside point : sum(ang2)=0 |
Parent Directory
| 
Revision Log
| 
Revision Graph