Making the Most of the Visualization Functions

Ulises Cervantes-Pimentel
Brett Champion

Wolfram Research, Inc.

Visualization in Mathematica

Visualization in Mathematica reaches a new level in V6.0. All the plotting operations, including function and data oriented operations, have been reimplemented from the ground up to take advantage of areas such as symbolics, numerics, geometric constructions, adaptive meshes, and spatial data structures. We have integrated a large number of pre- and post-processing operations such as function analysis, region selection, line overlays, etc., as an integral part of the geometry construction. This approach allows us to generate high quality graphics which were not possible before, not only in Mathematica, but in any other visualization system. With this technology, we put in the hands of our users the power to generate amazing visualizations as well as new tools for scientific explorations.

Function and Data Based

Function Evaluation Based

Plot

ParametricPlot

LogLinearPlot

LogLogPlot

Plot3D

ParametricPlot3D

ContourPlot

DensityPlot

ContourPlot3D

RegionPlot

RegionPlot3D

SphericalPlot3D

RevolutionPlo3D

Data Based

ListPlot

ListLinePlot

ListLogPlot

ListPlot3D

ListContourPlot

ListDensityPlot

ListContourPlot3D

ListSurfacePlot3D

ListPointPlot3D

ArrayPlot

ReliefPlot

Visualization Framework

For all the visualization functions, we can identify three major areas of optimization: geometric, numeric and symbolic.

Geometric

Adaptive Refinement

Mesh lines and points overlay

MeshFunctions

MeshShading

Arbitrary clipping by using a RegionFunction

Irregularities and Singularities removal

Extended Filling between Curves and Surfaces in Plot and Plot3D

Improved ContourPlot labeling

Irregular data plotting

Numeric

Arbitrary precision via a WorkingPrecision option

Evaluation monitor support

FiniteDifference and Symbolic vertex normals computation

Symbolic

Vector valued function plotting

Piecewise functions handling

Exclusions to handle Jumps and Singularities

Improved ContourPlot3D

Implicit Algebraic Surfaces

Adaptive Refinement

One of the basic changes for function based visualization functions is the use of adaptive refinement. This is a major improvement over the previous versions. It is possible to reduce the number of function evaluations by using lower values of PlotPoints, and still get smooth graphs and surfaces.

Plot[Sin[1/x],{x,0,1},Mesh->All]

ParametricPlot3D[{Sin[u],Sin[v],Sin[u+v]},{u,0,2 π},{v,0,2 π},Mesh->All]

This is another example, in this case for ParametricPlot3D, that demonstrates how the adaptivity helps to reduce the number of function evaluations.

Clear[r,KleinBottle] 
r[u_]:=4 (1-((Cos[u])/2)); 

KleinBottle[u_/;u≤π,v_]:={6Cos[u](1+Sin[u])+r[u]Cos[u]Cos[v],16Sin[u]+r[u]Sin[u]Cos[v],r[u]Sin[v]}; KleinBottle[u_/;u>π,v_]:={6 Cos[u] (1+Sin[u])+r[u] Cos[v+π],16 Sin[u],r[u] Sin[v]};

kbtable=Table[ParametricPlot3D[KleinBottle[u,v],{u,0,2π},{v,0,2π},ViewPoint->{-2,-8,-15},PlotPoints->5, MaxRecursion->rr,PlotStyle->Yellow,ImageSize->{220,Automatic},Evaluate@SemStyle1],{rr,0,3}];

Grid[Partition[kbtable,2]]

Computational Mesh and Visual Mesh

Manipulate[Block[{$PerformanceGoal="Quality",vpt=ViewPoint-> {0,0,Infinity},p1,p2},{p1,p2}={Plot3D[Sin[x + Cos[y]],{x,0,6},{y,-3,3},Mesh-> All,PlotStyle-> None,MaxRecursion-> n],Plot3D[Sin[x + Cos[y]],{x,0,6},{y,-3,3},PlotStyle-> None,MaxRecursion-> n]};GraphicsGrid[{{Show[p1,vpt],Show[p2,vpt]},{p1,p2}},ImageSize-> 450]],{n,Range[0,4]}]

Analytic Region Plotting

RegionPlot is used to visualize regions defined by boolean functions. Adaptivity is used consistently as with all the other visualization functions.

pstyle = {Mesh->All,PlotStyle->Blue};

tt1=Table[RegionPlot[Mod[Ceiling[x]+Ceiling[y],2]==0,{x,-3,3},{y,-3,3},MaxRecursion->rr,Frame->True,Evaluate@Sequence[pstyle]],{rr,0,3}];

tt2 = Table[RegionPlot[ (((x+(1/2)))^2)+(y^2)≤1⊻(((x-(1/2)))^2)+(y^2)≤1, {x,-1.61,1.6}, {y,-1.015,1},PlotPoints->9,MaxRecursion->rr,Evaluate@Sequence[pstyle]],{rr,0,3}];

tt3 = Table[RegionPlot[Im[ChebyshevT[12,x+I y]]>0,{x,-3,3},{y,-3,3},MaxRecursion->rr,PlotPoints->26,Evaluate@Sequence[pstyle]],{rr,0,3}];

GraphicsGrid[{tt1,tt2,tt3},ImageSize->{750,Automatic}]

Adaptive Contouring

For ContourPlot, the adaptivity is based on the contour lines in addition to the function itself.

fun[x_,y_]:=Arg[Tanh[(1/(x+I y))]];
t1=Table[ContourPlot[fun[x,y],{x,-(1/4),(1/4)},{y,-(1/2),(1/2)},MaxRecursion->rr,PlotPoints->15,ColorFunction->"BeachColors",Mesh->All,Exclusions->None,Contours->5],{rr,0,1}];
t2={ContourPlot[fun[x,y],{x,-(1/4),(1/4)},{y,-(1/2),(1/2)},MaxRecursion->3,ColorFunction->"BeachColors",Mesh->All,Exclusions->None,Contours->5],
ContourPlot[fun[x,y],{x,-(1/4),(1/4)},{y,-(1/2),(1/2)},MaxRecursion->3,ColorFunction->"BeachColors",MeshFunctions->{#1&,#2&},Mesh->6,Exclusions->None,Contours->5]};
GraphicsGrid[{t1,t2},ImageSize->{550,Automatic}]

Implicit Surfaces Based on an Adaptive Octree

For ContourPlot3D, we use adaptivity based on cubic cells.

Clear[x,y,z,pp1]
pp1=16 (x^3)-48 y (x^2)+24 z (x^2)-48 (y^2) x+16 (y^3)-31 (z^3)-93.5307 (z^2)+24 (y^2) z-72 z;

pstyle = {ImageSize->250,Boxed->False,Axes->None};

Grid[{Table[ContourPlot3D[pp1==0,{x,-2,2},{y,-2,2},{z,-2.3,2},Method->{"Refinement"->{Method->"OctreeCubes"}},PlotPoints->5,MaxRecursion->r,Evaluate@Sequence[pstyle]],{r,0,2}]}]

Grid[
{ContourPlot3D[pp1==0,{x,-2,2},{y,-2,2},{z,-2.3,2},Evaluate@Sequence[pstyle],Mesh->#]& /@ {All,None,8}}]

Performance Goal

Manipulate makes extensive use of an automatic feature of the plotting functions, that is, PerformanceGoal. During the dragging of sliders, Manipulate uses a "Speed" mode which will lower the quality requirements of the plot. This allows for faster refresh at the expense of generating low quality plots. It is possible to force a "Quality" mode or even a "Speed" mode explicitly.

Manipulate[
	ParametricPlot3D[{Cos[t]*(3 + aa*Cos[u]), Sin[t]*(3 + aa*Cos[u]), bb*Sin[u]}, {t, 0, cc}, {u, 0, 2*Pi}, Mesh -> mm, PerformanceGoal->pg,PlotRange -> {{-5.5, 5.5}, {-5.5, 5.5}, {-2.5, 2.5}}],

	{{mm,None,"Mesh: "},{None,All}},
	{{pg,Automatic,"PerformanceGoal: "},{Automatic,"Speed","Quality"}},
	{{aa, 15/8}, 1/4, 5/2,1/16,Appearance-> "Labeled"},
	{{bb, 3/4}, 1/4, 5/2,1/16,Appearance-> "Labeled"},
	{{cc, 3*(Pi/2)}, Pi, 2*Pi}
]

Arbitrary Mesh Functions over harmonic Functions

The use of MeshFunctions help to understand properties of the functions.

u[x_,y_]:=((2 ArcTan[((2 x y (Sqrt[(((x^2)+(y^2))/(y^2))]))/((1-((x^2)+(y^2))) (Sqrt[(x^2)+(y^2)])))])/π);

Plot3D[u[x,y],{x,-3,3},{y,-3,3},Exclusions->{x y==0,(x^2)+(y^2)==1},MeshFunctions->{(#1^2+#2^2)&,ArcTan[#1/#2]&},Mesh->{8,12},MeshStyle->{{Red,Opacity[.3]},{Blue,Opacity[.3]}},PlotPoints->40,PlotStyle->{{Yellow,Specularity[White,50]}},ViewPoint->{1,1,1.5},NormalsFunction->"Weighted",ImageSize->{Automatic,480},Boxed->False,Axes->False]

This example constructs the MeshFunction as the Gaussian curvature of a given function.

Clear[tt,x,y,z,t,K]

t=7.5;
tt=(x^2) Cos[x]+(y^2) Cos[y]+(z^2) Cos[z]+((37 Cos[(t/5)]+33)/2.);
K[x_, y_, z_] = -((1/((∂_x)tt (∂_x)tt+(∂_y)tt (∂_y)tt+(∂_z)tt (∂_z)tt)^2)*Simplify[(-(∂_z)tt)(∂_z)tt((∂_(x,x))tt (∂_(y,y))tt-(∂_(x,y))tt (∂_(x,y))tt)-(∂_y)tt (∂_y)tt ((∂_(x,x))tt (∂_(z,z))tt-(∂_(x,z))tt (∂_(x,z))tt)+2 ((∂_y)tt (∂_z)tt ((∂_(y,z))tt (∂_(x,x))tt-(∂_(x,y))tt (∂_(x,z))tt)+(∂_x)tt (∂_z)tt ((∂_(x,z))tt (∂_(y,y))tt-(∂_(x,y))tt (∂_(y,z))tt)+(∂_x)tt (∂_y)tt ((∂_(x,y))tt (∂_(z,z))tt-(∂_(x,z))tt (∂_(y,z))tt))-(∂_x)tt (∂_x)tt ((∂_(y,y))tt (∂_(z,z))tt-(∂_(y,z))tt (∂_(y,z))tt)])

ContourPlot3D[tt,{x,-5,5},{y,-5,5},{z,-5,5},PlotPoints->40,Contours->{0.0},Mesh->6,MeshStyle->Black,MeshFunctions->(K[#1,#2,#3]&),MeshShading->Table[{RGBColor[1,t,0]},{t,0,1,0.17}],MaxRecursion->0,ImageSize->{Automatic,480},Boxed->False,Axes->None]

Clear[r,KleinBottle]

r[u_]:=4 (1-((Cos[u])/2));
KleinBottle[u_/;u≤π,v_]:={6 Cos[u] (1+Sin[u])+r[u] Cos[u] Cos[v],16 Sin[u]+r[u] Sin[u] Cos[v],r[u] Sin[v]}
KleinBottle[u_/;u>π,v_]:={6 Cos[u] (1+Sin[u])+r[u] Cos[v+π],16 Sin[u],r[u] Sin[v]}

ParametricPlot3D[KleinBottle[u,v],{u,0,2 Pi},{v,0,2 Pi},ViewPoint->{-2,-8,-15},Mesh->{55,30},ImageSize->{Automatic,480},MeshShading->{{Orange,None,Blue,None},{Green,Pink,Yellow,Lighter[Brown]}},Evaluate@SemStyle1]

Working Precision

By default all functions evaluations are performed using machine precision. This will produce an incorrect plot for large numbers.

Plot[Sin[2^50+x],{x,0,Pi},Mesh->All]

This shows how to specify a given precision by using the WorkingPrecision option.

Plot[Sin[2^50+x],{x,0,Pi},WorkingPrecision->10,Mesh->All]

Irregularities Removal

Not real valued. The visualization functions in Mathematica are defined to visualize real valued functions only. Hence if a function becomes complex (or other) valued at some point we exclude those values. With the adaptive sampling scheme we will usually get fairly close to the real valued barrier.

There are many types of non-real values including: Complex, DirectedInfinity[c], Underflow[], Overflow[], Null etc. But literally anything can be returned.

Example of complex valued function.

Plot[(Sqrt[1-(x^2)]),{x,-2,2},PlotStyle->{{Blue,Thickness[0.01]}}]

Plot3D[(Sqrt[x y]),{x,-1,1},{y,-1,1},Mesh->All]

Plot3D[(Sqrt[64-(x^2)-(y^2)]),{x,-8,8},{y,-8,8},Mesh->All]

Automatic Processing of Piecewise Functions

CompareDensityPlot[(E^(Boole[(x^2)+(y^2) Boole[y>0]<1])), {x, -2, 2}, {y, -2, 2}]

CompareDensityPlot[Abs[x y] Boole[(x^2)+(y^2)≤(2^2)],{x,-5,5},{y,-5,5}]

CompareDensityPlot[(Sqrt[Abs[y-(x^2)]]) Boole[(x^2)≤y<4], {x, -5, 5}, {y, -2, 7},MaxRecursion->0]

Exclusions

Plot3D[Im[Sqrt[Cos[y]+I Sin[x]]],{x,-5,5},{y,-5,5},Mesh->None,ExclusionsStyle->Orange,Exclusions->None,PlotStyle->Orange,ImageSize->500]

Mathematica can automatically detect many cases of jumps and discontinuities.

Plot3D[Im[(Sqrt[Cos[y]+I Sin[x]])],{x,-5,5},{y,-5,5},Mesh->None,PlotStyle->Orange,ImageSize->500]

Exclusions as Level Sets

An exclusion is a region around a point or curve on which the function presents jumps.

Manipulate[
	DensityPlot[Im((Sqrt[x+I y])),{x,-1,1},{y,-1,1},
	ColorFunction->"BeachColors",ExclusionsStyle->{Red,None},Mesh->All,
	ImageSize->500,MaxRecursion->0,
	Exclusions->{Automatic,MaxRecursion->3,Offset->offset}],

	{{offset,0,"Offset"},{0.,0.05,0.1,0.15,0.2,0.25,0.3},ControlType->PopupMenu}
]

Manipulate[
	Plot3D[Im((Sqrt[x+I y])),{x,-1,1},{y,-1,1},
		ColorFunction->"BeachColors",ExclusionsStyle->{Red,None},Mesh->All,
		ImageSize->500,MaxRecursion->0,
		Exclusions->{Automatic,MaxRecursion->3,Offset->offset}],

	{{offset,0,"Offset"},{0.,0.05,0.1,0.15,0.2,0.25,0.3},ControlType->PopupMenu}
]

Manipulate[
	DensityPlot[Im((-4 (z^10)+4 (z^9)+4 (z^8)+6 (z^7)-7 (z^6)+7 (z^5)+5 (z^4)+(z^3)+6 (z^2)+4 z-5)^(1/7)/.z->x+I y),
		{x,-(3/2),3/2},{y,-(3/2),3/2},
		ColorFunction->"BeachColors",ExclusionsStyle->{Red,None},
		Evaluated->True,Mesh->All,ImageSize->500,MaxRecursion->1,
		Exclusions->{Automatic,MaxRecursion->3,Offset->offset,TimeConstraint->10}],

	{{offset,0,"Offset"},{0.,0.02,0.04,0.06,0.1,0.15,0.2},ControlType->PopupMenu}
]

Plot3D[Im((-4 (z^10)+4 (z^9)+4 (z^8)+6 (z^7)-7 (z^6)+7 (z^5)+5 (z^4)+(z^3)+6 (z^2)+4 z-5)^(1/7)/.z->x+I y),{x,-(3/2),3/2},{y,-(3/2),3/2},ColorFunction->"BeachColors",Mesh->None,ImageSize->500,MaxRecursion->2,Evaluated->True]

Algebraic Surfaces

ContourPlot3D can detect all the singular points for algebraic surfaces.

SetOptions[ContourPlot3D,Mesh->None,Boxed->False,BoundaryStyle->None,Axes->None,ImageSize->500,Background->GrayLevel[0.85]];

ContourPlot3D[x^2+y^2-(1-z) z^4==0.0,{x,-1,1},{y,-1,1},{z,-1,1.001},RegionFunction->(#1^2+#2^2<1&),PlotPoints->33,MaxRecursion->4,ContourStyle->{Specularity[White,50],RGBColor[0.15,0.45,0.94]}]

ContourPlot3D[4 x^2 (x^2+y^2+z^2)-y^2 (1-y^2-z^2)==0,{x,-1,1.001},{y,-1,1},{z,-1,1.02},PlotPoints->33,Mesh->1,MeshFunctions->{#2&},MaxRecursion->4,MeshShading->{{Specularity[White,50],RGBColor[0.15,0.45,0.94]},{Specularity[White,50],RGBColor[0.94`,0.90,0.16]}}]

ContourPlot3D[(x^4)-(x^2)+(y^4)+(z^4)-(y^2)-(z^2)+0.5==0,{x,-1,1},{y,-1,1},{z,-1,1},PlotPoints->33,ContourStyle->{Specularity[White,50],Blue},Lighting->{{White,ImageScaled[{1,1,1}]},{White,ImageScaled[{-1,-1,-1}]}}]

Surface Surface Intersections

Manipulate[Plot3D[{x^2+y^2,x+r y+a},{x,-10,10},{y,-10,10},Mesh->None,MaxRecursion->2,PlotRange->All,BoundaryStyle->{1->Red,2->Yellow, {1,2}->{Thickness[0.01],Blue}},PerformanceGoal->"Quality",PlotPoints->15],

{{a,10},10,30,0.25},
{{r,5},-20,20,0.5}]

Manipulate[Plot3D[{x+y+a,Im[Sqrt[x+y I]]},{x,-1,1},{y,-1,1},PerformanceGoal->"Quality",ClippingStyle->Yellow,BoundaryStyle->{Red,Green,{1,2}->{Thickness[0.01],Blue}},PlotRange->{-2,1.2}],

{{a,0},-2,2,0.1}]

Manipulate[ContourPlot3D[{x^2 + y^2 + z^2 == r, x^2 + y^3 == z^2}, {x, -2, 
  2}, {y, -2, 2}, {z, -2, 2}, Mesh -> None, 
 ContourStyle -> {Red, Yellow}, 
 BoundaryStyle -> {2 -> None, {1, 2} -> {Tube[0.12], Blue}}],
{{r,2},0.1,3,0.1}
]

Data Visualization

Several of the data based visualization functions support interpolation.

data1=Table[N[5-Log[x+1]]+Boole[Mod[x,3]==0]RandomReal[{-1,1}],{x,0,35,1}];

data1i=Interpolation[data1,InterpolationOrder->2];

Plot[data1i[x],{x,1,Length[data1]},PlotLabel->Style["Polynomial",FontColor->Black,FontFamily->"Trebuchet MS",FontSize->24],Mesh->{Range[Length[data1]]},MeshFunctions->{#1&},MeshStyle->{{Black,Opacity[0.6],PointSize[0.01]}}]

GraphicsRow[Table[ListLinePlot[data1 ,InterpolationOrder->n,ImageSize->300,Mesh->Full,MeshStyle->{{Black,Opacity[0.6],PointSize[0.01]}}],{n,0,2}],ImageSize->900]

Irregular Data

It is possible to plot irregular data as easily as regular data.

ReliefPlot[Helen,ColorFunction->"BeachColors",BoxRatios->{1,1,0.6}]

Manipulate[ListPlot3D[Helen,MaxPlotPoints->50,ColorFunction->"BeachColors",MeshFunctions->{#3&},Mesh->mmesh],{mmesh,{None,All,{ctrs}->"Multiple Contours"}}]

Data do not have to be on a regular grid, here we randomly sample data from elevation data.

rdata = SampleData[Helen,100];

Manipulate[ListPlot3D[rdata,ColorFunction->"BeachColors",MeshFunctions->{#3&},Mesh->mmesh],{mmesh,{None,All,{ctrs}->"Multiple Contours"}}]

For irregular data, the interpolation supported are : linear, Voronoi cells and Natural Neighbor.

Manipulate[GraphicsGrid[{
{
ListContourPlot[rdata,ColorFunction->bmap,Contours->ctrs,PlotRange->All,Mesh->All],ListContourPlot[rdata,ColorFunction->bmap,PlotRange->All,Contours->ctrs]},
{
ListContourPlot[rdata,ColorFunction->bmap,Contours->ctrs,InterpolationOrder->0,PlotRange->{{50,450},{50,350}},Mesh->All],
ListContourPlot[rdata,ColorFunction->bmap,Contours->ctrs,MaxPlotPoints->50,InterpolationOrder->1,PlotRange->{{50,450},{50,350}}]}},ImageSize->550],

	{{bmap,"BeachColors"},ColorData["Gradients"]}
]

Delaunay and Voronoi Diagrams

The basic operation needed to visualize irregular data is to construct a Delaunay triangulation.

data = Table[{x=RandomReal[{-1,1}],y=RandomReal[{-1,1}],Exp[-x^2-y^2]},{300}];

ListDensityPlot[data,InterpolationOrder->None,Mesh->All,ColorFunction->"BeachColors"]
ListDensityPlot[data,InterpolationOrder->0,Mesh->All,PlotRange->Automatic,ColorFunction->"BeachColors"]
ListPlot3D[data,Mesh->All,InterpolationOrder->0,ColorFunction->"BeachColors",Filling->Bottom,FillingStyle->Pink]
ListPlot3D[data,Mesh->{Full,Full},InterpolationOrder->0,ColorFunction->"BeachColors",Filling->None,FillingStyle->Pink]

The Delaunay and Voronoi Diagrams construction is highly optimized and can deal with highly complex data.

ll=N[stars[7,6]];
ll=Union[ll];
ll=LL[#]& /@ ll;

ListDensityPlot[ll,Mesh->All,ColorFunction->"BeachColors",ImageSize->600]
ListDensityPlot[ll,Mesh->All,InterpolationOrder->0,ColorFunction->"BeachColors",ImageSize->600]

ll=N[stars[17,4]];
ll=Union[ll];
ll=LL[#]& /@ ll;

ListDensityPlot[ll,Mesh->All,ColorFunction->"BeachColors",ImageSize->600]
ListDensityPlot[ll,Mesh->All,InterpolationOrder->0,ColorFunction->"BeachColors",ImageSize->600]

Summary

Same concepts and methods consistent in all visualization functions.

Great flexibility and a plethora of new features

High quality graphics with a minimum amount of coding

State of the art algorithms and spatial data structures

Combination of numerical and symbolic methods provide the new set of visualization functions with unique features not found in any other algebraic computational system.

Initialization

POrange = RGBColor[228./255,111./255,0./255];
POrangeO = RGBColor[228./255,111./255,0./255,0.3];
bmap="Breeze";

c1={143,79,186}/255.;
mypurple=RGBColor[c1];
c2 = {217,34,135}/255.;
myred = RGBColor[c2];
c3 = {188,243,0}/255.;
mygreen = RGBColor[c3];
c4={59,94,255}/255.;
myblue=RGBColor[c4];
c5={251,255,0}/255.;
myyellow=RGBColor[c5];
c6 = {0,168,245}/255.;
mylightblue=RGBColor[c6];

Clear[x,y]
ComparePlot[f_, x_, opts___]:=Grid[{{
Plot[f,x,opts, Method->{"pw"->False},PlotStyle->{POrange, Thickness[0.005]},ImageSize->400,Exclusions->None,FillingStyle->POrangeO],Plot[f,x,opts, PlotStyle->{POrange, Thickness[0.005]},FillingStyle->POrangeO,ImageSize->400,Exclusions->None,Method->{"pw"->True}]}}]

plotstyle = FaceForm[{Opacity[0.8],POrange}, {Opacity[0.8],Purple}];
ComparePlot3D[f_, x_, y_, opts___]:=Grid[{{
Plot3D[f,x,y,opts,PlotStyle->plotstyle,BoundaryStyle->None, ImageSize->280,Exclusions->None,Boxed->False,Axes->None,Method->{"pw"->False},Mesh->All],Plot3D[f,x,y,opts,PlotStyle->plotstyle,BoundaryStyle->None, ImageSize->280,Exclusions->None,Boxed->False,Axes->None,Mesh->All],Plot3D[f,x,y,opts,PlotStyle->plotstyle,BoundaryStyle->None, ImageSize->280,Exclusions->None,Boxed->False,Axes->None,Mesh->5]
}}]

CompareDensityPlot[f_, x_, y_, opts___]:=Grid[{{
DensityPlot[f,x,y,opts,ColorFunction->"BeachColors", BoundaryStyle->None, ImageSize->280,Exclusions->None,Frame->False,Axes->None,Method->{"pw"->False},Mesh->All],DensityPlot[f,x,y,opts,ColorFunction->"BeachColors", BoundaryStyle->None, ImageSize->280,Exclusions->None,Frame->False,Axes->None,Mesh->All,Method->{"pw"->True}],
DensityPlot[f,x,y,opts,ColorFunction->"BeachColors", BoundaryStyle->None, ImageSize->280,Exclusions->None,Frame->False,Axes->None,Mesh->5,Method->{"pw"->True}]
}}]

plotstyleexcl = FaceForm[{Opacity[0.8],POrange}, {Opacity[0.8],Purple}];
imagesizeexc = 200;
cfunctionexc[v_]:=RGBColor[v,v,0];
CompareExclusions[f_, x_, y_, opts___]:=Grid[{{
Plot3D[f,x,y,opts,PlotStyle->plotstyleexcl,ImageSize->imagesizeexc,Exclusions->None,Boxed->False,Axes->None,Mesh->All,Method->{pw->False}],
DensityPlot[f,x,y,opts,ImageSize->imagesizeexc,Exclusions->None,Frame->False,Axes->None,Mesh->None,ColorFunction->"BeachColors"],
Plot3D[f,x,y,opts,PlotStyle->plotstyleexcl,ImageSize->imagesizeexc,Exclusions->Automatic,Boxed->False,Axes->None,Mesh->4],
DensityPlot[f,x,y,opts,ImageSize->imagesizeexc,Exclusions->Automatic,Frame->False,Axes->None,Mesh->None,ColorFunction->"BeachColors"]

}}]

imagesizeexc2 = 270;
CompareExclusions2[f_, x_, y_, opts___]:=Grid[{{
Plot3D[f,{x,- Pi, Pi},{y,-2,2},Exclusions->None,opts,PlotStyle->plotstyleexcl,ViewPoint->{-3,-3,1},Mesh->All,BoxRatios->Automatic,Boxed->False,Axes->None,ImageSize->imagesizeexc2],
DensityPlot[f,{x,- Pi, Pi},{y,-2,2},Exclusions->None,opts,Mesh->None,Frame->False,Axes->None,ImageSize->imagesizeexc2,BoundaryStyle->Red,ColorFunction->"BeachColors"]},{
Plot3D[f,{x,- Pi, Pi},{y,-2,2},opts,PlotStyle->plotstyleexcl,ViewPoint->{-3,-3,1},BoxRatios->Automatic,Axes->None,Boxed->False,BoundaryStyle->Red,ImageSize->imagesizeexc2,Mesh->{4,10}],DensityPlot[f,{x,- Pi, Pi},{y,-2,2},opts,Frame->False,Axes->None,Exclusions->{{x==-Pi/2,Abs[y]>cut},{x==Pi/2,Abs[y]>cut}},BoundaryStyle->Red,ImageSize->imagesizeexc2,ColorFunction->"BeachColors"]
}}]

binimagesize={400,Automatic};
binplotstyle=FaceForm[Directive[{Opacity[0.85],POrange}],Directive[{Opacity[0.85],Purple}]];
Compare2Plot3D[f_, x_, y_, opts___]:=Grid[{{
Plot3D[f,{x,-3.5,1.5},{y,-2.5,1.5},opts,PlotStyle->binplotstyle,Boxed->False,Axes->None,PlotRange->{-6,6},ImageSize->binimagesize,Mesh->All,ViewPoint->{0.7,-2.3,1.5}],
Plot3D[Binomial[x,y],{x,-3.5,1.5},{y,-2.5,1.5},opts,PlotStyle->binplotstyle,Mesh->{Range[-6.001,6,1]},MeshFunctions->(#3&),MeshStyle->Black,Boxed->False,Axes->None,ClippingStyle->False,PlotRange->{-6,6},ViewPoint->{0.7,-2.3,1.5},ImageSize->binimagesize,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Exclusions->{x==-1,x==-2,x==-3},MaxRecursion->3],Graphics3D[{},Boxed->False,Axes->None]}
}]


bticks[xmin_,xmax_,step_,plen_,nlen_,size_,col_,op_]:=Table[{i,Style[i,FontColor->col,FontFamily->"Helvetica",FontSize->size],{plen,nlen},{col,Opacity[op]}},{i,xmin,xmax,step}];
xticks[xmin_,xmax_,step_]:=bticks[xmin,xmax,step,.02,.02,16,Black,.4];
yticks[ymin_,ymax_,step_]:=bticks[ymin,ymax,step,.98,.02,16,Black,.2];
fticks[xmin_,xmax_,step_]:=bticks[xmin,xmax,step,.01,0.,15,Black,1.];
cpts[data_,col_,op_,size_]:=Sequence[Mesh->{Range[Length[data]]},MeshFunctions->{#1&},MeshStyle->{{col,Opacity[op],PointSize[size]}}];
cpts3d[data_,col_,op_,size_]:=Graphics3D[{col,Opacity[op],PointSize[size],Flatten[Table[Point[{j,i,data[[i,j]]}],{i,Length[data]},{j,Length[First[data]]}]]}];
mesh3d[data_]:={Length[First[data]]-2,Length[data]-2};
terrainColor[v_]:=Blend[{Brown,Darker[Green,.75],White},v+RandomReal[.2]];

Clear[plot1style,image1style,label1style,plot2style];
plot1style[clist_]:=Table[{Thickness[.005],Opacity[.7],Darker[clist[[i]],.3]},{i,Length[clist]}];
image1style[min_,max_]:=Sequence[Filling->True,ImageSize->500,ImageMargins->15,Background->Lighter[Gray,.9],PlotRange->{min,max},AxesOrigin->{0,0}];
label1style=Sequence[FontColor->Black,FontFamily->"Trebuchet MS",FontSize->24];
plot2style=Sequence[PlotStyle->{Yellow,Specularity[White,50],Opacity[.7]},ImageSize->450,ImageMargins->15,Boxed->False,Axes->True];
plot3style=Sequence[ImageSize->400,ImageMargins->15,Frame->True,Axes->False];

Data

$DataDirectory=ToFileName[{NotebookDirectory[]}];
AppendTo[$Path,$DataDirectory];
rawHelen=Developer`ToPackedArray[Get["Helens.dat"]];
cx=600; cy=300; mm=500; nn=400;
Helen=Developer`ToPackedArray[Table[rawHelen[[cx+i,cy+j]],{i,0,mm},{j,0,nn}]];

SampleData[data_,n_]:=Module[{dx,dy, temp, temp2,i,j},
{dx,dy} = Dimensions[data];
temp = Table[{i=RandomInteger[{1,dx}],j=RandomInteger[{1,dy}],data[[i,j]]},{n}];
temp2 = {{1,1,data[[1,1]]},{i,dy,data[[1,dy]]},{dx,dy,data[[dx,dy]]},{dx,1,data[[dx,1]]}};
Union[Join[temp,temp2]]
]

ctrs = Table[i,{i,5000,7800,200}];
blendmap = "Topographic";

stars[side_,0]:={{0,0,1}};
stars[sides_,lvl_]:=Union[Flatten[Chop[(Map[starhelp[sides,#]&,stars[sides,lvl-1]])],1]]
starhelp[sides_,center_List]:=Table[center+{Re[#],Im[#],1}&[Exp[2 Pi I k/sides]],{k,0,sides-1}]
LL[{x_,y_,z_}]:={x,y,Norm[{x,y}]}

Default Styles.

SemStyle1 := Sequence[{Boxed->False,MeshStyle->Black,Mesh->All,Axes->None,BoundaryStyle->None}];

Lastly

$Line=0;