Re: I'd appreciate some help with some forulae for mathematical knots
- To: mathgroup at smc.vnet.net
- Subject: [mg5144] Re: [mg5115] I'd appreciate some help with some forulae for mathematical knots
- From: jpk at apex.mpe.FTA-Berlin.de (Jens-Peer Kuska)
- Date: Wed, 6 Nov 1996 01:33:10 -0500
- Sender: owner-wri-mathgroup at wolfram.com
Dear Alan, - for the smooth surfaces I had written a package - Mesh3D.m (attatched) and I had modifyed the POVray.m (I attatch it also) package by R. M"ader. Mesh3D introduces a new graphics object called Mesh3D. If You convert a graphics 3d object to a Mesh3D object the points are stored seperatly and for every point a averaged normal is calculated. The polygons gets triangulated. You can still convert the Mesh3D object back to a Graphics3D and check the triangulation. The modifyed POVray["filename",Mesh3D[graph3d]] writes out smoothtriangles with the normals at every point of the surface. This produces the smooth surfaces You want. I had used it for surfaces with 1000-2000 polygons. For sutch a large number of polygons it is a bit slow. But for Your knots I see no problem. You need POVray 3.0 from http://www.povray.org to get all working. Contact me if You need more information or a sample notebook. Hope that helps Jens ---------- BeginPackage["Mesh3D`",{"Geometry`Rotations`","Utilities`FilterOptions`"}] Mesh3D::usage= "Mesh3D[Graphics3D[...]] will produce a graphics 3d like object. The first argument of Mesh3D are the list of points, the second the surface normals at the points. The polygons, lines and points stored as index lists. Polygons are usualy triangulated when a graphics 3d object is convertet." TranslateMesh::usage= "TranslateMesh[mesh3d,vector] will translate the mesh3d object." RotateMesh::usage= "RotateMesh[mesh3d,phi,theta,psi] will rotate the mesh3d object by the Euler angles." Begin["`Private`"] (* triangulate polygons *) Clear[triangulate] triCount=0; triangulate[Polygon[pnts_ /; Length[pnts]<=3]]:=Polygon[pnts] triangulate[Polygon[{p1_,p2_,p3_,p4_}]]:= If[EvenQ[triCount++], {Polygon[{p1,p2,p3}],Polygon[{p3,p4,p1}]}, {Polygon[{p2,p3,p4}],Polygon[{p4,p1,p2}]} ] triangulate[Polygon[pnts_]]:= Module[{c}, c=(Plus @@ pnts)/Length[pnts]; Polygon[{c,#[[1]],#[[2]]}] & /@ Transpose[{pnts,RotateLeft[pnts]}] ] (**********************************************) (* Extract point data from graphic primitives *) (**********************************************) GetPointPool[Polygon[pnts_]]:=Point /@ pnts GetPointPool[Line[pnts_]]:=Point /@ pnts GetPointPool[p_Point]:=p GetPointPool[any_List ]:=GetPointPool /@ any GetPointPool[any_]:={} (**********************************************) (* Calculate the averaged normal vectors *) (**********************************************) (* Make the normal vector for the i-th point in the polygon list polys *) getNeigbours[i_Integer,lst:{_Integer..}]:= If[First[lst]==i, {Last[lst],lst[[2]]}, If[Last[lst]==i, {lst[[-2]],First[lst]}, (* else *) Part[lst,#] & /@ Flatten[ {#-1,#+1} & /@ Flatten[Position[lst,i]] ] ] ] cross[{ax_,ay_,az_},{bx_,by_,bz_}]:= {ay*bz - az*by,az*bx -ax*bz,ax*by-ay*bx} Clear[MakeNormalFor] (* Fluid from Mesh3D : pointPool *) MakeNormalFor[i_Integer, polys_]:= Module[{ipoly,thisp,vecs,norm}, ipoly=Select[polys,MemberQ[#,i] &]; If[ipoly==={}, Return[Null] ]; thisp=pointPool[[i]]; (* For the polygon k *----- i \ \ * l we calculate the pair {k,l} with getNeigbours and then the vectors v1=p[k]-p[i] and v2=p[l]-p[i]*) vecs= Map[ (Part[pointPool,#]-thisp) &, getNeigbours[i,#] & /@ ipoly, {2} ]; (* Normalise the vectors *) vecs= Map[#/Sqrt[Dot[#,#]] &,vecs,{2}]; vecs=Map[Apply[cross,#] &,vecs]; vecs=(Plus @@ vecs)/Length[vecs]; norm=Sqrt[Dot[vecs,vecs]]; vecs/norm ] Clear[Mesh3D] Format[_Mesh3D]:="-Mesh3D-" Mesh3D[surf_SurfaceGraphics,mopt___]:=Mesh3D[Graphics3D[surf],mopt] Mesh3D[Graphics3D[gdata_,opt___],mopt___]:= Block[{data,pointPool,idata,rules,i,pnts,lns,polys,normals}, triCount=0; data=gdata //. { Polygon[{a___,b_,b_,c___}] :> Polygon[{a,b,c}], Polygon[{a_,b___,a_}] :> Polygon[{a,b}] } /. p_Polygon :> triangulate[p]; pointPool=N[Flatten[GetPointPool[data]]]; pointPool=Sort[Union[pointPool]]; rules=MapIndexed[(First[#1]->First[#2]) &,pointPool]; idata=Flatten[N[data]] //. rules; pnts=First /@ Select[idata,Head[#]===Point &]; lns=First /@ Select[idata,Head[#]===Line &]; polys=First /@ Select[idata,Head[#]===Polygon &]; pointPool= First /@ pointPool; normals=MakeNormalFor[#,polys] & /@ Table[i,{i,1,Length[pointPool]}]; Mesh3D[pointPool,normals,pnts,lns,polys,Join[opt,{mopt}]] ] TranslateMesh[b_Mesh3D,{x_?NumberQ,y_?NumberQ,z_?NumberQ}]:= Module[{pointPool}, pointPool=(#+N[{x,y,z}]) & /@ First[b]; Mesh3D[pointPool,Sequence @@ Rest[b]] ] RotateMesh[b_Mesh3D,phi_ /; NumberQ[N[phi]] ,theta_ /; NumberQ[N[theta]],psi_ /; NumberQ[N[psi]]]:= Block[{mat,pointPool=b[[1]],i,normals}, mat=RotationMatrix3D[N[phi],N[theta],N[psi]]; pointPool=Map[mat.# &, pointPool]; normals=MakeNormalFor[#,Last[b]] & /@ Table[i,{i,1,Length[pointPool]}]; Mesh3D[pointPool,normals,Sequence @@ Rest[Rest[b]]] ] Mesh3D /: Graphics3D[Mesh3D[pool_,normals_,pnts_,lns_,polys_,opt_],gopt___]:= Module[{rules,gr}, gr={}; If[pnts=!={}, AppendTo[gr,Point[Part[pool,#]] & /@ pnts ]]; If[lns=!={}, AppendTo[ gr, Line /@ Map[ Part[pool,#] &, lns,{2} ] ] ]; If[polys=!={}, AppendTo[ gr, Polygon /@ Map[ Part[pool,#] &, polys,{2} ] ] ]; Graphics3D[gr,FilterOptions[Graphics3D,Sequence @@ Join[{gopt},opt]]] ] Mesh3D /: Show[mm_Mesh3D,opt___]:= Module[{gr}, Show[Graphics3D[mm],opt]; mm ] End[] EndPackage[] Null ---------- X-Sun-Data-Type: text X-Sun-Data-Name: POVray.m X-Sun-Charset: us-ascii X-Sun-Content-Lines: 374 (* :Title: POV-Ray format conversion *) (* :Author: Roman E. Maeder *) (* :Summary: Converts Mathematica graphics objects into POV-ray input format. *) (* :Context: POVray` *) (* :Package Version: 1.0 *) (* :Copyright: Copyright 1994, Roman E. Maeder. Permission is granted to use and distribute this file for any purpose except for inclusion in commercial software or program collections. This copyright notice must remain intact. *) (* :History: Version 1.0 for the Mathematica Journal, June 1994. modifyed for Mesh3D data structure by J.-P. Kuska August 1996 *) (* :Keywords: POV-Ray, ray tracing *) (* :Source: Maeder, Roman E. 1994. Ray Tracing and Graphics Extensions. The Mathematica Journal, 4(3). *) (* :Warning: only a small number of graphics directives and primitives is converted. *) (* :Mathematica Version: 2.2 *) (* :Discussion: may optionally be used together with the package SurfaceGraphics3D' to render smooth parametric surfaces. *) BeginPackage["POVray`", "Utilities`FilterOptions`"] POVray::usage = "POVray[filename, graphics] writes graphics in POVray format to file. POVray[stream, graphics] uses stream for output." POVHeader::usage = "POVHeader->\"string\" is an option of POVray that gives the string to write at the beginning of the output file." POVInformation::usage = "POVInformation->Automatic is an option of POVray that specifies which extra information is to be included. By default, the light sources and camera position are output. POVInformation->None prints no extra information." Options[POVray] = { POVHeader -> "#default { texture { pigment {color rgb <1,1,1>} finish {phong 0.8 ambient 0.4} }}", POVInformation->Automatic } Begin["`Private`"] Clear[writeobj,POVray]; POVray[filename_String, obj_, opts___] := Module[{file}, file = OpenWrite[filename, FormatType->OutputForm, PageWidth->Infinity]; If[ file === $Failed, Return[file] ]; POVray[file, obj, opts]; Close[file] ] POVray[file_OutputStream, obj_, opts___] := Module[{prolog, extra}, prolog = POVHeader /. {opts} /. Options[POVray]; If[ prolog =!= "", Write[file, prolog] ]; extra = POVInformation /. {opts} /. Options[POVray]; If[ extra === Automatic, extra = True ]; writeobj[file, obj, POVInformation->extra, opts]; file ] newline[file_] := Write[file, ""] (* fix Infix problem for functions with one argument *) protected = Unprotect[Infix] Infix[_[e_], h_:Null] := e Protect[ Evaluate[protected] ] (* optional text followed by numbers, FortranForm with 5 digits *) ndig = 5 nForm[r_] := NumberForm[FortranForm[r], ndig] writeNums[file_, txt_:"", r_?NumberQ] := Write[file, txt, " ", nForm[r]] writeNums[file_, txt_:"", r_List] := Write[file, txt, " ", Infix[nForm /@ r, " "]] (* POV-ray vector notation *) writeVector[file_, txt_:"", r_List,txt1_:""] /; Length[r] == 3 := Write[file, txt, "<", Infix[nForm /@ r, ", "], ">",txt1] (* write various graphics objects *) SetAttributes[writeobj, Listable] (* Graphics3D *) (* Old:: triangles passed as single objects, New:: triangles enclosed by a union { ..} *) writeobj[file_, g:Graphics3D[oblist_, opts_List:{}], myopts___] := Module[{pr, vc, vv, extra}, newline[file]; extra = POVInformation /. {myopts}; If[ extra, pr = FullOptions[g, PlotRange]; vc = FullOptions[g, ViewCenter]; vv = FullOptions[g, ViewVertical]; processOpts[file, pr, ViewCenter->vc, ViewVertical->vv, opts]; ]; newline[file]; Write[file, "/* begin Graphics3D */"]; Write[file," union {\n"]; writeobj[file, Chop[N[oblist]], myopts]; Write[file,"}\n"]; Write[file, "/* end Graphics3D */"] ] (* Original:: polygon points passed verbatim New:: multiple points will removed *) writeobj[file_, Polygon[vl_] , myopts___] := writeTriangle[file, #]& /@ toTriangles[vl /. {a___List,p_List,p_List,c___List} :> {a,p,c}] writeobj[file_, Cuboid[p1_], myopts___] := writeobj[file, Cuboid[p1, p1 + {1,1,1}]] writeobj[file_, Cuboid[p1_, p2_], myopts___] := ( newline[file]; Write[file, "box {"]; writeVector[file, "\t", p1]; writeVector[file, "\t", p2]; Write[file, "}"]; ) writeobj[file_,Point[p_],myopts___] := ( newline[file]; Write[file,"sphere {"]; writeVector[file,"\t",p,", MathPointSize"]; Write[file,"\ttexture{ MathPointTexture}"]; Write[file,"}"]; ) (* SurfaceGraphics3D from SurfaceGraphics3D.m *) writeobj[file_, g0:SurfaceGraphics3D`SurfaceGraphics3D[vl_, norms_, opts_], myopts___] := Module[{pr, vc, vv, extra, g}, newline[file]; extra = POVInformation /. {myopts}; If[ extra, g = Graphics3D[g0]; pr = FullOptions[g, PlotRange]; vc = FullOptions[g, ViewCenter]; vv = FullOptions[g, ViewVertical]; processOpts[file, pr, ViewCenter->vc, ViewVertical->vv, opts]; ]; newline[file]; Write[file, "/* begin SurfaceGraphics3D */"]; Scan[ writeSmoothTriangles[file, #[[1]], #[[2]]]&, Chop[Transpose[{makeMesh[vl], makeMesh[norms]}]] ]; Write[file, "/* end SurfaceGraphics3D */"] ] (* Mesh3D from Mesh3D.m *) writeobj[file_, g0:Mesh3D`Mesh3D[pointPool_,normals_,pnts_,_,poly_,opts_], myopts___] := Module[{pr, vc, vv, extra, g,plst,}, newline[file]; extra = POVInformation /. {myopts}; If[ extra, g = Graphics3D[g0]; pr = FullOptions[g, PlotRange]; vc = FullOptions[g, ViewCenter]; vv = FullOptions[g, ViewVertical]; processOpts[file, pr, ViewCenter->vc, ViewVertical->vv, opts]; ]; newline[file]; Write[file, "/* begin Mesh3D */"]; Write[file, " mesh {"]; plst=poly; While[plst=!={}, p=First[plst]; plst=Rest[plst]; writeSmoothTriangle[ file, Sequence @@ (Flatten[Map[ {pointPool[[#]],normals[[#]]} &, p],1]) ]; ]; Write[file, " }"]; Write[file, "/* end Mesh3D */"] ] (* without normals: turn into Graphics3D without loss of information *) writeobj[file_, g0:_SurfaceGraphics3D`SurfaceGraphics3D, myopts___] := writeobj[file, Graphics3D[g0], myopts] (* SurfaceGraphics: turn into Graphics3D without loss of information *) writeobj[file_, g0:_SurfaceGraphics, myopts___] := writeobj[file, Graphics3D[g0], myopts] (* add other graphics types/primitives here *) (* catch all *) writeobj[file_, _, myopts___] := Null (* unkown object *) makeMesh[vl_List] := Module[{l = Drop[#, -1]& /@ vl, l1 = Drop[#, 1]& /@ vl, mesh}, mesh = {Drop[l, -1], Drop[l1, -1], Drop[l1, 1], Drop[l, 1]}; Transpose[ Flatten[#, 1]& /@ mesh ] ] (* turn n-gon into list of triangles *) toTriangles[vlist_List] /; Length[vlist]<=2 := Null toTriangles[vlist_List] /; Length[vlist] == 3 := {vlist} (* treat 4gons specially for efficiency *) toTriangles[vlist_List] /; Length[vlist] == 4 := {vlist[[{1, 2, 3}]], vlist[[{3, 4, 1}]]} (* general case: use center of gravity as additional vertex *) (* in this way, some nonconvex polygons can also be rendered correctly *) toTriangles[vlist_List] := Module[{bary = (Plus @@ vlist)/Length[vlist], circ}, circ = Partition[ Append[vlist, First[vlist]], 2, 1 ]; Apply[ {#1, #2, bary}&, circ, {1} ] ] wrimmateSmoothTriangles[file_, {a_, b_, c_, d_}, {na_, nb_, nc_, nd_}] := ( writeSmoothTriangle[file, d, nd, c, nc, a, na]; writeSmoothTriangle[file, b, nb, a, na, c, nc]; ) Clear[writeTriangle] writeTriangle[file_,Null]:=Null writeTriangle[file_, vl_List] /; Length[vl]==3 := ( Write[file, "triangle {"]; writeVector[file, "\t", #]& /@ vl; Write[file, "}"] ) writeSmoothTriangle[file_, p1_, n1_, p2_, n2_, p3_, n3_] := ( Write[file, "smooth_triangle {"]; writeVector[file, "\t ", p1]; Scan[ writeVector[file, "\t, ", #]&, {n1, p2, n2, p3, n3} ]; Write[file, "}"]; ) (* box to world coordinates *) scale[{min_, max_}, rat_] := (min+max)/2 + (max-min)/2 rat boxScale[pr_, pos_] := MapThread[ scale, {pr, pos} ] cam = 3.0 (* approximate field of view *) processOpts[file_, pr_, opts___] := Module[{optlist = Flatten[{opts}], prlow, prhigh, vp, vc, vv, svp, svc, svv, fow, lights,spherical}, Write[file, "/* plot range:"]; writeNums[file, " * ", #]& /@ pr; Write[file, " */"]; {prlow, prhigh} = Transpose[pr]; vp = ViewPoint /. optlist /. Options[Graphics3D]; svp = boxScale[ pr, vp ]; vc = ViewCenter /. optlist /. Options[Graphics3D]; svc = prlow + vc (prhigh - prlow); vv = ViewVertical /. optlist /. Options[Graphics3D]; svv = boxScale[ pr, vv ]; fow = Sqrt[Plus @@ (vp^2)]; (* focal length *) Write[file, "camera {"]; writeVector[file, "\tlocation ", svp]; writeVector[file, "\tdirection ", {0, fow, 0}]; writeVector[file, "\tup ", {0,0,cam}]; writeVector[file, "\tright ", {cam,0,0}]; writeVector[file, "\tsky ", vv]; writeVector[file, "\tlook_at ", svc]; Write[file, "}"]; Write[file,""]; spherical=toSpherical[svp]; Write[file,"/* Sperical camera position */"]; Write[file,"#declare MPi = 3.1415926"]; sp=writeSpherical[file,"Cam",spherical]; Write[file, "/* camera {"]; Write[file, "\tlocation ", sp]; writeVector[file, "\tdirection ", {0, fow, 0}]; writeVector[file, "\tup ", {0,0,cam}]; writeVector[file, "\tright ", {cam,0,0}]; writeVector[file, "\tsky ", vv]; writeVector[file, "\tlook_at ", svc]; Write[file, "} */"]; Write[file,"#declare MathPointSize = 0.1"]; Write[file,"#declare MathPointTexture = texture {pigment {color rgbf<1, 1, 1, 1>}}"]; lights = LightSources /. optlist /. Options[Graphics3D]; i=0; Scan[doLight[file, pr, #,++i]&, lights]; ] toSpherical[{x_,y_,z_}]:= Module[{r,theta,phi}, r=Sqrt[x*x+y*y+z*z]; phi=ArcTan[x,y]; theta=ArcCos[z/r]; Chop /@ {r,theta,phi} ] writeSpherical[file_,prefix_String,{r_,theta_,phi_}]:= Module[{rname,tname,pname}, rname=prefix <> "R"; tname=prefix <> "T"; pname=prefix <> "P"; Write[file,""]; Write[file,"#declare ",rname," = ",r]; Write[file,"#declare ",tname," = ",theta]; Write[file,"#declare ",pname," = ",phi]; Write[file,""]; StringJoin["<",rname,"*cos(",pname,")*sin(",tname,"), ", rname,"*sin(",pname,")*sin(",tname,"), ", rname,"*cos(",tname,")>" ] ] Clear[doLight] doLight[file_, pr_, {pos_, RGBColor[r_, g_, b_]},n_:1] := Module[{scaled,sp}, scaled = boxScale[ pr, pos ]; Write[file, "light_source { " ]; writeVector[file, "\t", scaled]; writeVector[file, "\tcolor rgb ", {r, g, b} ]; Write[file, "}"]; Write[file,""]; Write[file,"/* light in spherical coordinates */"]; sp=toSpherical[scaled]; sp=writeSpherical[file,"Light"<> ToString[n],sp]; Write[file, "/* light_source { " ]; Write[file, "\t", sp]; writeVector[file, "\tcolor rgb ", {r, g, b} ]; Write[file, "} */"]; Write[file,""] ] End[] Protect["POVray`*"] EndPackage[]