Edit the Nadir and Zenith of an Equirectangular Panorama

Often the most difficult part of an equirectangular panorama (360° wide, 180° vertically) is the nadir, i.e. the bottom. This is where the tripod is, so no clear shot can be had without some sort of acrobatics to tilt the tripod away. So an important part of creating a nice equirectangular panorama is editing the nadir (and sometimes, the zenith).
The procedure I had been using so far is:

1. Create the equirectangular image (in Gimp)
   
2. Use hugin to create the rectilinear projection of the nadir, with a 90° field-of-view
3. Edit this rectilinear image (in Gimp)
   
4. Use hugin to output this edited rectilinear image into an equirectangular panorama
5. Manually blend the extra nadir unto the equirectangular image (in Gimp).

As you can see, this needs going back and forth between Gimp and hugin, and is not very practical. Spurred by the Photoshop tool mentioned here I bring to you a Gimp/Mathmap script to create the zenith and nadir of an equirectangular panorama, edit in place and transform back.

You should first make a copy of the layer with the equirectangular image, since the transformation to and back from the zenith/nadir does not keep the part around the horizon.

Filter used to create the Zenith and Nadir:

filter ToNadirZenith (image in)
# Filter created by Seb Przd
# Licensed under the GPL

if x<0 then
sinphi1=1;
xx=x+X/2;
else
sinphi1=-1;
xx=x-X/2;
end;

yy=y;

rr=sqrt(xx^2+yy^2);

c=atan(rr/Y);

phi =
if rr == 0 then
 0
else
 asin(cos(c)*sinphi1)
end;

xxx=atan(xx,-yy*sinphi1)*X/pi;
yyy=phi*Y/(pi/2);

in(xy:[xxx,yyy])
end

Filter used to transform back the Zenith and Nadir:

filter FromNadirZenith (image in)
# Filter created by Seb Przd
# Licensed under the GPL

output=1;
if y>Y/4 then
sinphi1=1;
xc=-X/2;
else if y<-Y/4 then
sinphi1=-1;
xc=X/2;
else
 output=0;
end;
end; 

cosc=sinphi1*sin(y/Y*pi/2);
xx=cos(y/Y*pi/2)*sin(x/X*pi)/cosc;
yy=-sinphi1*cos(y/Y*pi/2)*cos(x/X*pi)/cosc; 

if abs(xx)>1 then output=0; end;

if output then
in(xy:[xx*X/2+xc,yy*Y])
else
rgbaColor(0,0,0,0)
end

end