Equirectangular rotation, G'mic / GIMP / Krita

# G'mic script to be used with the G'mic plugin in GIMP or Krita.
# Paste this into a text-file. ".gmic" in your home directory on *nix or "user.gmic"
# in /Users/%username%/AppData/Roaming (maybe) on Windows 10 at least.



#@gui ____<b>Map Projection</b>
#-----------------------------
#@gmic rotate_equirectangular_map : roll angle, pitch angle, yaw angle
#@gmic : Take a map or equirectangular panorama
#@gmic : and rotate it around three axises

#@gui Rotate Equirectangular Map : rotate_equirectangular_map
#@gui : note = note("Rotate an equirectangular map or panorama around three axises in order in degrees.")
#@gui : Roll = float(0,-360,360)
#@gui : Pitch ("north and south" = float(0,-360,360)
#@gui : Yaw ("east and west") = float(0,-360,360)
#@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-23</i>.</small>")


rotate_equirectangular_map :

# Default to 0 degrees
-skip ${1=0},${2=0},${3=0}

# Black magic spell to make the preview look better???
_fF_zoomaccurate=1
_fF_zoomfactor=1


# Create new image, the ""-block is the mathematical 'function' that
# determines the value of the output pixel wherever we are iterating

100%,100%,1,100%,"

# Command parameters from the command line or the sliders.
# Variables from outside the ""-block have to be outside the ""s.
  yaw="{$3}";
  pitch="{$2}";
  roll="{$1}";

# Converting degrees into radians for the matrices later
# Those minus signs are just a direction correction.
  alpha=-roll * pi / 180;
  beta=-pitch * pi / 180;
  gamma=-yaw * pi / 180;

# Equirectangular coordinates into spherical radians coordinates
# phi <- height
# theta <- width
# radius rho=1, we are always working with a unit sphere / globe with radius 1
# The extra pi in "otheta=PI..." is to offset the axis towards the center of the image
    otheta=pi + (x/w)*2*pi;
    ophi=(y/h)*pi;
    rho=1;

# Convert to Cartesian x, y, z, we can leave out rho=1
#    cx=rho*sin(ophi)*cos(otheta);
#    cy=rho*sin(ophi)*sin(otheta);
#    cz=rho*cos(ophi);

    cx=sin(ophi)*cos(otheta);
    cy=sin(ophi)*sin(otheta);
    cz=cos(ophi);

# Do the transformation in Cartesian using rotation matrix maths magic

rotz=cz*( cos(alpha)*cos(beta) )
	+ cy*( cos(alpha)*sin(beta)*sin(gamma) - sin(alpha)*cos(gamma) )
	+ cx*( cos(alpha)*sin(beta)*cos(gamma) + sin(alpha)*sin(gamma) );

roty=cz*( sin(alpha)*cos(beta) )
    + cy*( sin(alpha)*sin(beta)*sin(gamma) + cos(alpha)*cos(gamma) )
    + cx*( sin(alpha)*sin(beta)*cos(gamma) - cos(alpha)*sin(gamma) );

rotx=cz*(-sin(beta))
    + cy*(cos(beta)*sin(gamma))
    + cx*(cos(beta)*cos(gamma));

# Convert the rotated point BACK into spherical coordinates
# Atan2 is a special programming thing to do arctan, except better
# It's got its own Wikipedia page

   ntheta=atan2(roty, rotx);
   nphi=acos( rotz / sqrt(rotx*rotx + roty*roty + rotz*rotz) );

# Convert spherical coordinates almost straight back into
# equirectangular pixel coordinates.

  x2=( (ntheta/pi/2)*w+w/2 );
  y2=(nphi/pi)*h;

# I: from image #0 take the vector(pixel value) of x2,y2,z(=0)
# "2": Cubic interpolation of pixels
# "0": Boundary condition: if you hit "nothing", what do you do? N/A here.

   I(#0,x2,y2,0,2,0)"

# Keep the last image generated ("k." = keep[-1] image index)

keep[-1]



#@gmic sinusoidal_map : yaw angle, pitch angle, roll angle, bgred, bggreen, bgblue, slices
#@gmic : Take a map or equirectangular panorama
#@gmic : and rotate it around z, y and x angles,
#@gmic : then project it as a sliced sinusoidal map

#@gui Sinusoidal Map : sinusoidal_map
#@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a sinusoidal projection.")
#@gui : Roll = float(0,-360,360)
#@gui : Pitch ("north and south" = float(0,-360,360)
#@gui : Yaw ("east and west") = float(0,-360,360)
#@gui : Background color = color(128,128,128)
#@gui : Slices = int(1,1,72)
#@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-23</i>.</small>")


sinusoidal_map :

# $1-3 Default to 0 degrees
# $4-6 color(r,g,b) gives you three separate parameters instead of, say, a vector
-skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=0} -check ${7=1}>=1

slices=$7

_fF_zoomaccurate=1
_fF_zoomfactor=1

# First rotate the map with the appropriate custom command
    rotate_equirectangular_map $1,$2,$3

# Slice the rotated map image into vertical sections
    split x,$slices

# Loop over each slice. $! = number of images in the image list
repeat $!

# No keyword: create an image. The ""-block is the mathematical function
# whose result/return value at the end determines the value of the pixel
# x, y we are iterating over.
	100%,100%,1,100%,"

        # Background color
        bgr="{$4}";
        bgg="{$5}";
        bgb="{$6}";

        # Half-accidental magical formula. Half-sin is used to center the coordinates somehow??
        halfsin=sin( (y/h)*pi )*w / 2;
        x2=(x+halfsin-w/2)/sin( y/h*pi);
        # y = y in this case, just use original y.

        # Check whether we're coloring inside or outside the sinusoidal shape
        if ( (x >= (w/2-halfsin) && x <= (w/2+halfsin) ),
            I(#"$>",x2,y,0,2,0),
            color=[bgr, bgg, bgb]);"
# repeat-loop end
done

# Delete the old map, which was sliced up
remove[0-{$slices-1}]

# join all the rest (new) slices together in order on the x axis
append x

# Keep the last (. = [-1] image index) image generated
keep[-1]