API tools faq
paste
Naeddyr

G'mic projection scripts - 2021-01-13

Naeddyr
Aug 27th, 2020 (edited)
326
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 60.46 KB | None | 0 0
raw download clone embed print report
  1. #@gui ____<b>Map Projection</b>
  2. #-----------------------------
  3. ###############################
  4. ## EQUIRECTANGULAR ROTATION
  5. #@gmic rotate_equirectangular_map : roll angle, pitch angle, yaw angle, revert bool
  6. #@gmic : Take a map or equirectangular panorama
  7. #@gmic : and rotate it around three axises
  8.  
  9. #@gui Rotate Equirectangular Map : rotate_equirectangular_map
  10. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order in degrees.")
  11. #@gui : Roll = float(0,-360,360)
  12. #@gui : Pitch ("north and south") = float(0,-360,360)
  13. #@gui : Yaw ("east and west") = float(0,-360,360)
  14. #@gui : Revert rotations = bool(false)
  15. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  16.  
  17.  
  18. rotate_equirectangular_map :
  19.  
  20. # Default to 0 degrees
  21. -skip ${1=0},${2=0},${3=0},${4=0}
  22.  
  23. # Black magic spell to make the preview look better???
  24. _fF_zoomaccurate=1
  25. _fF_zoomfactor=1
  26.  
  27.  
  28. # Create new image, the ""-block is the mathematical 'function' that
  29. # determines the value of the output pixel wherever we are iterating
  30.  
  31. -input 100%,100%,1,100%,"
  32.  
  33. # Command parameters from the command line or the sliders.
  34. # Variables from outside the ""-block have to be outside the ""s.
  35. yaw="{$3}";
  36. pitch="{$2}";
  37. roll="{$1}";
  38. revert="{$4}";
  39.  
  40. # Converting degrees into radians for the matrices later
  41. # Those minus signs are just a direction correction.
  42. alpha=-roll * pi / 180;
  43. beta=-pitch * pi / 180;
  44. gamma=-yaw * pi / 180;
  45.  
  46. # Equirectangular coordinates into spherical radians coordinates
  47. # phi <- height
  48. # theta <- width
  49. # radius rho=1, we are always working with a unit sphere / globe with radius 1
  50. # The extra pi in "otheta=PI..." is to offset the axis towards the center of the image
  51. otheta=pi + (x/w)*2*pi;
  52. ophi=(y/h)*pi;
  53. rho=1;
  54.  
  55. # Convert to Cartesian x, y, z, we can leave out rho=1
  56. # cx=rho*sin(ophi)*cos(otheta);
  57. # cy=rho*sin(ophi)*sin(otheta);
  58. # cz=rho*cos(ophi);
  59.  
  60. cx=sin(ophi)*cos(otheta);
  61. cy=sin(ophi)*sin(otheta);
  62. cz=cos(ophi);
  63.  
  64. # Do the transformation in Cartesian using rotation matrix maths magic
  65. if( revert==0,
  66.  
  67. rotz=cz*( cos(alpha)*cos(beta) )
  68. + cy*( cos(alpha)*sin(beta)*sin(gamma) - sin(alpha)*cos(gamma) )
  69. + cx*( cos(alpha)*sin(beta)*cos(gamma) + sin(alpha)*sin(gamma) );
  70.  
  71. roty=cz*( sin(alpha)*cos(beta) )
  72. + cy*( sin(alpha)*sin(beta)*sin(gamma) + cos(alpha)*cos(gamma) )
  73. + cx*( sin(alpha)*sin(beta)*cos(gamma) - cos(alpha)*sin(gamma) );
  74.  
  75. rotx=cz*(-sin(beta))
  76. + cy*(cos(beta)*sin(gamma))
  77. + cx*(cos(beta)*cos(gamma)); ,
  78.  
  79. # !!!!!!!!!!!!!!!! Reversion matrix??
  80. # small letter x = cos(x)
  81. # capital letter X = sin(x)
  82. # {{ a*b, -A*b, B }
  83. # { B*G*a+A*g, -A*B*G+a*g, -G*b }
  84. # { A*G-B*a*g, G*a+A*B*g, b*g }}
  85.  
  86. ## Testing reversion
  87.  
  88. rotz=cz*( cos(-alpha)*cos(-beta) )
  89. + cy*( -sin(-alpha)*cos(-beta) )
  90. + cx*( sin(-beta) );
  91.  
  92. roty=cz*( sin(-beta)*sin(-gamma)*cos(-alpha) + sin(-alpha)*cos(-gamma) )
  93. + cy*( -sin(-alpha)*sin(-beta)*sin(-gamma) + cos(-alpha)*cos(-gamma) )
  94. + cx*( -sin(-gamma)*cos(-beta) );
  95.  
  96. rotx=cz*(sin(-alpha)*sin(-gamma) - sin(-beta)*cos(-alpha)*cos(-gamma) )
  97. + cy*(sin(-gamma)*cos(-alpha) + sin(-alpha)*sin(-beta)*cos(-gamma) )
  98. + cx*(cos(-beta)*cos(-gamma));
  99. );
  100.  
  101. # Convert the rotated point BACK into spherical coordinates
  102. # Atan2 is a special programming thing to do arctan, except better
  103. # It's got its own Wikipedia page
  104.  
  105. ntheta=atan2(roty, rotx);
  106. nphi=acos( rotz / sqrt(rotx*rotx + roty*roty + rotz*rotz) );
  107.  
  108. # Convert spherical coordinates almost straight back into
  109. # equirectangular pixel coordinates.
  110.  
  111. x2=( (ntheta/pi/2)*w+w/2 );
  112. y2=(nphi/pi)*h;
  113.  
  114. # I: from image #0 take the vector(pixel value) of x2,y2,z(=0)
  115. # "2": Cubic interpolation of pixels
  116. # "0": Boundary condition: if you hit "nothing", what do you do? N/A here.
  117.  
  118. I(#0,x2,y2,0,2,1)"
  119.  
  120. # Keep the last image generated ("k." = keep[-1] image index)
  121.  
  122. keep[-1]
  123.  
  124.  
  125. ########################
  126. ## SINUSOIDAL
  127. #@gmic sinusoidal_map : yaw angle, pitch angle, roll angle, bgred, bggreen, bgblue, slices
  128. #@gmic : Take a map or equirectangular panorama
  129. #@gmic : and rotate it around z, y and x angles,
  130. #@gmic : then project it as a sliced sinusoidal map
  131.  
  132. #@gui Sinusoidal Map : sinusoidal_map
  133. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a sinusoidal projection.")
  134. #@gui : Roll = float(0,-360,360)
  135. #@gui : Pitch ("north and south") = float(0,-360,360)
  136. #@gui : Yaw ("east and west") = float(0,-360,360)
  137. #@gui : Background color = color(128,128,128)
  138. #@gui : Background opacity (%) = float(100,0,100)
  139. #@gui : Slices = int(1,1,72)
  140. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  141.  
  142.  
  143. sinusoidal_map :
  144.  
  145. # $1-3 Default to 0 degrees
  146. # $4-6 color(r,g,b) gives you three separate parameters instead of, say, a vector
  147. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=0},${7=100} -check ${8=1}>=1
  148.  
  149. slices=$8
  150. # Give input alpha-channel to prevent errors with transparency
  151. to_a[0]
  152.  
  153. _fF_zoomaccurate=1
  154. _fF_zoomfactor=1
  155.  
  156. # First rotate the map with the appropriate custom command
  157. if $1!=0" || "$2!=0" || "$3!=0
  158. rotate_equirectangular_map $1,$2,$3
  159. endif
  160.  
  161. # Slice the rotated map image into vertical sections
  162. split x,$slices
  163.  
  164. # Loop over each slice. $! = number of images in the image list
  165. repeat $!
  166.  
  167. # No keyword: create an image. The ""-block is the mathematical function
  168. # whose result/return value at the end determines the value of the pixel
  169. # x, y we are iterating over.
  170. -input 100%,100%,1,100%,"
  171.  
  172. # Background color
  173. bgr="{$4}";
  174. bgg="{$5}";
  175. bgb="{$6}";
  176. bgopacity="{$7}/100*255";
  177.  
  178. # Half-accidental magical formula. Half-sin is used to center the coordinates somehow??
  179. halfsin=sin( (y/h)*pi )*w / 2;
  180. x2=(x+halfsin-w/2)/sin( y/h*pi);
  181. # y = y in this case, just use original y.
  182.  
  183. # Check whether we're coloring inside or outside the sinusoidal shape
  184. if ( (x >= (w/2-halfsin) && x <= (w/2+halfsin) ),
  185. I(#"$>",x2,y,0,2,0),
  186. color=[bgr, bgg, bgb, bgopacity]);"
  187. # repeat-loop end
  188. done
  189.  
  190. # Delete the old map, which was sliced up
  191. remove[0-{$slices-1}]
  192.  
  193. # join all the rest (new) slices together in order on the x axis
  194. append x
  195.  
  196. # Keep the last (. = [-1] image index) image generated
  197. keep[-1]
  198.  
  199.  
  200.  
  201. #######################
  202. ## TRIANGULAR
  203. #@gmic triangular_projection : yaw angle, pitch angle, roll angle, bgred, bggreen, bgblue, slices
  204. #@gmic : Take a map or equirectangular panorama
  205. #@gmic : and rotate it around z, y and x angles,
  206. #@gmic : then project it as a sliced triangular map.
  207. #@gmic : This is not the Collignon projection.
  208. #@gmic : It's a naive triangular interpolation.
  209.  
  210. #@gui Triangular projection : triangular_projection
  211. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a naive rectilinear projection. This is not a Collignon projection, but a naive linear interpolation.")
  212. #@gui : note = note("When the width at the poles is 50%, this becomes the Eckert I projection. These projections are not equal-area nor do they have any other useful cartographical property.")
  213. #@gui : Roll = float(0,-360,360)
  214. #@gui : Pitch ("north and south") = float(0,-360,360)
  215. #@gui : Yaw ("east and west") = float(0,-360,360)
  216. #@gui : Background color = color(128,128,128)
  217. #@gui : Background opacity (%) = float(100,0,100)
  218. #@gui : Slices = int(1,1,72)
  219. #@gui : Width at the poles % = float(50,0,100)
  220. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  221.  
  222.  
  223. triangular_projection :
  224.  
  225.  
  226. # $1-3 Default to 0 degrees
  227. # $4-7 color(r,g,b) gives you three separate parameters instead of, say, a vector
  228. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=0},${7=100} -check ${8=1}>=1" && "${9=5}>=0
  229.  
  230. slices=$8
  231. magic=$9
  232.  
  233. # Give input alpha-channel to prevent errors with transparency
  234. to_a[0]
  235.  
  236.  
  237. # Background color
  238. bgr=$4;
  239. bgg=$5;
  240. bgb=$6;
  241. bgopacity={$7/100*255};
  242.  
  243. # Give input alpha-channel to prevent errors with transparency
  244. to_a[0]
  245.  
  246. _fF_zoomaccurate=1
  247. _fF_zoomfactor=1
  248.  
  249. # First rotate the map with the appropriate custom command
  250. if $1!=0" || "$2!=0" || "$3!=0
  251. rotate_equirectangular_map $1,$2,$3
  252. endif
  253.  
  254. # Slice the rotated map image into vertical sections
  255. split x,$slices
  256.  
  257. # Input: create an image. The ""-block is the mathematical function
  258. # whose result/return value at the end determines the value of the pixel
  259. # x, y we are iterating over.
  260. -input 100%,100%,1,4,"
  261.  
  262.  
  263. startwidth="$magic";
  264.  
  265. if( (y<=h/2),
  266. my=y,
  267. my=h-y);
  268.  
  269. sc=startwidth/100;
  270. squish=w / lerp((w*sc), w, my/(h/2));
  271. shift=(w-squish*w)/2;
  272. X = squish * x +shift;
  273. Y=y;
  274.  
  275. #debug normalisation values: cap at w so everything normalizes
  276. #result=[min(w,X),min(Y,w),0,w];"
  277.  
  278. result=[X,Y,0,0];"
  279.  
  280. # Loop over each slice. $! = number of images in the image list
  281. # Because the warp map is identical for each "slice" in this version
  282. # of the script / process, it doesn't need to generated for each slice,
  283. # so it's just reused. The sinusoidal script does an in-place operation
  284. # for each slice.
  285.  
  286. repeat $slices
  287. # the $> is a loop index incrementing forwards (think i++).
  288. # if there was $< it would loop backwards (like i--)
  289. # but here we only need to use the one warp image at the very
  290. # end "[-1]" of the image stack
  291. -warp[$>] [-1],0,2,0
  292. done
  293.  
  294. #DISPLACEMENT WARP
  295. # Delete the warp map, which was sliced up
  296. remove[-1]
  297.  
  298. # join all the rest (new) slices together in order on the x axis
  299.  
  300. #for debugging by normalizing and appending the warp field image
  301. #-normalize[-1] 0,255
  302.  
  303. append x
  304.  
  305. # generate blank image, then fill_color with bg-color
  306. -input 100%,100%,1,4,0
  307. -fill_color[-1] $bgr,$bgg,$bgb,$bgopacity
  308.  
  309. reverse
  310. -blend alpha
  311.  
  312. # Keep the last (. = [-1] image index) image generated
  313. -keep[-1]
  314.  
  315.  
  316. #######################
  317. ## CYLINDRICAL EQUAL AREA
  318. #@gmic cylindrical_equal_area : yaw angle, pitch angle, roll angle, standard_parallel angle
  319. #@gmic : Take a map or equirectangular panorama
  320. #@gmic : and rotate it around z, y and x angles,
  321. #@gmic : then project it as a cylindrical equal area map.
  322.  
  323. #@gui Cylindrical equal-area projection : cylindrical_equal_area
  324. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a cylindrical equal-area projection.")
  325. #@gui : Roll = float(0,-360,360)
  326. #@gui : Pitch ("north and south") = float(0,-360,360)
  327. #@gui : Yaw ("east and west") = float(0,-360,360)
  328. #@gui : Standard parallels = float(45,0,80)
  329. #@gui : Named specializations = choice("[Standard parallels ↑]","Lambert 0° π:1","Behrmann 30°","Smyth/Craster ≈37°04\′17\″ 2:1","Trystan Edwards 37°24\′","Hobo-Dyer 37°30\′","Gall-Peters 45°","Balthasar 50°","Tobler ≈55°39\′14\″ 1:1","[Custom formula ↓]")
  330. #@gui : Custom formula = text{"acos(sqrt(1/pi))"}
  331. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  332.  
  333.  
  334. cylindrical_equal_area :
  335.  
  336. # $1-3 Default to 0 degrees
  337. -skip ${1=0},${2=0},${3=0} -check ${4=45}<=80
  338.  
  339. # Give input alpha-channel to prevent errors with transparency
  340. to_a[0]
  341.  
  342. # Standard parallells custom
  343. phi_input=$4
  344.  
  345. # G'mic-style "array"
  346. if $5==1 phi0=0
  347. elif $5==2 phi0={30*pi/180}
  348. elif $5==3 phi0={acos(sqrt(2/pi))}
  349. elif $5==4 phi0={37.4*pi/180}
  350. elif $5==5 phi0={37.5*pi/180}
  351. elif $5==6 phi0={45*pi/180}
  352. elif $5==7 phi0={50*pi/180}
  353. elif $5==8 phi0={acos(sqrt(1/pi))}
  354. elif $5==9 phi0={$6}
  355. # use Standards parallels input
  356. else phi0={$phi_input*pi/180}
  357. endif
  358. # phi0 is the input parameter Standard parallels in radians
  359.  
  360.  
  361.  
  362.  
  363. _fF_zoomaccurate=1
  364. _fF_zoomfactor=1
  365.  
  366. # First rotate the map with the appropriate custom command
  367. if $1!=0" || "$2!=0" || "$3!=0
  368. rotate_equirectangular_map $1,$2,$3
  369. endif
  370.  
  371.  
  372. # width of original map
  373. nw={0,w}
  374.  
  375. # height of original map
  376. nh={0,h}
  377.  
  378. # Determine the dimensions of the output projection by using
  379. # the ratio of w:h as give by Wikipedia
  380. # https://en.wikipedia.org/wiki/Cylindrical_equal-area_projection
  381. # w:h = pi*(cos phi_0)²
  382. W={$nw*pi*cos($phi0)^2}
  383. H={$nw}
  384.  
  385. # Initialize the image using the above precalculated width and height
  386. # The "" block contains the mathematical 'formula' that determines
  387. # what each x,y pixel is.
  388.  
  389. -input $W,$H,1,4,"
  390.  
  391. phi0="$phi0";
  392. pi2=2*pi;
  393.  
  394. # Use centered coordinates. 0,0 is in the middle of the
  395. # projection
  396. cx=x-w/2;
  397. cy=y-h/2;
  398.  
  399. # Turn them into radians??
  400. cx = (cx/w *pi2);
  401. cy = (cy/w *pi2);
  402.  
  403. # Somehow I cancelled out the formula for nX, it's
  404. # now part of the output image width formula above????
  405. nX = cx;
  406. # I don't know! I don't know! But this formula works!!
  407. nY = asin(cy * cos(phi0)^2);
  408.  
  409. # Here's the magic trick: the X and Y coordinates we are
  410. # calculating are not about where something goes, it's about
  411. # where something COMES FROM. Because of the way G'mic makes
  412. # you think about 'filters', you have to invert the transformation
  413. # direction: you look at the OUTPUT images' x and y coordinates,
  414. # then figure out where *from* that values has to come from.
  415. # Basically, we're taking the output and calculating it in reverse,
  416. # so that we get equirectangular (=spherical) coordinates back.
  417. # That's why instead of calculating with SINE, which is the way you
  418. # do it forward, we calculate with ARCSINE, which is backwards.
  419. # This is the number one reason why doing this stuff in G'mic is torture.
  420. # Check out the formulas in
  421. # https://mathworld.wolfram.com/CylindricalEqual-AreaProjection.html
  422.  
  423. # EQUIRECTANGULAR COORDINATES
  424. # Get the width and height of the original equirectangular input
  425. ow="{0,w}";
  426. oh="{0,h}";
  427.  
  428. # Pretty simply translate radians coordinates straight into
  429. # equirectangular pixel coordinates.
  430. X=ow/2 + (ow * nX/pi2) ;
  431. Y=oh/2 + (oh * nY/pi) ;
  432.  
  433. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2),
  434. # TAKE the value of X,Y FROM the original equirectangular source map
  435. # and input the value into x,y in the output.
  436. result=I(#0,X%ow,Y,0,2,0),
  437. # else output transparent pixel
  438. result=[0,0,0,0]
  439. );
  440.  
  441. result"
  442.  
  443. -autocrop[1]
  444.  
  445. # Keep the last (. = [-1] image index) image generated
  446. -keep[-1]
  447.  
  448.  
  449.  
  450. #######################
  451. ## LAMBERT AZIMUTHAL
  452. #@gmic lambert_azimuthal_projection : yaw angle, pitch angle, roll angle
  453. #@gmic : Take a map or equirectangular panorama
  454. #@gmic : and rotate it around z, y and x angles,
  455. #@gmic : then project it as a Lambert azimuthal equal area map.
  456.  
  457. #@gui Lambert Azimuthal Equal-Area : lambert_azimuthal_projection
  458. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Lambert Azimuthal Equal-Area map.")
  459. #@gui : Roll = float(0,-360,360)
  460. #@gui : Pitch = float(0,-360,360)
  461. #@gui : Yaw = float(0,-360,360)
  462. #@gui : Central longitude = float(0,-180,180)
  463. #@gui : Standard parallel = float(0,-90,90)
  464. #@gui : Cut-off at = float(180,0,180)
  465. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  466.  
  467.  
  468. lambert_azimuthal_projection :
  469.  
  470. # $1-3 Default to 0 degrees
  471. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=180}
  472.  
  473. phi0=0*pi/180
  474. cutoff={$6*pi/180}
  475.  
  476. # Give input alpha-channel to prevent errors with transparency
  477. to_a[0]
  478.  
  479.  
  480. _fF_zoomaccurate=1
  481. _fF_zoomfactor=1
  482.  
  483. # First rotate the map with the appropriate custom command
  484. if $1!=0" || "$2!=0" || "$3!=0
  485. rotate_equirectangular_map $1,$2,$3
  486. endif
  487.  
  488. # There is a bug in the naive formulas I put into this
  489. # script. Changing phi0 causes the map to get a literal
  490. # rip hole into it. If that was the only way for me to
  491. # change the angle of the globe view, this script
  492. # would be basically pretty useless and would only make
  493. # maps in one orientation. Happily, you can subsitute
  494. # those calculations with the usual equirectangular
  495. # rotation command. In this case, the two below are
  496. # needed to make "parallel" and "longitude" make sense.
  497. # The equirectangular_projection command applies the
  498. # rotations in a specific order. If you try to use it
  499. # *like* longitudes and parallels, you have to find the
  500. # equivalent values somehow. Below, to easily get the
  501. # "right" orientation, first East-West movement is
  502. # applied, then North-South.
  503. if $4!=0
  504. rotate_equirectangular_map 0,0,$4
  505. endif
  506.  
  507. if $5!=0
  508. rotate_equirectangular_map 0,{-$5},0
  509. endif
  510.  
  511.  
  512. # width of original map
  513. nw={0,w}
  514.  
  515. # height of original map
  516. nh={0,h}
  517.  
  518. # Initialize the image using the above precalculated width and height
  519. # The "" block contains the mathematical 'formula' that determines
  520. # what each x,y pixel is.
  521.  
  522. -input {$nw+1},{$nw+1},1,4,"
  523.  
  524. #inside the ""-block it's easier to pre-translate dollar-variables
  525. phi0="$phi0";
  526. cutoff="$cutoff";
  527. pi2=2*pi;
  528.  
  529.  
  530. # Use centered coordinates. 0,0 is in the middle of the
  531. # projection
  532. cx=x-w/2;
  533. cy=y-h/2;
  534.  
  535.  
  536. # Turn them into radians??
  537. cx = (cx/w *pi2);
  538. cy = (cy/w *pi2);
  539.  
  540. # For convenience sake. Rho is obviously, in retrospect
  541. # distance (from center)... At rho==sqrt/2 is where the
  542. # switch in the atan-function should happen.
  543. # These are from
  544. # https://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html
  545. rho=sqrt(cx^2 + cy^2);
  546. C=2*asin(1/2*rho);
  547.  
  548. # atan2 can be used here! The formula looks complicated, but just
  549. # take the parts of the fraction and replace the / with , for the function
  550. nX = atan2( (cx*sin(C)), ( (rho*cos(phi0)*cos(C))-(cy*sin(phi0)*sin(C))) );
  551.  
  552. nY = asin( cos(C)*sin(phi0)+(cy*sin(C)*cos(phi0)) / rho);
  553. # the sin(phi0)'s cancel out with multiplying by zero, should remove...
  554.  
  555. # For the logic behind inverting the formula, check out the
  556. # Cylindrical Equal-Area filter for a rant
  557.  
  558. # EQUIRECTANGULAR COORDINATES
  559. # Get the width and height of the original equirectangular input
  560. ow="{0,w}";
  561. oh="{0,h}";
  562.  
  563. # Pretty simply translate radians coordinates straight into
  564. # equirectangular pixel coordinates.
  565. X=ow/2 + (ow * nX/pi2);
  566. Y=oh/2 + (oh * nY/pi);
  567.  
  568. # Calculate the map's rendering radius based on user's cutoff angle
  569. # This is a simplified Great Circle distance formula
  570. # where half the formula disappears because sin(0)=0
  571. deltaS=acos(cos(nY)*cos(nX));
  572.  
  573. # check whether the pixel's geographical nX, nY coordinates
  574. # are within the normal range AND inside the deltaS cutoff
  575. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2 && abs(deltaS)<=cutoff),
  576. # TAKE the value of X,Y FROM the original equirectangular source map
  577. # and input the value into x,y in the output.
  578. result=I(#0,X,Y,0,2,0),
  579. # else output transparent pixel
  580. result=[0,0,0,0]
  581. );
  582. # return the last value
  583. result;"
  584.  
  585. # New trick: if you can't calculate the width and height of the image
  586. # in advance, doing it on a "blank" square image can work as well.
  587. # In the ""-block of the -input above the last pixel assignment
  588. # checks if the geographical coordinates we're looking for make
  589. # any sense, that is, no weird values above or below plus-minus 180 or
  590. # 90 degrees. If the pixel x,y coordinates would result in something
  591. # nonsensical, output a transparent pixel. This is necessary, because it's
  592. # not automatic: you get weird blocks of white *and* transparent
  593. # lines. Then you can use autocrop from g'mic to do some magic and remove
  594. # the transparent areas.
  595. -autocrop[-1]
  596.  
  597.  
  598. # Keep the last (. = [-1] image index) image generated
  599. -keep[-1]
  600.  
  601.  
  602. #######################
  603. ## MOLLWEIDE
  604. #@gmic mollweide_projection : yaw angle, pitch angle, roll angle, central_longitude angle, standard_parallel angle
  605. #@gmic : Take a map or equirectangular panorama
  606. #@gmic : and rotate it around z, y and x angles,
  607. #@gmic : then project it as a Mollweide projection map.
  608.  
  609. #@gui Mollweide : mollweide_projection
  610. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Mollweide equal-area pseudocylindrical map.")
  611. #@gui : Roll = float(0,-360,360)
  612. #@gui : Pitch = float(0,-360,360)
  613. #@gui : Yaw = float(0,-360,360)
  614. #@gui : Central longitude = float(0,-180,180)
  615. #@gui : Standard parallel = float(0,-90,90)
  616. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  617.  
  618.  
  619. mollweide_projection :
  620.  
  621. # $1-3 Default to 0 degrees
  622. -skip ${1=0},${2=0},${3=0},${4=0},${5=0}
  623.  
  624. # Give input alpha-channel to prevent errors with transparency
  625. to_a[0]
  626.  
  627. _fF_zoomaccurate=1
  628. _fF_zoomfactor=1
  629.  
  630. # First rotate the map with the appropriate custom command
  631. if $1!=0" || "$2!=0" || "$3!=0
  632. rotate_equirectangular_map $1,$2,$3
  633. endif
  634.  
  635. # To apply more standard "longitude" and "standard parallel
  636. # translations, use the same rotate command, except in an order
  637. # where longitude (yaw) is applied first and latitude (pitch)
  638. # is applied second, in two different operations.
  639. if $4!=0
  640. rotate_equirectangular_map 0,0,$4
  641. endif
  642.  
  643. if $5!=0
  644. rotate_equirectangular_map 0,{-$5},0
  645. endif
  646.  
  647.  
  648. # width and height of original map
  649. nw={0,w}
  650. nh={0,h}
  651.  
  652. # Initialize the image using the above precalculated width and height
  653. # The "" block contains the mathematical 'formula' that determines
  654. # what each x,y pixel is.
  655.  
  656. -input {$nw+1},{$nw+1},1,4,"
  657.  
  658. #inside the ""-block it's easier to pre-translate dollar-variables
  659. pi2=2*pi;
  660.  
  661.  
  662. # Use centered coordinates. 0,0 is in the middle of the
  663. # projection
  664. cx=x-w/2;
  665. cy=y-h/2;
  666.  
  667.  
  668. # Turn them into radians
  669. cx = (cx/w *pi2);
  670. cy = (cy/w *pi2);
  671.  
  672.  
  673. # constants for the formulas. R=1 is the Radius of the globe,
  674. # I've put 1 here just because i couldn't find a convenient
  675. # fraction of pi, and it doesn't change the shape, just scale.
  676. # Theta is separated as a constant following the example of
  677. # the mathworld equations I found:
  678. # https://mathworld.wolfram.com/MollweideProjection.html
  679. R=1;
  680. theta=asin(cy/(sqrt(2)*R));
  681.  
  682. nX = (pi*cx) / (2*R*sqrt(2)*cos(theta));
  683. nY = asin( ((2*theta)+sin(2*theta)) / pi);
  684.  
  685.  
  686. # For the logic why I'm inverting the formula, check out the
  687. # Cylindrical Equal-Area filter for a rant
  688.  
  689. # EQUIRECTANGULAR COORDINATES
  690. # Get the width and height of the original equirectangular input
  691. ow="{0,w}";
  692. oh="{0,h}";
  693.  
  694. # Pretty simply translate radians coordinates straight into
  695. # equirectangular pixel coordinates.
  696. X=ow/2 + (ow * nX/pi2);
  697. Y=oh/2 + (oh * nY/pi);
  698.  
  699. # check whether the pixel's geographical nX, nY coordinates
  700. # are within the normal range
  701. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2),
  702. # then TAKE the value of X,Y FROM the original equirectangular
  703. # source map and input the value into x,y in the output.
  704. result=I(#0,X,Y,0,2,0),
  705. # else output a transparent pixel
  706. result=[0,0,0,0]
  707. );
  708. # return the last value
  709. result;"
  710.  
  711. # New trick: if you can't calculate the width and height of the image
  712. # in advance, doing it on a "blank" square image can work as well.
  713. # In the ""-block of the -input above the last pixel assignment
  714. # checks if the geographical coordinates we're looking for make
  715. # any sense, that is, no weird values above or below plus-minus 180 or
  716. # 90 degrees. If the pixel x,y coordinates would result in something
  717. # nonsensical, output a transparent pixel. This is necessary, because it's
  718. # not automatic: you get weird blocks of white *and* transparent
  719. # lines. Then you can use autocrop from g'mic to do some magic and remove
  720. # the transparent areas.
  721. -autocrop[-1]
  722.  
  723.  
  724. # Keep the last (. = [-1] image index) image generated
  725. -keep[-1]
  726.  
  727.  
  728. #######################
  729. ## ORTHOGRAPHIC
  730. #@gmic orthographic_projection : yaw angle, pitch angle, roll angle, central_longitude angle, standard_parallel angle
  731. #@gmic : Take a map or equirectangular panorama
  732. #@gmic : and rotate it around z, y and x angles,
  733. #@gmic : then project it as an orthographic projection map.
  734.  
  735. #@gui Orthographic : orthographic_projection
  736. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an orthographic hemisphere map. This map projection is basically the view you'd get if you looked at earth from an infinite distance.")
  737. #@gui : Roll = float(0,-360,360)
  738. #@gui : Pitch = float(0,-360,360)
  739. #@gui : Yaw = float(0,-360,360)
  740. #@gui : Central longitude = float(0,-180,180)
  741. #@gui : Standard parallel = float(0,-90,90)
  742. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  743.  
  744.  
  745. orthographic_projection :
  746.  
  747. # $1-3 Default to 0 degrees
  748. -skip ${1=0},${2=0},${3=0},${4=0},${5=0}
  749.  
  750. # Give input alpha-channel to prevent errors with transparency
  751. to_a[0]
  752.  
  753. _fF_zoomaccurate=1
  754. _fF_zoomfactor=1
  755.  
  756. # First rotate the map with the appropriate custom command
  757. if $1!=0" || "$2!=0" || "$3!=0
  758. rotate_equirectangular_map $1,$2,$3
  759. endif
  760.  
  761. # To apply more standard "longitude" and "standard parallel
  762. # translations, use the same rotate command, except in an order
  763. # where longitude (yaw) is applied first and latitude (pitch)
  764. # is applied second, in two different operations.
  765. if $4!=0
  766. rotate_equirectangular_map 0,0,$4
  767. endif
  768.  
  769. if $5!=0
  770. rotate_equirectangular_map 0,{-$5},0
  771. endif
  772.  
  773.  
  774. # width and height of original map
  775. nw={0,w}
  776. nh={0,h}
  777.  
  778. # Initialize the image using the above precalculated width and height
  779. # The "" block contains the mathematical 'formula' that determines
  780. # what each x,y pixel is.
  781.  
  782. -input {$nw+1},{$nw+1},1,4,"
  783.  
  784. #inside the ""-block it's easier to pre-translate dollar-variables
  785. pi2=2*pi;
  786.  
  787.  
  788. # Use centered coordinates. 0,0 is in the middle of the
  789. # projection
  790. cx=x-w/2;
  791. cy=y-h/2;
  792.  
  793.  
  794. # Turn them into radians
  795. cx = (cx/w *pi2);
  796. cy = (cy/w *pi2);
  797.  
  798.  
  799. # Formulas and constants:
  800. # https://en.wikipedia.org/wiki/Orthographic_projection_in_cartography
  801. # R, radius of globe seems to work best with R=pi/2
  802. R=pi/2;
  803. rho=sqrt(cx^2 + cy^2);
  804. C=asin(rho/R);
  805. phi0=0;
  806. l0=0;
  807.  
  808. # the formulas are simplified a bit because this script uses
  809. # phi0=0 and l0=0, leading to anything with sin(phi0) to
  810. # multiply with 0. The standard parallel and central longitude
  811. # are calculated earlier with actual global rotations.
  812. # Happily, here we can use atan2.
  813. nY = asin((cy*sin(C)*cos(phi0))/rho);
  814. nX = atan2((cx*sin(C)), (rho*cos(C)*cos(phi0)));
  815.  
  816. # For the logic why I'm inverting the formula, check out the
  817. # Cylindrical Equal-Area filter for a rant
  818.  
  819. # EQUIRECTANGULAR COORDINATES
  820. # Get the width and height of the original equirectangular input
  821. ow="{0,w}";
  822. oh="{0,h}";
  823.  
  824. # Pretty simply translate radians coordinates straight into
  825. # equirectangular pixel coordinates.
  826. X=ow/2 + (ow * nX/pi2);
  827. Y=oh/2 + (oh * nY/pi);
  828.  
  829. # check whether the pixel's geographical nX, nY coordinates
  830. # are within the normal range
  831. # Why does this break if the condition is set to true manually???
  832. if( (nX>-pi && nX<pi && nY>-pi/2 && nY<pi/2),
  833. # then TAKE the value of X,Y FROM the original equirectangular
  834. # source map and input the value into x,y in the output.
  835. result=I(#0,X,Y,0,2,0),
  836. # else output a transparent pixel
  837. result=[0,0,0,0]
  838. );
  839. # return the last value
  840. result;"
  841.  
  842. # New trick: if you can't calculate the width and height of the image
  843. # in advance, doing it on a "blank" square image can work as well.
  844. # In the ""-block of the -input above the last pixel assignment
  845. # checks if the geographical coordinates we're looking for make
  846. # any sense, that is, no weird values above or below plus-minus 180 or
  847. # 90 degrees. If the pixel x,y coordinates would result in something
  848. # nonsensical, output a transparent pixel. This is necessary, because it's
  849. # not automatic: you get weird blocks of white *and* transparent
  850. # lines. Then you can use autocrop from g'mic to do some magic and remove
  851. # the transparent areas.
  852. -autocrop[-1]
  853.  
  854.  
  855. # Keep the last (. = [-1] image index) image generated
  856. -keep[-1]
  857.  
  858. #######################
  859. ## AZIMUTHAL EQUIDISTANT
  860. #@gmic azimuthal_equidistant_projection : yaw angle, pitch angle, roll angle, central_longitude angle, standard_parallel angle
  861. #@gmic : Take a map or equirectangular panorama
  862. #@gmic : and rotate it around z, y and x angles,
  863. #@gmic : then project it as an azimuthal equidistant map.
  864.  
  865. #@gui Azimuthal equidistant : azimuthal_equidistant_projection
  866. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an azimuthal equidistant map. All points on the map are at proportionally correct distances from the center point, and that all points on the map are at the correct azimuth (direction) from the center point.")
  867. #@gui : Roll = float(0,-360,360)
  868. #@gui : Pitch = float(0,-360,360)
  869. #@gui : Yaw = float(0,-360,360)
  870. #@gui : Central longitude = float(0,-180,180)
  871. #@gui : Standard parallel = float(0,-90,90)
  872. #@gui : Cut-off at = float(180,0,180)
  873. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  874.  
  875.  
  876. azimuthal_equidistant_projection :
  877.  
  878. # $1-3 Default to 0 degrees
  879. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=180}
  880.  
  881. # Give input alpha-channel to prevent errors with transparency
  882. to_a[0]
  883.  
  884. cutoff={$6*pi/180}
  885.  
  886. _fF_zoomaccurate=1
  887. _fF_zoomfactor=1
  888.  
  889. # First rotate the map with the appropriate custom command
  890. if $1!=0" || "$2!=0" || "$3!=0
  891. rotate_equirectangular_map $1,$2,$3
  892. endif
  893.  
  894. # To apply more standard "longitude" and "standard parallel
  895. # translations, use the same rotate command, except in an order
  896. # where longitude (yaw) is applied first and latitude (pitch)
  897. # is applied second, in two different operations.
  898. if $4!=0
  899. rotate_equirectangular_map 0,0,$4
  900. endif
  901.  
  902. if $5!=0
  903. rotate_equirectangular_map 0,{-$5},0
  904. endif
  905.  
  906.  
  907. # width and height of original map
  908. nw={0,w}
  909. nh={0,h}
  910.  
  911. # Initialize the image using the above precalculated width and height
  912. # The "" block contains the mathematical 'formula' that determines
  913. # what each x,y pixel is.
  914.  
  915. -input {$nw+1},{$nw+1},1,4,"
  916.  
  917. #inside the ""-block it's easier to pre-translate dollar-variables
  918. pi2=2*pi;
  919. cutoff="$cutoff";
  920.  
  921.  
  922. # Use centered coordinates. 0,0 is in the middle of the
  923. # projection
  924. cx=x-w/2+1;
  925. cy=y-h/2+1;
  926.  
  927.  
  928. # Turn them into radians
  929. cx = (cx/w *pi2);
  930. cy = (cy/w *pi2);
  931.  
  932.  
  933. # Formulas and constants:
  934. # https://en.wikipedia.org/wiki/Orthographic_projection_in_cartography
  935. # R, radius of globe seems to work best with R=pi/2
  936. C=sqrt(cx^2 + cy^2);
  937. phi0=0;
  938. l0=0;
  939.  
  940. # the formulas are simplified a bit because this script uses
  941. # phi0=0 and l0=0, leading to anything with sin(phi0) to
  942. # multiply with 0. The standard parallel and central longitude
  943. # are calculated earlier with actual global rotations.
  944. # Happily, here we can use atan2.
  945. nY = asin( (cy*sin(C)) / C );
  946. nX = atan2((cx*sin(C)), (C*cos(C)));
  947.  
  948. # For the logic why I'm inverting the formula, check out the
  949. # Cylindrical Equal-Area filter for a rant
  950.  
  951. # EQUIRECTANGULAR COORDINATES
  952. # Get the width and height of the original equirectangular input
  953. ow="{0,w}";
  954. oh="{0,h}";
  955.  
  956. # Pretty simply translate radians coordinates straight into
  957. # equirectangular pixel coordinates.
  958. X=ow/2 + (ow * nX/pi2);
  959. Y=oh/2 + (oh * nY/pi);
  960.  
  961. # Calculate the map's rendering radius based on user's cutoff angle
  962. # This is a simplified Great Circle distance formula
  963. # where half the formula disappears because sin(0)=0
  964. deltaS=acos(cos(nY)*cos(nX));
  965.  
  966. # check whether the pixel's geographical nX, nY coordinates
  967. # are within the normal range AND inside the deltaS cutoff
  968. if( (nX>=-pi && nX<=pi && nY>=-pi/2 && nY<=pi/2 && abs(deltaS)<=cutoff),
  969. # TAKE the value of X,Y FROM the original equirectangular source map
  970. # and input the value into x,y in the output.
  971. result=I(#0,X,Y,0,2,0),
  972. # else output transparent pixel
  973. result=[0,0,0,0]
  974. );
  975. result"
  976.  
  977.  
  978. # generate blank image, then fill_color with bg-color
  979. -input 100%,100%,1,4,0
  980. -fill_color[-1] 0,0,0,0
  981. -shape_circle {-1,w}
  982. -normalize[-1] 0,255
  983. -negate[-1]
  984. +blend_fade[-2,-3] [-1]
  985.  
  986.  
  987. # New trick: if you can't calculate the width and height of the image
  988. # in advance, doing it on a "blank" square image can work as well.
  989. # In the ""-block of the -input above the last pixel assignment
  990. # checks if the geographical coordinates we're looking for make
  991. # any sense, that is, no weird values above or below plus-minus 180 or
  992. # 90 degrees. If the pixel x,y coordinates would result in something
  993. # nonsensical, output a transparent pixel. This is necessary, because it's
  994. # not automatic: you get weird blocks of white *and* transparent
  995. # lines. Then you can use autocrop from g'mic to do some magic and remove
  996. # the transparent areas.
  997. -autocrop[-1]
  998.  
  999.  
  1000. # Keep the last (. = [-1] image index) image generated
  1001. -keep[-1]
  1002.  
  1003.  
  1004.  
  1005.  
  1006.  
  1007.  
  1008.  
  1009. #######################
  1010. ## CONIC EQUIDISTANT
  1011. #@gmic conic_equidistant_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle
  1012. #@gmic : Take a map or equirectangular panorama
  1013. #@gmic : and rotate it around z, y and x angles,
  1014. #@gmic : then project it as a conic equidistant map.
  1015.  
  1016. #@gui Conic equidistant : conic_equidistant_projection
  1017. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a conic equidistant map. Neither equal-area or conformal.")
  1018. #@gui : Roll = float(0,-360,360)
  1019. #@gui : Pitch = float(0,-360,360)
  1020. #@gui : Yaw = float(0,-360,360)
  1021. #@gui : Central longitude = float(0,-180,180)
  1022. #@gui : Reference parallel = float(0,-90,90)
  1023. #@gui : Standard parallel #1 = float(20,-90,90)
  1024. #@gui : Standard parallel #2 = float(60,-90,90)
  1025. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1026.  
  1027.  
  1028. conic_equidistant_projection :
  1029.  
  1030. # $1-3 Default to 0 degrees
  1031. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=20},${7=60}
  1032.  
  1033. cutoff={$6*pi/180}
  1034. l0=$4
  1035. phi0=$5
  1036. phi1=$6
  1037. phi2=$7
  1038.  
  1039. # Give input alpha-channel to prevent errors with transparency
  1040. to_a[0]
  1041.  
  1042. _fF_zoomaccurate=1
  1043. _fF_zoomfactor=1
  1044.  
  1045. # First rotate the map with the appropriate custom command
  1046. if $1!=0" || "$2!=0" || "$3!=0
  1047. rotate_equirectangular_map $1,$2,$3
  1048. endif
  1049.  
  1050. # To apply more standard "longitude" and "standard parallel
  1051. # translations, use the same rotate command, except in an order
  1052. # where longitude (yaw) is applied first and latitude (pitch)
  1053. # is applied second, in two different operations.
  1054. # if $4!=0
  1055. # rotate_equirectangular_map 0,0,$4
  1056. # endif
  1057. #
  1058. # if $5!=0
  1059. # rotate_equirectangular_map 0,{-$5},0
  1060. # endif
  1061.  
  1062.  
  1063. # width and height of original map
  1064. nw={0,w}
  1065. nh={0,h}
  1066.  
  1067. # Initialize the image using the above precalculated width and height
  1068. # The "" block contains the mathematical 'formula' that determines
  1069. # what each x,y pixel is.
  1070.  
  1071. -input {$nw+1},{$nw+1},1,4,"
  1072.  
  1073. #inside the ""-block it's easier to pre-translate dollar-variables
  1074. pi2=2*pi;
  1075. phi0="$phi0"*pi/180;
  1076. l0="$l0"*pi/180;
  1077.  
  1078. phi1="$phi1"*pi/180;
  1079. phi2="$phi2"*pi/180;
  1080.  
  1081.  
  1082.  
  1083.  
  1084. # Use centered coordinates. 0,0 is in the middle of the
  1085. # projection
  1086. cx=x-w/2+1;
  1087. cy=y-h/2+1;
  1088.  
  1089.  
  1090. # Turn them into radians
  1091. cx = (cx/w *pi2*4/3);
  1092. cy = -(cy/w *pi2*4/3);
  1093.  
  1094.  
  1095. # Formulas and constants:
  1096. # https://mathworld.wolfram.com/ConicEquidistantProjection.html
  1097.  
  1098. if( (phi1==phi2),
  1099. n=sin(phi1),
  1100. n=(cos(phi1)-cos(phi2)) / (phi2-phi1)
  1101. );
  1102. G=cos(phi1)/n + phi1;
  1103. rho0=G-phi0;
  1104.  
  1105. rho=sign(n)*sqrt( cx^2 + (rho0-cy)^2);
  1106.  
  1107. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1108. # page 114, if using ATAN2, the signs must be reversed.
  1109. # I KNEW THERE WAS SOMETHING WRONG WITH THE FORMULA
  1110. # I had to dig deep.
  1111. if( (n>0),
  1112. (theta=atan2(cx, (rho0-cy))),
  1113. (theta=atan2(-cx, (-rho0+cy)))
  1114. );
  1115.  
  1116. nY = G-rho;
  1117. nX = l0 + theta/n;
  1118.  
  1119. # For the logic why I'm inverting the formula, check out the
  1120. # Cylindrical Equal-Area filter for a rant
  1121.  
  1122. # EQUIRECTANGULAR COORDINATES
  1123. # Get the width and height of the original equirectangular input
  1124. ow="{0,w}";
  1125. oh="{0,h}";
  1126.  
  1127. # Pretty simply translate radians coordinates straight into
  1128. # equirectangular pixel coordinates.
  1129. X=ow/2 + (ow * nX/pi2);
  1130. Y=oh/2 + (oh * -nY/pi);
  1131.  
  1132.  
  1133. # check whether the pixel's geographical nX, nY coordinates
  1134. # are within the normal range, here adjusted for l0 and phi0
  1135. # which are kind of superfluous here, they don't do much anything.
  1136.  
  1137. L=nX-l0;
  1138. P=nY-phi0;
  1139. if( (L>=-pi && L<=pi && P>=-pi/2 && P<=pi/2),
  1140. # TAKE the value of X,Y FROM the original equirectangular source map
  1141. # and input the value into x,y in the output.
  1142. result=I(#0,X%ow,Y,0,2,0),
  1143. # else output transparent pixel
  1144. result=[0,0,0,0]
  1145. );
  1146.  
  1147. result"
  1148.  
  1149.  
  1150. # New trick: if you can't calculate the width and height of the image
  1151. # in advance, doing it on a "blank" square image can work as well.
  1152. # In the ""-block of the -input above the last pixel assignment
  1153. # checks if the geographical coordinates we're looking for make
  1154. # any sense, that is, no weird values above or below plus-minus 180 or
  1155. # 90 degrees. If the pixel x,y coordinates would result in something
  1156. # nonsensical, output a transparent pixel. This is necessary, because it's
  1157. # not automatic: you get weird blocks of white *and* transparent
  1158. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1159. # the transparent areas.
  1160. autocrop[-1]
  1161.  
  1162.  
  1163. # Keep the last (. = [-1] image index) image generated
  1164. -keep[-1]
  1165.  
  1166.  
  1167.  
  1168. #######################
  1169. ## LAMBERT CONFORMAL CONIC
  1170. #@gmic lambert_conformal_conic_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1171. #@gmic : Take a map or equirectangular panorama
  1172. #@gmic : and rotate it around z, y and x angles,
  1173. #@gmic : then project it as a Lambert conic conformal map.
  1174.  
  1175. #@gui Lambert Conformal Conic : lambert_conformal_conic_projection
  1176. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Lambert conic conformal map.")
  1177. #@gui : Roll = float(0,-360,360)
  1178. #@gui : Pitch = float(0,-360,360)
  1179. #@gui : Yaw = float(0,-360,360)
  1180. #@gui : Central longitude = float(0,-180,180)
  1181. #@gui : Reference parallel = float(0,-90,90)
  1182. #@gui : Standard parallel #1 = float(20,-90,90)
  1183. #@gui : Standard parallel #2 = float(60,-90,90)
  1184. #@gui : Crop latitude to = float(180,0,180)
  1185. #@gui : Crop longitude to = float(360,0,360)
  1186. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1187.  
  1188.  
  1189. lambert_conformal_conic_projection :
  1190.  
  1191. # $1-3 Default to 0 degrees
  1192. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=20},${7=60},${8=180},${9=360}
  1193.  
  1194. latcutoff={$8*pi/180}
  1195. longcutoff={$9*pi/180}
  1196.  
  1197. l0=$4
  1198. phi0=$5
  1199. phi1=$6
  1200. phi2=$7
  1201.  
  1202. # Give input alpha-channel to prevent errors with transparency
  1203. to_a[0]
  1204.  
  1205. _fF_zoomaccurate=1
  1206. _fF_zoomfactor=1
  1207.  
  1208. # First rotate the map with the appropriate custom command
  1209. if $1!=0" || "$2!=0" || "$3!=0
  1210. rotate_equirectangular_map $1,$2,$3
  1211. endif
  1212.  
  1213. # width and height of original map
  1214. nw={0,w}
  1215. nh={0,h}
  1216.  
  1217. # Initialize the image using the above precalculated width and height
  1218. # The "" block contains the mathematical 'formula' that determines
  1219. # what each x,y pixel is.
  1220.  
  1221. -input {$nw+1},{$nw+1},1,4,"
  1222.  
  1223. #inside the ""-block it's easier to pre-translate dollar-variables
  1224. pi2=2*pi;
  1225. phi0="$phi0"*pi/180;
  1226. l0="$l0"*pi/180;
  1227.  
  1228. phi1="$phi1"*pi/180;
  1229. phi2="$phi2"*pi/180;
  1230.  
  1231. latcutoff="$latcutoff";
  1232. longcutoff="$longcutoff";
  1233.  
  1234.  
  1235.  
  1236.  
  1237. # Use centered coordinates. 0,0 is in the middle of the
  1238. # projection
  1239. cx=x-w/2+1;
  1240. cy=y-h/2+1;
  1241.  
  1242.  
  1243. # Turn them into radians
  1244. cx = (cx/w *pi2*pi/2);
  1245. cy = -(cy/w *pi2*pi/2);
  1246.  
  1247.  
  1248. # Formulas and constants:
  1249. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1250.  
  1251. if((phi1==phi2),
  1252. n=sin(phi1),
  1253. n= log(cos(phi1)/cos(phi2)) / log(tan((1/4)*pi+(1/2)*phi2)/tan((1/4)*pi+(1/2)*phi1))
  1254. );
  1255.  
  1256. F=(cos(phi1)*tan((1/4)*pi+(1/2)*phi1)^n) / n;
  1257.  
  1258. rho0=F/tan(1/4*pi+1/2*phi0)^n;
  1259.  
  1260. rho=sign(n)*sqrt( cx^2 + (rho0-cy)^2);
  1261.  
  1262. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1263. # page 114, if using ATAN2, the signs must be reversed.
  1264.  
  1265. if( (n>0),
  1266. (theta=atan2(cx, (rho0-cy))),
  1267. (theta=atan2(-cx, (-rho0+cy)))
  1268. );
  1269.  
  1270. # When the reference parallel is in the southern hemisphere,
  1271. # flip the projection and map
  1272. if( (phi0<0),
  1273. (n=-n)
  1274. );
  1275.  
  1276.  
  1277. if( (rho==0),
  1278. (nY = sign(n)*90*pi/180),
  1279. (nY = 2*atan((F/rho)^(1/n))-1/2*pi)
  1280. );
  1281. nX = l0 + theta/n;
  1282.  
  1283. # For the logic why I'm using inverse formulas, check out the
  1284. # Cylindrical Equal-Area filter
  1285.  
  1286. # EQUIRECTANGULAR COORDINATES
  1287. # Get the width and height of the original equirectangular input
  1288. ow="{0,w}";
  1289. oh="{0,h}";
  1290.  
  1291. # Pretty simply translate radians coordinates straight into
  1292. # equirectangular pixel coordinates.
  1293. X=ow/2 + (ow * nX/pi2);
  1294. Y=oh/2 + (oh * -nY/pi);
  1295.  
  1296.  
  1297. # check whether the pixel's geographical nX, nY coordinates
  1298. # are within the normal range, then further crop to the angles
  1299. # given by the crop values.
  1300.  
  1301. L=nX;
  1302. P=nY;
  1303. if( (L>=l0-longcutoff/2 && L<=l0+longcutoff/2 && P>=phi0-latcutoff/2 && P<=phi0+latcutoff/2),
  1304. # TAKE the value of X,Y FROM the original equirectangular source map
  1305. # and input the value into x,y in the output.
  1306. result=I(#0,X%ow,Y,0,2,0),
  1307. # else output transparent pixel
  1308. result=[0,0,0,0]
  1309. );
  1310.  
  1311. result"
  1312.  
  1313. # If the reference parallel is in the southern hemisphere,
  1314. # we've flipped the projection. Here we need to further rotate
  1315. # the image so it's not upside down.
  1316. if $phi0<0
  1317. rotate[-1] 180
  1318. endif
  1319.  
  1320.  
  1321. # New trick: if you can't calculate the width and height of the image
  1322. # in advance, doing it on a "blank" square image can work as well.
  1323. # In the ""-block of the -input above the last pixel assignment
  1324. # checks if the geographical coordinates we're looking for make
  1325. # any sense, that is, no weird values above or below plus-minus 180 or
  1326. # 90 degrees. If the pixel x,y coordinates would result in something
  1327. # nonsensical, output a transparent pixel. This is necessary, because it's
  1328. # not automatic: you get weird blocks of white *and* transparent
  1329. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1330. # the transparent areas.
  1331. autocrop[-1]
  1332.  
  1333.  
  1334. # Keep the last (. = [-1] image index) image generated
  1335. -keep[-1]
  1336.  
  1337.  
  1338. #######################
  1339. ## MERCATOR
  1340. #@gmic mercator : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1341. #@gmic : Take a map or equirectangular panorama
  1342. #@gmic : and rotate it around z, y and x angles,
  1343. #@gmic : then project it as a Mercator projection.
  1344. #@gui Mercator projection : mercator
  1345. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as a Mercator projection. You can get a Transverse Mercator by rolling the globe by 90")
  1346. #@gui : Roll = float(0,-360,360)
  1347. #@gui : Pitch = float(0,-360,360)
  1348. #@gui : Yaw = float(0,-360,360)
  1349. #@gui : Central longitude = float(0,-180,180)
  1350. #@gui : Crop latitude to = float(85,0,85)
  1351. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1352.  
  1353.  
  1354. mercator :
  1355.  
  1356. # $1-3 Default to 0 degrees
  1357. -skip ${1=0},${2=0},${3=0},${4=0},${5=85}
  1358.  
  1359. latcutoff={$5*pi/180}
  1360.  
  1361. l0=$4
  1362.  
  1363. # Give input alpha-channel to prevent errors with transparency
  1364. to_a[0]
  1365.  
  1366. _fF_zoomaccurate=1
  1367. _fF_zoomfactor=1
  1368.  
  1369. # First rotate the map with the appropriate custom command
  1370. if $1!=0" || "$2!=0" || "$3!=0
  1371. rotate_equirectangular_map $1,$2,$3
  1372. endif
  1373.  
  1374. # width and height of original map
  1375. nw={0,w}
  1376. nh={0,h}
  1377.  
  1378. # Initialize the image using the above precalculated width and height
  1379. # The "" block contains the mathematical 'formula' that determines
  1380. # what each x,y pixel is.
  1381.  
  1382. -input {$nw+1},{$nw+1},1,4,"
  1383.  
  1384. #inside the ""-block it's easier to pre-translate dollar-variables
  1385. pi2=2*pi;
  1386. l0="$l0"*pi/180;
  1387.  
  1388. latcutoff="$latcutoff";
  1389.  
  1390. # Use centered coordinates. 0,0 is in the middle of the
  1391. # projection
  1392. cx=x-w/2+1;
  1393. cy=y-h/2+1;
  1394.  
  1395.  
  1396. # Turn them into radians
  1397. cx = (cx/w *pi2);
  1398. cy = -( (cy/w *pi2) );
  1399.  
  1400.  
  1401. # Formulas and constants:
  1402. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1403.  
  1404. nY=atan(sinh(cy));
  1405. nX=l0+cx;
  1406.  
  1407. # For the logic why I'm using inverse formulas, check out the
  1408. # Cylindrical Equal-Area filter
  1409.  
  1410. # EQUIRECTANGULAR COORDINATES
  1411. # Get the width and height of the original equirectangular input
  1412. ow="{0,w}";
  1413. oh="{0,h}";
  1414.  
  1415. # Pretty simply translate radians coordinates straight into
  1416. # equirectangular pixel coordinates.
  1417. X=ow/2 + (ow * nX/pi2);
  1418. Y=oh/2 + (oh * -nY/pi);
  1419.  
  1420.  
  1421. # check whether the pixel's geographical nX, nY coordinates
  1422. # are within the normal range, then further crop to the angles
  1423. # given by the crop values.
  1424.  
  1425. L=nX-l0;
  1426. P=nY;
  1427. if( (L>=-pi && L<=pi && P>=-latcutoff && P<=+latcutoff),
  1428. # TAKE the value of X,Y FROM the original equirectangular source map
  1429. # and input the value into x,y in the output.
  1430. # X%ow, when combined with an l0 longitude shift, allows for
  1431. # rotating around the globe, otherwise it just cuts off at 180.
  1432. result=I(#0,X%ow,Y,0,2,0),
  1433. # else output transparent pixel
  1434. result=[0,0,0,0]
  1435. );
  1436.  
  1437. result"
  1438.  
  1439. autocrop[-1]
  1440.  
  1441.  
  1442. # Keep the last (. = [-1] image index) image generated
  1443. -keep[-1]
  1444.  
  1445.  
  1446. #######################
  1447. ## ALBERS EQUAL AREA CONIC
  1448. #@gmic alberts_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1449. #@gmic : Take a map or equirectangular panorama
  1450. #@gmic : and rotate it around z, y and x angles,
  1451. #@gmic : then project it as an Albers equal area conic map.
  1452.  
  1453. #@gui Albers Equal Area Conic : albers_projection
  1454. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an Albers equal area conic projection.")
  1455. #@gui : Roll = float(0,-360,360)
  1456. #@gui : Pitch = float(0,-360,360)
  1457. #@gui : Yaw = float(0,-360,360)
  1458. #@gui : Central longitude = float(0,-180,180)
  1459. #@gui : Reference parallel = float(0,-90,90)
  1460. #@gui : Standard parallel #1 = float(20,-90,90)
  1461. #@gui : Standard parallel #2 = float(60,-90,90)
  1462. #@gui : Crop latitude to = float(180,0,180)
  1463. #@gui : Crop longitude to = float(360,0,360)
  1464. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1465.  
  1466.  
  1467. albers_projection :
  1468.  
  1469. # $1-3 Default to 0 degrees
  1470. -skip ${1=0},${2=0},${3=0},${4=0},${5=0},${6=20},${7=60},${8=180},${9=360}
  1471.  
  1472. latcutoff={$8*pi/180}
  1473. longcutoff={$9*pi/180}
  1474.  
  1475. l0=$4
  1476. phi0=$5
  1477. phi1=$6
  1478. phi2=$7
  1479.  
  1480. # Give input alpha-channel to prevent errors with transparency
  1481. to_a[0]
  1482.  
  1483. _fF_zoomaccurate=1
  1484. _fF_zoomfactor=1
  1485.  
  1486. # First rotate the map with the appropriate custom command
  1487. if $1!=0" || "$2!=0" || "$3!=0
  1488. rotate_equirectangular_map $1,$2,$3
  1489. endif
  1490.  
  1491. # width and height of original map
  1492. nw={0,w}
  1493. nh={0,h}
  1494.  
  1495. # Initialize the image using the above precalculated width and height
  1496. # The "" block contains the mathematical 'formula' that determines
  1497. # what each x,y pixel is.
  1498.  
  1499. -input {$nw+1},{$nw+1},1,4,"
  1500.  
  1501. #inside the ""-block it's easier to pre-translate dollar-variables
  1502. pi2=2*pi;
  1503. phi0="$phi0"*pi/180;
  1504. l0="$l0"*pi/180;
  1505.  
  1506. phi1="$phi1"*pi/180;
  1507. phi2="$phi2"*pi/180;
  1508.  
  1509. latcutoff="$latcutoff";
  1510. longcutoff="$longcutoff";
  1511.  
  1512.  
  1513.  
  1514.  
  1515. # Use centered coordinates. 0,0 is in the middle of the
  1516. # projection
  1517. cx=x-w/2+1;
  1518. cy=y-h/2+1;
  1519.  
  1520.  
  1521. # Turn them into radians
  1522. cx = (cx/w *pi2*8/7);
  1523. cy = -(cy/w *pi2*8/7);
  1524.  
  1525.  
  1526. # Formulas and constants:
  1527. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1528.  
  1529. if((phi1==phi2),
  1530. n=sin(phi1),
  1531. n= (sin(phi1)+sin(phi2))/2
  1532. );
  1533.  
  1534. C= cos(phi1)^2 + 2*n*sin(phi1);
  1535.  
  1536. rho0=sqrt((C-2*n*sin(phi0))) / n;
  1537.  
  1538. rho=sqrt(cx^2 + (rho0-cy)^2);
  1539.  
  1540. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1541. # page 114, if using ATAN2, the signs must be reversed.
  1542.  
  1543. if( (n>0),
  1544. (theta=atan2(cx, (rho0-cy))),
  1545. (theta=atan2(-cx, (-rho0+cy)))
  1546. );
  1547.  
  1548. # When the reference parallel is in the southern hemisphere,
  1549. # flip the projection and map
  1550. if( (phi0<0),
  1551. (n=-n)
  1552. );
  1553.  
  1554.  
  1555. if( (rho==0),
  1556. (nY = asin( (C-(rho*n)^2) / (2*n)) ),
  1557. (nY = asin( (C-(rho*n)^2) / (2*n)) )
  1558. );
  1559.  
  1560. nX = l0 + theta/n;
  1561.  
  1562. # For the logic why I'm using inverse formulas, check out the
  1563. # Cylindrical Equal-Area filter
  1564.  
  1565. # EQUIRECTANGULAR COORDINATES
  1566. # Get the width and height of the original equirectangular input
  1567. ow="{0,w}";
  1568. oh="{0,h}";
  1569.  
  1570. # Pretty simply translate radians coordinates straight into
  1571. # equirectangular pixel coordinates.
  1572. X=ow/2 + (ow * nX/pi2);
  1573. Y=oh/2 + (oh * -nY/pi);
  1574.  
  1575.  
  1576. # check whether the pixel's geographical nX, nY coordinates
  1577. # are within the normal range, then further crop to the angles
  1578. # given by the crop values.
  1579.  
  1580. L=nX;
  1581. P=nY;
  1582. if( (L>=l0-longcutoff/2 && L<=l0+longcutoff/2 && P>=phi0-latcutoff/2 && P<=phi0+latcutoff/2),
  1583. # TAKE the value of X,Y FROM the original equirectangular source map
  1584. # and input the value into x,y in the output.
  1585. result=I(#0,X%ow,Y,0,2,0),
  1586. # else output transparent pixel
  1587. result=[0,0,0,0]
  1588. );
  1589.  
  1590. result"
  1591.  
  1592. # If the reference parallel is in the southern hemisphere,
  1593. # we've flipped the projection. Here we need to further rotate
  1594. # the image so it's not upside down.
  1595. if $phi0<0
  1596. rotate[-1] 180
  1597. endif
  1598.  
  1599.  
  1600. # New trick: if you can't calculate the width and height of the image
  1601. # in advance, doing it on a "blank" square image can work as well.
  1602. # In the ""-block of the -input above the last pixel assignment
  1603. # checks if the geographical coordinates we're looking for make
  1604. # any sense, that is, no weird values above or below plus-minus 180 or
  1605. # 90 degrees. If the pixel x,y coordinates would result in something
  1606. # nonsensical, output a transparent pixel. This is necessary, because it's
  1607. # not automatic: you get weird blocks of white *and* transparent
  1608. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1609. # the transparent areas.
  1610. autocrop[-1]
  1611.  
  1612.  
  1613. # Keep the last (. = [-1] image index) image generated
  1614. -keep[-1]
  1615.  
  1616.  
  1617. #######################
  1618. ## ECKERT IV
  1619. #@gmic eckert_iv_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1620. #@gmic : Take a map or equirectangular panorama
  1621. #@gmic : and rotate it around z, y and x angles,
  1622. #@gmic : then project it as an Eckert IV map.
  1623.  
  1624. #@gui Eckert IV : eckert_iv_projection
  1625. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an Eckert IV projection.")
  1626. #@gui : Roll = float(0,-360,360)
  1627. #@gui : Pitch = float(0,-360,360)
  1628. #@gui : Yaw = float(0,-360,360)
  1629. #@gui : Central longitude = float(0,-180,180)
  1630. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1631.  
  1632.  
  1633. eckert_iv_projection :
  1634.  
  1635. # $1-3 Default to 0 degrees
  1636. -skip ${1=0},${2=0},${3=0},${4=0}
  1637.  
  1638. l0=$4
  1639.  
  1640. # Give input alpha-channel to prevent errors with transparency
  1641. to_a[0]
  1642.  
  1643. _fF_zoomaccurate=1
  1644. _fF_zoomfactor=1
  1645.  
  1646. # First rotate the map with the appropriate custom command
  1647. if $1!=0" || "$2!=0" || "$3!=0
  1648. rotate_equirectangular_map $1,$2,$3
  1649. endif
  1650.  
  1651. # width and height of original map
  1652. nw={0,w}
  1653. nh={0,h}
  1654.  
  1655. # Initialize the image using the above precalculated width and height
  1656. # The "" block contains the mathematical 'formula' that determines
  1657. # what each x,y pixel is.
  1658.  
  1659. -input {$nw+1},{$nw+1},1,4,"
  1660.  
  1661. #inside the ""-block it's easier to pre-translate dollar-variables
  1662. pi2=2*pi;
  1663.  
  1664. l0="$l0"*pi/180;
  1665.  
  1666.  
  1667. # Use centered coordinates. 0,0 is in the middle of the
  1668. # projection
  1669. cx=x-w/2+1;
  1670. cy=y-h/2+1;
  1671.  
  1672.  
  1673. # Turn them into radians
  1674. cx = (cx/w *pi*sqrt(3));
  1675. cy = -(cy/w *pi*sqrt(3));
  1676.  
  1677.  
  1678. # Formulas and constants:
  1679. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1680.  
  1681. theta=asin(cy*sqrt(4+pi) / (2*sqrt(pi)));
  1682.  
  1683. # When the reference parallel is in the southern hemisphere,
  1684. # flip the projection and map
  1685. nY=asin( (theta + sin(theta)*cos(theta) + 2*sin(theta)) / (2+pi/2));
  1686. nX=l0 + sqrt(pi*(4+pi)) * cx/(2*(1+cos(theta)));
  1687.  
  1688. # For the logic why I'm using inverse formulas, check out the
  1689. # Cylindrical Equal-Area filter
  1690.  
  1691. # EQUIRECTANGULAR COORDINATES
  1692. # Get the width and height of the original equirectangular input
  1693. ow="{0,w}";
  1694. oh="{0,h}";
  1695.  
  1696. # Pretty simply translate radians coordinates straight into
  1697. # equirectangular pixel coordinates.
  1698. X=ow/2 + (ow * nX/pi2);
  1699. Y=oh/2 + (oh * -nY/pi);
  1700.  
  1701.  
  1702. # check whether the pixel's geographical nX, nY coordinates
  1703. # are within the normal range, then further crop to the angles
  1704. # given by the crop values.
  1705.  
  1706. L=nX;
  1707. P=nY;
  1708. if( (L>=l0-pi && L<=l0+pi && P>=-pi/2 && P<=pi/2),
  1709. # TAKE the value of X,Y FROM the original equirectangular source map
  1710. # and input the value into x,y in the output.
  1711. result=I(#0,X%ow,Y,0,2,0),
  1712. # else output transparent pixel
  1713. result=[0,0,0,0]
  1714. );
  1715.  
  1716. result"
  1717.  
  1718. # If the reference parallel is in the southern hemisphere,
  1719. # we've flipped the projection. Here we need to further rotate
  1720. # the image so it's not upside down.
  1721. if $phi0<0
  1722. rotate[-1] 180
  1723. endif
  1724.  
  1725.  
  1726. # New trick: if you can't calculate the width and height of the image
  1727. # in advance, doing it on a "blank" square image can work as well.
  1728. # In the ""-block of the -input above the last pixel assignment
  1729. # checks if the geographical coordinates we're looking for make
  1730. # any sense, that is, no weird values above or below plus-minus 180 or
  1731. # 90 degrees. If the pixel x,y coordinates would result in something
  1732. # nonsensical, output a transparent pixel. This is necessary, because it's
  1733. # not automatic: you get weird blocks of white *and* transparent
  1734. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1735. # the transparent areas.
  1736. autocrop[-1]
  1737.  
  1738.  
  1739. # Keep the last (. = [-1] image index) image generated
  1740. -keep[-1]
  1741.  
  1742. #######################
  1743. ## ECKERT VI
  1744. #@gmic eckert_vi_projection : yaw angle, pitch angle, roll angle, central_longitude angle, reference_parallel angle, standard_parallel1 angle, standard_parallel2 angle, latcutoff angle, longcutoff angle
  1745. #@gmic : Take a map or equirectangular panorama
  1746. #@gmic : and rotate it around z, y and x angles,
  1747. #@gmic : then project it as an Eckert VI map.
  1748.  
  1749. #@gui Eckert VI : eckert_vi_projection
  1750. #@gui : note = note("Rotate an equirectangular map or panorama around three axises in order and output it as an Eckert VI projection.")
  1751. #@gui : Roll = float(0,-360,360)
  1752. #@gui : Pitch = float(0,-360,360)
  1753. #@gui : Yaw = float(0,-360,360)
  1754. #@gui : Central longitude = float(0,-180,180)
  1755. #@gui : note = note("<small>Author: <i>Kristian Järventaus</i>.      Latest Update: <i>2020-08-31</i>.</small>")
  1756.  
  1757.  
  1758. eckert_vi_projection :
  1759.  
  1760. # $1-3 Default to 0 degrees
  1761. -skip ${1=0},${2=0},${3=0},${4=0}
  1762.  
  1763. l0=$4
  1764.  
  1765. # Give input alpha-channel to prevent errors with transparency
  1766. to_a[0]
  1767.  
  1768. _fF_zoomaccurate=1
  1769. _fF_zoomfactor=1
  1770.  
  1771. # First rotate the map with the appropriate custom command
  1772. if $1!=0" || "$2!=0" || "$3!=0
  1773. rotate_equirectangular_map $1,$2,$3
  1774. endif
  1775.  
  1776. # width and height of original map
  1777. nw={0,w}
  1778. nh={0,h}
  1779.  
  1780. # Initialize the image using the above precalculated width and height
  1781. # The "" block contains the mathematical 'formula' that determines
  1782. # what each x,y pixel is.
  1783.  
  1784. -input {$nw+1},{$nw+1},1,4,"
  1785.  
  1786. #inside the ""-block it's easier to pre-translate dollar-variables
  1787. pi2=2*pi;
  1788.  
  1789. l0="$l0"*pi/180;
  1790.  
  1791.  
  1792. # Use centered coordinates. 0,0 is in the middle of the
  1793. # projection
  1794. cx=x-w/2+1;
  1795. cy=y-h/2+1;
  1796.  
  1797.  
  1798. # Turn them into radians
  1799. cx = (cx/w *pi*2);
  1800. cy = -(cy/w *pi*2);
  1801.  
  1802.  
  1803. # Formulas and constants:
  1804. # https://epic.awi.de/id/eprint/39585/1/USGS_Bulletin_1532.pdf
  1805.  
  1806. theta=sqrt(2+pi)*cy/2;
  1807.  
  1808. # When the reference parallel is in the southern hemisphere,
  1809. # flip the projection and map
  1810. nY=asin( (theta + sin(theta)) / (1+pi/2));
  1811. nX=l0 + sqrt(2+pi) * cx/(1+cos(theta));
  1812.  
  1813. # For the logic why I'm using inverse formulas, check out the
  1814. # Cylindrical Equal-Area filter
  1815.  
  1816. # EQUIRECTANGULAR COORDINATES
  1817. # Get the width and height of the original equirectangular input
  1818. ow="{0,w}";
  1819. oh="{0,h}";
  1820.  
  1821. # Pretty simply translate radians coordinates straight into
  1822. # equirectangular pixel coordinates.
  1823. X=ow/2 + (ow * nX/pi2);
  1824. Y=oh/2 + (oh * -nY/pi);
  1825.  
  1826.  
  1827. # check whether the pixel's geographical nX, nY coordinates
  1828. # are within the normal range, then further crop to the angles
  1829. # given by the crop values.
  1830.  
  1831. L=nX;
  1832. P=nY;
  1833. if( (L>=l0-pi && L<=l0+pi && P>=-pi/2 && P<=pi/2),
  1834. # TAKE the value of X,Y FROM the original equirectangular source map
  1835. # and input the value into x,y in the output.
  1836. result=I(#0,X%ow,Y,0,2,0),
  1837. # else output transparent pixel
  1838. result=[0,0,0,0]
  1839. );
  1840.  
  1841. result"
  1842.  
  1843. # If the reference parallel is in the southern hemisphere,
  1844. # we've flipped the projection. Here we need to further rotate
  1845. # the image so it's not upside down.
  1846. if $phi0<0
  1847. rotate[-1] 180
  1848. endif
  1849.  
  1850.  
  1851. # New trick: if you can't calculate the width and height of the image
  1852. # in advance, doing it on a "blank" square image can work as well.
  1853. # In the ""-block of the -input above the last pixel assignment
  1854. # checks if the geographical coordinates we're looking for make
  1855. # any sense, that is, no weird values above or below plus-minus 180 or
  1856. # 90 degrees. If the pixel x,y coordinates would result in something
  1857. # nonsensical, output a transparent pixel. This is necessary, because it's
  1858. # not automatic: you get weird blocks of white *and* transparent
  1859. # lines. Then you can use autocrop from g'mic to do some magic and remove
  1860. # the transparent areas.
  1861. autocrop[-1]
  1862.  
  1863.  
  1864. # Keep the last (. = [-1] image index) image generated
  1865. -keep[-1]
  1866.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!