Message Boards

GROUPS:

Staff Picks Data Science Mathematics Recreation Visual Arts Geometry Graphics and Visualization Wolfram Language Optimization Machine Learning

Silvia Hao

[GIF] Ribbons of three materials: splines, neural nets & contour surfaces

Silvia Hao, CTO at Glimscape Technology

Posted 25 days ago

1960 Views

7 Replies

23 Total Likes

Follow this post

( ^ Click the GIF for a 3840x1280 wallpaper! :)

— or how to fit B-spline onto contour surface

Happy New Year to everyone who reached this post! Today we are going to explore a mixed power of B-spline and neural network!

When plotting contour surface, have you, like me, done things like increasing PlotPoints so to have smooth looking surface and mesh lines, in the cost of bloomed up notebook size and/or less responsive Graphics3D result? If the answer is yes, then please do take a look at the approach of B-spline. As a nice function approximator, we may not need that much control points to have a applause B-spline if it’s just for meeting the eyes. And for RAM efficiency, we are going to do the computing with differentiable programming, through Wolfam Language’s neural network framework.

Helper functions

ClearAll[FindDivisionsExact]FindDivisionsExact[range:{_,_},n_Integer?(#>=2&)]:=Rescale[Range[n]//N//Rescale,{0,1},range]FindDivisionsExact[n_Integer?(#>=2&)]:=FindDivisionsExact[{0,1},n]

ClearAll[pipe,branch,branchSeq]pipe=RightComposition;branch=Through@*{##}&;branchSeq=pipe[branch@##,Apply@Sequence]&;

The contour surface model — a well-known example

For demonstration I’ll use one of my favourite example, the transcendental periodic implicit surface in ContourPlot3D’s documentation, constrained in a ball region:

ClearAll[,ℛ]=Function[{x,y,z},Cos[x]Sin[y]+Cos[y]Sin[z]+Cos[z]Sin[x]];ℛ=Ball[{0,0,0},π];

contourGraph=ContourPlot3D[[x,y,z]0,{x,y,z}∈ℛ,PlotPoints40,MaxRecursion0,MeshStyleGrayLevel[.7],RegionBoundaryStyleNone,BoxedFalse,AxesFalse,Lighting"Accent",ContourStyleWhite]

We are going to find mesh lines (curves) on  with their terminals being randomly selected from the surface of ball ℛ:

ClearAll[bℛ]bℛ=And[[x,y,z]0,{x,y,z}∈RegionBoundary[ℛ]]//pipe[ImplicitRegion[#,{x,y,z}]&,DiscretizeRegion[#,BoxedTrue,AxesTrue]&]

The B-spline curve model

As we have talked about in past posts, a uniform B-spline curve is simply a linear combination of its control points with BSplineBasis as the weights:

C(t):=

∑

i=1

i,d

(t),t∈[0,1]

(

)

For our purpose here, we discretize the curve

C(t)

by a naïve uniform sample over

t∈[0,1]

, so it becomes:

c×s×d

:=

s×n

n×c×d

,

(

)

where

is the discrete B-spline basis matrix,

is the array of control points. The meaning of their dimensions are

Symbol	Meaning	Variable name	Typical value
c	number of curves	$curveNum	11
s	sample size	$sampleNum	100
d	embedded dimension	$embeddedDim	3
n	number of control points	$CPtsNum	5

Prototype for the curve

As we have shown, the curves are represented as

C=N·P

. But we can’t just throw a LinearLayer for it because we want the terminals predefined and fixed while keeping the “internal” control points trainable. So instead we separate them as

c×s×d

:=

s×n

·Join

startingterminal

1×c×d

middle

(n-2)×c×d

endingterminal

1×c×d



(

)

<<NeuralNetworks`

ClearAll[curvePrototype]curvePrototype=NetGraph[(*s*n*)"basis"PlaceholderLayer[InputPorts{},OutputPorts{MatrixT[]}],"startCpts"ReplicateLayer[1],"endCpts"ReplicateLayer[1],(*(n-2)*c*d*)"midCpts"NetArrayLayer["Output"TensorOfRankT[3]],(*n*c*d*)"cpts"CatenateLayer[],"curve"DotLayer[],(*c*s*d*)"curveSet"TransposeLayer[],{Thread[(NetPort/@{"start","end"}){"startCpts","endCpts"}],{"startCpts","midCpts","endCpts"}"cpts",{"basis","cpts"}"curve""curveSet"}//Flatten,(*c*d*)"start"MatrixT[],(*c*d*)"end"MatrixT[]]

Replacing “basis” and “midCpts” layers with concrete NetArrayLayer will instantize the prototype:

Module[{$sampleNum=5,$CPtsNum=3,$curveNum=11,$embeddedDim=2,data,net},net=NetReplacePart[curvePrototype,{"basis"NetArrayLayer["Array"RandomReal[1,{$sampleNum,$CPtsNum}]],{"midCpts","Array"}RandomReal[1,{$CPtsNum-2,$curveNum,$embeddedDim}]}];data=RandomReal[1,{2,$curveNum,$embeddedDim}]//AssociationThread[{"start","end"},#]&;data//pipe[net,MapThread[{#2,Line@#1}&,{#,#//Length//FindDivisionsExact//Map[ColorData@"Rainbow"]}]&,Graphics]]

B-spline basis

To specify a set of B-spline basis through BSplineBasis, we need to firstly decide how “smooth” we want the curves to be (the degree of the spline $bsDegree), and how capable of fitting wiggle shapes we want the curves to be (the number of control points for each spline $CPtsNum). Then the uniform knots can be determined unambiguously. To discretize the basis we need to specify the sample size $sampleNum. The basis set will then be generated and stored in an un-trainable NetArrayLayer:

ClearAll[basisNet]basisNet[$sampleNum_Integer?Positive,{$CPtsNum_Integer?Positive,$bsDegree_Integer?Positive,$fullKnot_List}]:=Table[BSplineBasis[{$bsDegree,$fullKnot},i-1,#],{i,$CPtsNum}]//pipe[Map@Function,Apply@branch,Map,branchSeq[Identity,FindDivisionsExact[$sampleNum]&],Construct,NetArrayLayer["Array"#,LearningRateMultipliersNone]&]

Clear[$bsDegree,$CPtsNum,$bsKnot,$sampleNum]$bsDegree=3;$CPtsNum=4;If[Not[$CPtsNum-$bsDegree+1≥2],Echo@"$CPtsNum and $bsDegree is not compatible."];$bsKnot={ConstantArray[0,$bsDegree],Range[$CPtsNum-$bsDegree+1]//Rescale,ConstantArray[1,$bsDegree]}//Flatten;$sampleNum=100;

basisNet[$sampleNum,{$CPtsNum,$bsDegree,$bsKnot}]//pipe[Construct,Transpose,ListLinePlot]

Instance of curves

Replacing “basis” layer with basisNet will instantize a group of B-spline curves:

Module{$curveNum=5,$embeddedDim=2,net,data},net=NetReplacePart[curvePrototype,{"basis"basisNet[$sampleNum,{$CPtsNum,$bsDegree,$bsKnot}],{"midCpts","Array"}RandomReal[1,{$CPtsNum-2,$curveNum,$embeddedDim}],"start"{$curveNum,$embeddedDim}}];Echo[net];data=RandomReal[1,{2,$curveNum,$embeddedDim}]//AssociationThread[{"start","end"},#]&;net//pipebranch[Identity,pipe[NetTake[#,"cpts"],Transpose]&],Apply@branch,#[data]&,MapThread(*GroundtruthB-splinecurve:*){#2,BSplineCurve[#1〚2〛,SplineDegree$bsDegree,SplineKnots$bsKnot]},(*NNcomputedcurve:*){#2,Dashed,AbsoluteThickness[3],Line@#1〚1〛},(*Controllineandpoints:*){Lighter@#2,Line@#1〚2〛},EdgeForm[Lighter@#2],FaceForm

,Translate[Disk[{0,0},.01],#]&@#12&,{#,#//Length//FindDivisionsExact//Map[ColorData@"Rainbow"]}&,Graphics

The contour surface as a neural network

On the surface side, we are going to construct a neural network that takes a group of

(d=$embeddedDim=3)

dimensional points as input, gives a total score on how well those points fit on the surface:

ClearAll[$embeddedDim,surfaceMatchNet]$embeddedDim=3;surfaceMatchNet=Module[{samplePoints},samplePoints=StringTemplate["x``"]/@Range[$embeddedDim];NetGraph[(*Destructinputcurvearraytodports:*)"curve2pts"NetGraph[Flatten@{PartLayer/@Range[$embeddedDim],TransposeLayer[1-1]},Flatten@{($embeddedDim+1)Range[$embeddedDim],Thread[Range[$embeddedDim](NetPort/@samplePoints)]}],(*Feedintoformulaofthesurface:*)"surface"ThreadingLayer@(*Fitgoodnessasameansquaredloss:*),"fit goodness"{ElementwiseLayer[

&],SummationLayer[]},{NetPort["Input"]"curve2pts",(NetPort["curve2pts",#]&/@samplePoints)"surface""fit goodness"},"Input"TensorWithMinRankT[2]]]

Our dimension setting TensorWithMinRankT[2] is general enough, so surfaceMatchNet works as expected for both

s×d

and

c×s×d

sample point set:

isoSurf=NetTake[surfaceMatchNet,"surface"];

Example of

points:

RandomPoint[ℛ,

]//pipe[branchSeq[Identity,pipe[EchoTiming[isoSurf@#,"isoSurf timing"]&,Thread[

<0.01]&]],Pick,Graphics3D[{AbsolutePointSize[1],Point@#}]&,AbsoluteTiming]

⌚

isoSurf timing0.0397876

0.249484,



Example of 8 groups with

points per group:

Tuples[{-1,1},3]//pipe[Map@pipe[{{0,0,0},#π}&,N,Map@Sort,Transpose,Apply@Cuboid,Function[octant,RegionIntersection[ℛ,octant]//RandomPoint[#,

]&]],Developer`ToPackedArray,branchSeq[Identity,pipe[EchoTiming[isoSurf@#,"isoSurf timing"]&,UnitStep[.01-

]&,Transpose]],Pick[##,1]&,MapThread[{#2,Point@#1}&,{#,#//Length//FindDivisionsExact//Map[ColorData@"Rainbow"]}]&,Graphics3D[{AbsolutePointSize[1],#}]&,AbsoluteTiming]

⌚

isoSurf timing0.0251432

1.64324,



Solve the fitting problem with neural network optimizer

Now we have a curve model curvePrototype that outputs

curves in

dimensional space with sample size

, in a shape of

c×s×d

array, it’s easy to pipe its output to the contour surface measurement net surfaceMatchNet, so to get a whole training target on how well the B-spline curves lie on the surface:

NetGraph[{curvePrototype,surfaceMatchNet},{12}]//Information[#,"FullSummaryGraphic"]&

We select

2×c

random points from the boundary region

bℛ

as the terminal points of the

curves, and ask for a fit-goodness equals zero:

Clear[$curveNum,terminalPointSet,trainingData]$curveNum=100;terminalPointSet=RandomPoint[bℛ,{2,$curveNum}];trainingData=terminalPointSet//pipe[Append[0],AssociationThread[{"start","end","fit goodness"},#]&,Map[Developer`ToPackedArray]];

Generally, higher spline degree and more control points result in better fitting:

Clear[$bsDegree,$CPtsNum,$bsKnot,$sampleNum]$bsDegree=5;$CPtsNum=7;If$CPtsNum-$bsDegree+1≥2,Echo@Style"$CPtsNum and $bsDegree are compatible.",Bold,

,Echo@Style"$CPtsNum and $bsDegree are not compatible.",Bold,

;$bsKnot={ConstantArray[0,$bsDegree],Range[$CPtsNum-$bsDegree+1]//Rescale,ConstantArray[1,$bsDegree]}//Flatten;$sampleNum=100;

$CPtsNum and $bsDegree are compatible.

nnResult=Module[{midCpts,net,trained},midCpts=FindDivisionsExact[$CPtsNum]〚2;;-2〛//{1-#,#}.{trainingData@"start",trainingData@"end"}&//N;net=Module[{curvenet},curvenet=NetReplacePart[curvePrototype,{"basis"basisNet[$sampleNum,{$CPtsNum,$bsDegree,$bsKnot}],{"midCpts","Array"}midCpts,"start"{$curveNum,$embeddedDim}}];NetGraph["curve"curvenet,"surface"surfaceMatchNet,{"curve""surface"NetPort["fit goodness"]}]];Off[NetTrain::novalidation];trained=NetTrain[net,{trainingData},All,LossFunction"fit goodness",BatchSize1,TimeGoalQuantity[10,"Minutes"],TrainingStoppingCriterion"Criterion""Loss",WorkingPrecision"Real64",Method"RMSProp"];On[NetTrain::novalidation];trained];//AbsoluteTiming

{0.782218,Null}

Mean fit-goodness per point:

nnResult["RoundLoss"]

$curveNum$sampleNum

0.0024313

We can see the result indeed fits well to the contour surface:

nnResult["TrainedNet"]//pipe[NetPart[#,1]&,NetTake[#,"cpts"]&,#@trainingData&,Transpose,MapThread[{#2,BSplineCurve[#1,SplineDegree$bsDegree,SplineKnots$bsKnot]}&,{#,#//Length//FindDivisionsExact//Map@ColorData@"Rainbow"}]&,Graphics3D[{AbsoluteThickness[2],#}]&,Show[contourGraph,#]&]

It takes two splines to make a ribbon

Last section gives us a convenient way to fit BSpline curve onto contour surface. And one pair of such curves with sufficient close shape will make a nice ribbon as tensor product BSplineSurface. (Please refer to my last community post.)

Say we are creating

group

200

ribbon

group

, then we can firstly initialize it with

group

200

curve

group



1000

curves

, which are

group

200

curve

group

terminal point

curve



2000

terminal points

$curveNum=5×200×2;

terminalPointSet1=RandomPoint[bℛ,{5,200,2}];

Graphics3D[{AbsolutePointSize[1],Point@Flatten[terminalPointSet1,2]},ImageSizeSmall]

For each terminal point

in terminalPointSet1, we look for its accompanying point

′

such that

′

∈bℛ∧

′

-P=

width

ribbon

=0.1

ribbonWidth=.1;terminalPointSet2=Map[Function[tp,{x,y,z}/.FindRoot[{[x,y,z]0,Norm[{x,y,z}]π,Norm[{x,y,z}-tp]ribbonWidth},{{x,y,z},tp+.01}]],terminalPointSet1,{3}];

{terminalPointSet1,terminalPointSet2}//pipe[Transpose[#,{4,1,2,3,5}]&,Flatten[#,2]&,MapThread[{#2,Line[#1]}&,{#,#//Length//FindDivisionsExact//Map@ColorData["Rainbow"]}]&,{AbsoluteThickness[20],CapForm[None],#}&,Flatten,Graphics3D[#,ImageSizeMedium]&]

Then we combine both set into single one suitable for previous neural net:

terminalPointSet={terminalPointSet1,terminalPointSet2}//Transpose[#,{3,1,2,4,5}]&//Flatten[#,2]&//Transpose;

terminalPointSet//Dimensions

{2,2000,3}

Here the dimension is changed from



curve

ribbon

group

200

ribbon

group

terminal point

curve

dimension

point





group

200

ribbon

group

curve

ribbon

terminal point

curve

dimension

point





group

200

ribbon

group

curve

ribbon

2000

curve

terminal point

curve

dimension

point





terminal point

curve

2000

curve

dimension

point



, which fits the required data shape

2×c×d

for the training data of the final neural net in last section.

A similar training process with sufficiently small learning rate (

∼

-5

) should produce something like

$bsDegree=3;$CPtsNum=4;$bsKnot={ConstantArray[0,$bsDegree],Range[$CPtsNum-$bsDegree+1]//Rescale,ConstantArray[1,$bsDegree]}//Flatten;

resultForRibbons=

result for ribbons

Head: List

Dimensions: {2000,4,3}

Byte count: 192216

;

which has a shape as



2000

curve

control point

curve

dimension

point



Visualizations

As we mentioned before, ribbons are visualized as BSplineSurface:

resultForRibbons//pipe(*partitioncontrol-pointsetintogroupsandribbons:*)ArrayReshape[#,{5,200,2,4,3}]&,pipe[branchSeq[Identity,(*computedistanceofthe2ctrl-pointssetofthe2edge-curvesofthesameribbon:*)Map[pipe[{1,-1}.#&,Map@Norm,Total,#<5&],#,{2}]&],(*selectribbonswithsufficientlynarrowwidth:*)Pick],Map@pipe[(*visualizeasinglegroup:*)Map[(*visualizeasingleribbon:*)BSplineSurface[#,SplineDegree{Automatic,$bsDegree},SplineKnots{Automatic,$bsKnot}]&]],(*placewherethestylinghappens:*)EdgeForm[],ToonShading

,#&,Graphics3D[#,Lighting"Neutral",BoxedFalse,ViewVector{{-28.5,-12.5,5},{0,0,-0.2}},ViewVertical{-0.42,0.21,0.88}]&

Playing with the different shadings and Lighting etc. gives all kinds of interesting visualizations:

resultForRibbons//pipe(*partitioncontrol-pointsetintogroupsandribbons:*)ArrayReshape[#,{5,200,2,4,3}]&,pipe[branchSeq[Identity,(*computedistanceofthe2ctrl-pointssetofthe2edge-curvesofthesameribbon:*)Map[pipe[{1,-1}.#&,Map@Norm,Total,#<5&],#,{2}]&],(*selectribbonswithsufficientlynarrowwidth:*)Pick],Map@pipe[(*visualizeasinglegroup:*)MapThread[(*visualizeasingleribbon:*)(*placewherethestylinghappens:*){EdgeForm[],FaceForm[{Specularity[1,10],Opacity[.3],#2},None],BSplineSurface[#1,SplineDegree{Automatic,$bsDegree},SplineKnots{Automatic,$bsKnot}]}&,{#,#//Length//FindDivisionsExact//Map[ColorData["DeepSeaColors"]]}]&],Graphics3D#,Lighting"Ambient",

,Background

,BoxedFalse,ViewVector{{-21.7,16.2,-3.37},{0,0,-0.2}},ViewVertical{0.05,0.89,0.46}&

Or finner ribbons for StippleShading or HalftoneShading:

resultForStipple=

result for stippling and halft