(* Copyright (c) Dror Bar-Natan 1991 *)

BeginPackage["PDI`"]

PDIPlot::usage ="\n
  PDIPlot[exp,{x,xmin,xmax},{y,ymin,ymax}] returns a perceived-depth-image\n
  of the graph of the function given by exp. Options:\n
    PlotPoints (can be an integer or pair of integers; defaults to 60);\n
    Periods (6); PlotRange ({-1,1}); Density (.5); ApparentDepth (1.);\n
    Guides (True); BasicBlock (Rectangle[#1-#2,#1+#2]&).\n
    Other options are passed on to Graphics."

SetAttributes[PDIPlot,HoldFirst];

Options[PDIPlot] := {PlotPoints -> 60, Periods -> 6, Guides -> True,
  PlotRange -> {-1,1}, Density -> .5, ApparentDepth -> 1.,
  BasicBlock :> (Rectangle[#1-#2,#1+#2]&)}

Begin["`private`"]

(* A new Graphics primitive *)
Unprotect[Display]
Display[channel_,graphics_?(!FreeQ[#,PDIArray]&)] :=
  (Display[channel,graphics /. 
    (PDIArray[basicblock_,size_,pts_] :> 
      (basicblock[#,size]& /@ 
        (Flatten[Outer[List,First[#],{Last[#]}],1])& /@ pts))];
  graphics)
Protect[Display]

PDIPlot[expr_,{x_Symbol,xmin_?NumberQ,xmax_?NumberQ},
	      {y_Symbol,ymin_?NumberQ,ymax_?NumberQ},opts___] :=
  Block[{plotpts    =(PlotPoints     /. {opts} /. Options[PDIPlot]),
         periods    =(Periods        /. {opts} /. Options[PDIPlot]),
         zrange     =(PlotRange      /. {opts} /. Options[PDIPlot]),
         density    =(Density        /. {opts} /. Options[PDIPlot]),
         depth      =(ApparentDepth  /. {opts} /. Options[PDIPlot]),
         basicblock =(BasicBlock     /. {opts} /. Options[PDIPlot]),
         guides     =(Guides         /. {opts} /. Options[PDIPlot]),
         xres,xpts,xv,xsample,xsteps,xfull,xshift,
         ypts,i,yv,zmin,zmax,zmid,zscale},

(* Some scaling calculations *)
    {zmin,zmax} = zrange ; xres=Floor[xpts/periods];
    zmid = (zmax+zmin)/2 ; zscale = .06depth(xmax-xmin)/(zmax-zmin);
    {xpts,ypts} = If[Length[plotpts]==2,plotpts,{plotpts,plotpts}];
(* The strips *)
    xsteps=xres*Table[i,{i,-.5,periods-1.5}];
(* The main program *)
    Show[Graphics[{PDIArray[basicblock,{(xmax-xmin)/xpts,(ymax-ymin)/ypts}/2,
(* Loop over the values of yv *)
      Table[
(* Filling the current row of the left-most strip randomly *)
        xsample=Flatten[Position[Table[Random[]<density,{xres}],True]];
(* The x values to be used in this row *)
        xv=N[xmin+Outer[Plus,xsteps,xsample-1]*(xmax-xmin)/xpts];
(* Computing the left/right shift of every block *)
        z=N[Map[(expr /. {y->yv,x->#})&,xv,{2}]]-zmid;
        xshift=zscale Accumulate[Plus,z];
(* Final positioning of the blocks *)
        xfull=xv-((#-Last[xshift]/2)& /@ xshift);
        {Flatten[xfull],yv}, {yv,ymin,ymax,(ymax-ymin)/ypts}]],
(* If guides==True, add guiding rectangles *)
      If[guides, (Rectangle @@ #)& /@ Map[
          ({.5xmax+.5xmin,1.1ymin-.1ymax}+{xmax-xmin,ymax-ymin}#/2/periods)&,
          {{{-2.1,-.1},{-1.9,.1}},{{-.1,-.1},{.1,.1}}},{2}], {}]},
      Sequence@@Select[{opts},!MemberQ[First /@ Options[PDIPlot],First[#]]&]
  ]] ]

End[]
EndPackage[]