#############################################################################################################################
#How to use this script
#
#  1. Load a suitable equirectangular image in Gimp (with four walls in a rectangle
# - this script does not cover trapezoidal rooms)
#   2. Open the MathMap plug-in
#   3. Copy & paste the code that you will find below
#  4. Select "Wrap" in the settings
#   5. In the User values, change lambda0 until a corner of the room in your panorama is at the border of the image.
#   6. Still in the User values, change t3 until the leftmost wall is exactly covered by
#     the first rectilinear projection
#   7. Change t2 until the second leftmost wall is covered exactly by the second rectilinear
#      projection
#   8. Finally, change t3 until the last corner of the room is split between the last two
#      rectilinear projections
#   9. Change size (the vertical field of view) and yshift until the image looks fine
#  10. Press OK, wait, and save.
#####################################################
#    Four rectilinear projections covering 360°     #
#         Mathmap script by Seb Przd (2007)         #
#          http://flickr.com/photos/sbprzd          #
#            Licensed under the GPLv3               #
#####################################################
filter fourrectilinear (image in,float size: 1-2 (1),
float lambda0:-1-1 (0),float yshift:-1-1(0),float t3:
0-1(0.5),float t2:0-1(0.5),float t1:0-1(0.5))
il=0.01;ir=0.99;epsl=-100;epsr=100;iter=0;eps=1;while
(abs(eps)>0.0001)&&(iter<20) do xu=(il+ir)/2;tt1=tan(
pi/2-t1*pi)/2; tt3=tan(pi/2-t3*pi);yu=tt1+sqrt(xu*(1-
xu)+tt1^2); llmyy=tt3+sqrt(xu*(1-xu)+tt3^2);ll=llmyy+
yu;eps=tan(pi/2-t2*pi)-(-yu*llmyy+(1-xu)^2)/(ll*(1-xu
));if eps*epsl<0 then ir=xu;epsr=eps;else il=xu;epsl=
eps;end;iter=iter+1;end;l2=(1-xu)/(1+ll);l4=xu/(1+ll)
;l1=llmyy/(1+ll);l3=yu/(1+ll);xmvec=[-1,-1+l2+l4,0,l2
+l4,1]; xcvec=[-1+l4,-l3,l2,1-l1];lvec=[l1,l2,l3,l4];
lpv=[-3,-1,1,3]*pi/4+atan(l4,l1);xx=x*size;yy=y*size;
co=0; while co<=3 do xc=xcvec[co];xt=xx-xc*X;yt=yy-X*
yshift;if xx<=xmvec[co+1]*X&&xx>= xmvec[co]*X then lp
=lpv[co]; x0=xt;y0=yt;lc=lvec[co];end;co=co+1;end;rr=
sqrt(x0^2+y0^2);c=atan(rr/Y/2/lc);phi=asin(y0*sin(c)/
rr);in(xy:[atan(x0*tan(c),rr)*X/pi+lp*X/pi+lambda0*X,
phi*Y/(pi/2)]) end