transforms.ps

%%BeginResource: procset cass
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        This Postscript is Copyright (c) 2016, Peter J Billam          %
% Permission is granted  to any individual or institution to use, copy, %
% modify or redistribute this software, so long as it is not resold for %
% profit, and provided this notice is retained. It is provided "as is", %
% without any express or implied warranty.    http://www.pjb.com.au     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% http://www.math.ubc.ca/~cass/graphics/manual/  chapter 8
% ~/ps/cass/08_nonlinear_transformations.pdf

% 20160719 all these parameters ... it's ugly.
% My _xy2uv routines could take 3 arguments, the first of which
% is a dict containing the parameters, eg:
%  <<(ax 9  (wy) 300  (py) 0  (ay) 6  (wx) 200  (px) 90>>
%    /sine_xy2uv ctransform fill  pop
% This makes sine_xy2uv get each param every invocation, and needs a pop;
% so it might be faster to modify ctransform into c_transform
% which copies the params into a module-local dict and then pops.
% Could also react cleverly if there's no dict,
%   by re-using previous (or default) param-values, and not popping.
% Or, it could test which _xy2uv name it has been given,
%   and thence extract the relevant params in the right order
%   and pass them to the _xy2uv routine on the stack ?
%   But _xy2uv routine will still put them into its local dict...
% All these options are only implementation-details,
%   compared with a c_transform which just invokes {ctransform pop}
%   so could start there, though it still means rewriting all the _xy2uv's
%
% And/Or could handle an array,  containing ordered, unnamed params,
%   by leaving it on the stack for the _xy2uv to interpret? (then popping)
%
% OK :-) first, use the array; then many of the _xy2uv's could get directly,
%   which would be fast.  c_transform would check types, invoke, then pop.
% Then later, accept also named-params in a dict, as an alternative:
%   c_transform would check types, and construct the correct array
%   according to the name of the _xy2uv (or print helpful error messages)
%   The _xy2uv routines would not notice the change.
% Then in the API the params could be either named, or ordered.

% 20170621
% I also need somehow to do reflections, or other one-to-many transforms
% for example, a kaleidoscope needs to drop 5/6 of its x,y input
% and for the remaining 1/6 needs to return six u,v pairs !
% return an array either  [ ]  or  [ u1 v1 u2 v2 u3 v3 u4 v4 u5 v5 u6 v6 ]  ?
% PS offers rotate, but not reflect :-(     x0 y0 angle reflect
% should put together the necessary matrix then invoke transform !
% Better draw stroke and fill the whole page, clip to a 60-degree segment,
% then redraw :-( (eg: having saved any rangdom values) clip, reflect, rotate

/ctransform {   % the argument is the named function f(x,y) -> u,v
  load          % first load the procedure onto the stack
  1 dict begin  /f exch def  % and give it a local name
  % could keep track of x,y and weed out zero-length stuff ?
  % or break "long" line-segments into subsections, or into curves ?
  [  % build an array from the current path
    {   % x y
      [ 3 1 roll f {moveto} ]
    }
    {   % x y
      [ 3 1 roll f {lineto} ]
    }
    {   % x1 y1 x2 y2 x3 y3
      [ 7 1 roll    % [ P1 P2 P3
        f 6 2 roll  % [ U3 P1 P2
        f 6 2 roll  % [ U2 U3 P1
        f 6 2 roll  % [ U1 U2 U3
        { curveto }
      ]
    }
    { [ {closepath} ] }
    pathforall
  ]
  newpath  % and then replace the current path
  { aload pop exec } forall
  end
} def

/c_transform {
  exch /::p exch def   % should protect this in a local dictionary...
  ctransform    % first stub-prototype
} def

% ------------------   The *_xy2uv procedures ....  ----------------

/collapse_xy2uv {   % [ x0 y0 r0 range ]
  10 dict begin [ /y /x ] {exch def } forall
  /x0 ::p 0 get def  /y0    ::p 1 get def
  /r0 ::p 2 get def  /range ::p 3 get def
  /r  x x0 sub dup mul  y y0 sub dup mul  add  sqrt  def
  r r0 le {
    /u x0 def  /v y0 def
  } {
    r r0 range add ge {
      /u x def /v y def
    } {
      /theta  y y0 sub  x x0 sub  atan def
      /new_r  r r0 sub range div 0.35 exp  r mul  def
      /u  theta cos  new_r mul  x0 add  def
      /v  theta sin  new_r mul  y0 add  def
    } ifelse
  } ifelse
  u v
  end
} def

/perspective_xy2uv { 6 dict begin [ /y /x ] { exch def } forall
  % parameters:  [xvp yvp x0 theta]
  % xvp and yvp are the location of the vanishing-point.
  % x0 is the x at which the heights are left at 1:1 scale
  %   and at which x->u will be left unchanged.
  % theta is the angle at which the viewer sees the x0 point in the surface
  %   so if theta=90 the viewer is face-on to the x=u=x0 point
  %   this means solving du/dx = sin(u/x0) . u/x0     :-(
  % https://en.wikipedia.org/wiki/Lists_of_integrals#Lists_of_integrals
  % https://en.wikipedia.org/wiki/List_of_integrals_of_trigonometric_functions
  ::p length 3 eq {  % no theta is present
    /r 2.718281828  x ::p 0 get sub ::p 2 get ::p 0 get sub div 1 sub  exp def
  } {
    /sintheta ::p 3 get sin def % amateurish kludge: use a constant angle :-(
    /r 2.718281828  x ::p 0 get sub ::p 2 get ::p 0 get sub div
      1 sub sintheta mul  exp def
  } ifelse
  /u ::p 2 get ::p 0 get sub r mul  ::p 0 get add  def
  /v  y ::p 1 get sub  r mul  ::p 1 get add  def
  u v
end } def

/perspective_flat_xy2uv { 6 dict begin [/y /x] {exch def} forall
  % parameters:  [xvp yvp y0 theta]
  % as perspective_xy2uv, except that the surface is flat rather than vertical
  ::p length 3 eq {  % no theta is present
    /r 2.718281828  y ::p 1 get sub ::p 2 get ::p 1 get sub div 1 sub  exp def
  } {
    /sintheta ::p 3 get sin def % kludge: use a constant angle :-(
    /r 2.718281828  y ::p 1 get sub ::p 2 get ::p 1 get sub div
      1 sub sintheta mul  exp  def
  } ifelse
  /u  x ::p 0 get sub  r mul  ::p 0 get add  def
  /v ::p 2 get ::p 1 get sub r mul  ::p 1 get add def
  u v
end } def

/sine_xy2uv {   % needs the parameters:
  % [x_amplitude y_wavelength y_phase  y_amplitude x_wavelength x_phase]
  % The wavelength is in the x dimension if the amplitude is in the y dimension
  % u = x + x_amplitude.sin(360y/y_wavelength + y_phase)
  % v = y + y_amplitude.sin(360x/x_wavelength + x_phase)
  3 dict begin  [ /y /x ] { exch def } forall
  360 y mul ::p 1 get div  ::p 2 get add  sin  ::p 0 get mul  x add
  360 x mul ::p 4 get div  ::p 5 get add  sin  ::p 3 get mul  y add
  end
} def

/hardlens_xy2uv {   % params: [ xcentre ycentre r magnification ]
  10 dict begin  [ /y /x ] { exch def } forall
  /x0 ::p 0 get def  /y0  ::p 1 get def
  /r0 ::p 2 get def  /mag ::p 3 get def
  /r  x x0 sub dup mul  y y0 sub dup mul  add  sqrt  def
  r r0 ge {
    /u x def  /v y def
  } {
    /stretch  r r0 div 90 mul cos  mag 1 sub  mul  1 add  def
    % problem: if mag>2 then some points inside get moved to outside :-(
    r r0 div stretch mul 1.0 gt { /stretch r0 r div def } if   % kludge
    /u  x x0 sub stretch mul  x0 add  def
    /v  y y0 sub stretch mul  y0 add  def
  } ifelse
  u v
  end
} def

/softlens_xy2uv {   % params: [ xcentre ycentre r magnification ]
  10 dict begin  [ /y /x ] { exch def } forall
  /x0 ::p 0 get def  /y0  ::p 1 get def
  /r0 ::p 2 get def  /mag ::p 3 get def
  /r  x x0 sub dup mul  y y0 sub dup mul  add  sqrt  def
  r r0 ge {
    /u x def  /v y def
  } {
    /stretch  r r0 div 180 mul cos  1 add  0.5 mul mag 1 sub mul  1 add  def
    % problem: if mag>2 then some points inside get moved to outside :-(
    r r0 div stretch mul 1.0 gt { /stretch r0 r div def } if   % kludge
    /u  x x0 sub stretch mul  x0 add  def
    /v  y y0 sub stretch mul  y0 add  def
  } ifelse
  u v
  end
} def

/twist_xy2uv {   % params: [ xcentre ycentre r angle ]
  10 dict begin  [ /y /x ] { exch def } forall
  /x0 ::p 0 get def  /y0  ::p 1 get def
  /r0 ::p 2 get def  /ang ::p 3 get def
  /r  x x0 sub dup mul  y y0 sub dup mul  add  sqrt  def
  r 0 eq {
    x y
  } {
    /theta  y y0 sub  x x0 sub  atan def
    /newtheta 2.718281828  0 r r0 div sub  exp  ang  mul theta add  def
    /u  r newtheta cos mul  x0 add  def
    /v  r newtheta sin mul  y0 add  def
    u v
  } ifelse
  end
} def

/grand_xy2uv {  % usage:  [stdev_x stdev_y]  /grand_xy2uv  c_transform
  10 dict begin  [ /y /x ] { exch def } forall
  /stdev_x ::p 0 get def  /stdev_y ::p 1 get def
  /u x stdev_x grand def  /v y stdev_y grand def
  u v
  end
} def

%%EndResource