MathGroup Archive 1996

[Date Index] [Thread Index] [Author Index]

Search the Archive

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[]



  • Prev by Date: Re: [Q] How to read (and analyze) .WAV files?
  • Next by Date: Re: arguments [Manipulation of FindRoot variable lists]
  • Previous by thread: Re: [Q] How to read (and analyze) .WAV files?
  • Next by thread: Re: arguments [Manipulation of FindRoot variable lists]