-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSnails for Mathematica 1.0.m
203 lines (154 loc) · 9.85 KB
/
Snails for Mathematica 1.0.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
(* ::Package:: *)
(* ::Title:: *)
(*Snail functions for Mathematica*)
(* ::Text:: *)
(*Version 1.0*)
(*12 July 2022*)
(**)
(*This notebook contains functions for simulating and processing virtual gastropod (and other mollusc) shells. *)
(**)
(*P. David Polly*)
(*Earth & Atmospheric Sciences*)
(*Indiana University*)
(*pdpolly@indiana.edu*)
(*https://pollylab.indiana.edu/ *)
(* This function prints the version number for the installed verison of this package. *)
SnailsVersion[]:=Print["Snails 1.0\n(c) P. David Polly, 12 July 2022\n"];
SnailsVersion[];
(* This function generates a three-dimensional rendered shell using the functions
developed by David Raup (1966).
Usage: RaupCoil3D[shape, W, T, D, TurnNum]
where shape is a matrix of 3D semilandmarks that circumscribe the aperture
W is Raup's whorl expansion rate parameter
T is Raup's vertical translation rate parameter
D is Raup's parameter for distance of aperture from columella
TurnNum is the desired number of whorls
Created 16 March 2016.
*)
RaupCoil3D[S_,W_,T_,Di_, turns_,CoordsOnly_:False]:=Quiet[Block[{r,z,x,CoiledCoords,\[Theta],Sprime,S3D,yaspect,rc},
Sprime=#-Mean[S]&/@S;
Sprime=(Sprime/(Max[Sprime[[1;;,1]]]-Min[Sprime[[1;;,1]]]));
Sprime[[1;;,1]]=Sprime[[1;;,1]]-Min[Sprime[[1;;,1]]]+Di-0.3;
Sprime[[1;;,3]]=Sprime[[1;;,3]]+0.33*(Max[Sprime[[1;;,3]]]-Min[Sprime[[1;;,3]]]);
yaspect=Max[Sprime[[1;;,3]]]-Min[Sprime[[1;;,3]]];
rc=Mean[Sprime[[1;;,1]]];
S3D=Sprime;
CoiledCoords=Table[Table[
r=S3D[[x,1]]*(W^(\[Theta]/(2Pi)));
z=S3D[[x,3]]*(W^(\[Theta]/(2Pi)))+rc*T*yaspect*((W^(\[Theta]/(2Pi)))-1);
{r*Cos[\[Theta]],r*Sin[\[Theta]]+S3D[[x,2]]*(W^(\[Theta]/(2Pi))),z}
,{x,Length[S3D]}]
,{\[Theta],0,turns*2Pi,5Degree}];
If[CoordsOnly==True,Return[CoiledCoords],
Return[Graphics3D[Table[{Table[{RGBColor[0.997681, 0.980789, 0.842313],EdgeForm[None],Opacity[1],Polygon[{CoiledCoords[[x,y]],CoiledCoords[[x-1,y]],CoiledCoords[[x-1,y+1]],CoiledCoords[[x,y+1]]}]},{y,Length[CoiledCoords[[x]]]-1}],{GrayLevel[0.5],Opacity[0],AbsolutePointSize[1.75],Point[CoiledCoords[[x]]]}},{x,2,Length[CoiledCoords]}],Boxed->False]]];
]];
(* This function generates a three-dimensional rendered shell using the functions
developed by David Raup (1966) as a three-dimensional mesh suitable for printing
with a 3D preinter. It is nearly idential ato RaupCoil3D execpt it creates a shell
with two layers that is capped at apex and aperture.
Usage: RaupCoil3DForPrint[shape, W, T, D, TurnNum]
where shape is a matrix of 3D semilandmarks that circumscribe the aperture
W is Raup's whorl expansion rate parameter
T is Raup's vertical translation rate parameter
D is Raup's parameter for distance of aperture from columella
TurnNum is the desired number of whorls
Created 16 March 2016.
*)
RaupCoil3DForPrint[S_,W_,T_,Di_, turns_]:=Quiet[Module[{r,z,x,CoiledCoords,\[Theta],Sprime,S3D,yaspect,rc},
Sprime=#-Mean[S]&/@S;
Sprime=Partition[Flatten[{Sprime,Sprime*0.85}],3];
Sprime=(Sprime/(Max[Sprime[[1;;,1]]]-Min[Sprime[[1;;,1]]]));
Sprime[[1;;,1]]=Sprime[[1;;,1]]-Min[Sprime[[1;;,1]]]+Di;
Sprime[[1;;,3]]=Sprime[[1;;,3]]+.33*(Max[Sprime[[1;;,3]]]-Min[Sprime[[1;;,3]]]);
yaspect=Max[Sprime[[1;;,3]]]-Min[Sprime[[1;;,3]]];
rc=Mean[Sprime[[1;;,1]]];
S3D=Sprime;
CoiledCoords=Table[Table[
r=S3D[[x,1]]*(W^(\[Theta]/(2Pi)));
z=S3D[[x,3]]*(W^(\[Theta]/(2Pi)))+rc*T*yaspect*((W^(\[Theta]/(2Pi)))-1);
{r*Cos[\[Theta]],r*Sin[\[Theta]]+S3D[[x,2]]*(W^(\[Theta]/(2Pi))),z}
,{x,Length[S3D]}]
,{\[Theta],0,turns*2Pi,10Degree}];
Return[CoiledCoords];
]];
(* This function tilts the aperture shape at a specified number of degrees.
Usage: TiltAperture[shape, angle]
where shape is a matrix of 3D semilandmarks that circumscribe the aperture
and angle is an angle from vertical given in degrees
Created 16 March 2016.
*)
TiltAperture[aperture_,angle_]:=Module[{tiltedaperture},
tiltedaperture=aperture . RotationMatrix[angle Degree,{1,0,0}];
tiltedaperture[[1;;,2]]=tiltedaperture[[1;;,2]]+Mean[tiltedaperture[[1;;,2]]];
Return[tiltedaperture];
];
(* This function generates three-dimensional elliptical Fourier coefficiences from
a closed curve (like an aperture shape). It returns Fourier cofficients for a
single shape.
Usage: EllipticalFourierCoefficients[points]
where points is a matrix of 3D semilandmarks that circumscribe the aperture
Created 16 March 2016.
*)
EllipticalFourierCoefficients[standardizedoutline_]:=Module[{deltaXY,deltaT,TotalT,FourierCoeffs},
deltaXY=Table[standardizedoutline[[k+1]]-standardizedoutline[[k]],{k,Length[standardizedoutline]-1}];
deltaT=Sqrt[Plus@@(#^2)]&/@deltaXY;
TotalT=Plus@@deltaT;
FourierCoeffs=Table[{TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,1]]/deltaT[[k]]*(Cos[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Cos[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}]),TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,1]]/deltaT[[k]]*(Sin[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Sin[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}]),TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,2]]/deltaT[[k]]*(Cos[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Cos[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}]),TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,2]]/deltaT[[k]]*(Sin[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Sin[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}])},{n,Length[deltaXY]/2}];
Return[{FourierCoeffs,deltaT,TotalT}];
];
(* This function generates a shape in X, Y coordinates from elliptical Fourier coefficients.
Usage: HarmonicReconstruction[coefficients, deltaT, T, n]
where coefficients is a matrix of Fourier coefficients
deltaT is the offset
T is the number
n is
Created 16 March 2016.
*)
HarmonicReconstruction[coeffs_,deltaT_,T_,n_]:=Module[{xy},
xy=Table[Flatten[{Plus@@Table[coeffs[[j,1]]*Cos[(2Pi*j*Plus@@deltaT[[1;;i]])/T]+coeffs[[j,2]]*Sin[(2Pi*j*Plus@@deltaT[[1;;i]])/T],{j,n}],Plus@@Table[coeffs[[j,3]]*Cos[(2Pi*j*Plus@@deltaT[[1;;i]])/T]+coeffs[[j,4]]*Sin[(2Pi*j*Plus@@deltaT[[1;;i]])/T],{j,n}]}],{i,Length[deltaT]}];
Return[xy]
];
(* This function saves a 3D shell suitable for printing as an OBJ mesh file.
Usage: SaveShell[shellobject, parameters, title, path]
where shellobject is a shell generated with the RaupCoil3DForPrint function
parameters are the shell coiling parameters used to generate the shell(which are included as metadata in the OBJ file)
title is a title for the shell object (included as metadata in the OBJ file)
path is the path to the location where the file will be saved.
Created 16 March 2016.
*)
SaveShell[myShell_,parameters_,title_,path_]:=Module[{s,vertexnumbers,verts,outer,inner,outerfaces,innerfaces,endfaces,innerprotoconchcenter,innerprotofaces, outerprotofaces,outerprotoconchcenter,W,T,shp,turn,Di},
{W,T,Di,turn}=parameters;
(* calculating the faces *)
vertexnumbers=Partition[Table[x,{x,Dimensions[myShell][[1]]*Dimensions[myShell][[2]]}],Dimensions[myShell][[2]]];
verts=StringRiffle[Prepend[#,"v"]]&/@Flatten[myShell,1];
outer={1,Dimensions[vertexnumbers][[2]]/2};
inner={(Dimensions[vertexnumbers][[2]]/2)+1,Dimensions[vertexnumbers][[2]]};
outerfaces=Flatten[{Table[Table[{StringRiffle[Prepend[
{vertexnumbers[[i+1,j]],vertexnumbers[[i,j+1]],vertexnumbers[[i,j]]},"f"]],StringRiffle[Prepend[{vertexnumbers[[i+1,j+1]],vertexnumbers[[i,j+1]],vertexnumbers[[i+1,j]]},"f"]]},{j,outer[[1]],outer[[2]]-1}],{i,Dimensions[vertexnumbers][[1]]-1}],Table[StringRiffle[Prepend[{vertexnumbers[[i,outer[[1]]]],vertexnumbers[[i,outer[[2]]]],vertexnumbers[[i+1,outer[[2]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}],
Table[StringRiffle[Prepend[{vertexnumbers[[i,outer[[1]]]],vertexnumbers[[i+1,outer[[2]]]],vertexnumbers[[i+1,outer[[1]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}]}];
innerfaces=Flatten[{Table[Table[{StringRiffle[Prepend[
{vertexnumbers[[i,j]],vertexnumbers[[i,j+1]],vertexnumbers[[i+1,j]]},"f"]],StringRiffle[Prepend[{vertexnumbers[[i+1,j]],vertexnumbers[[i,j+1]],vertexnumbers[[i+1,j+1]]},"f"]]},{j,inner[[1]],inner[[2]]-1}],{i,Dimensions[vertexnumbers][[1]]-1}],Table[StringRiffle[Prepend[{vertexnumbers[[i+1,inner[[2]]]],vertexnumbers[[i,inner[[2]]]],vertexnumbers[[i,inner[[1]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}],
Table[StringRiffle[Prepend[{vertexnumbers[[i+1,inner[[1]]]],vertexnumbers[[i+1,inner[[2]]]],vertexnumbers[[i,inner[[1]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}]}];
endfaces=Flatten[{Table[{StringRiffle[{"f",vertexnumbers[[-1,i+1]],vertexnumbers[[-1,i]],vertexnumbers[[-1,i+inner[[1]]-1]]}],StringRiffle[{"f",vertexnumbers[[-1,i+inner[[1]]]],vertexnumbers[[-1,i+1]],vertexnumbers[[-1,i+inner[[1]]-1]]}]},{i,outer[[2]]-1}],
StringRiffle[{"f",vertexnumbers[[-1,outer[[1]]]],vertexnumbers[[-1,outer[[2]]]],vertexnumbers[[-1,inner[[2]]]]}],
StringRiffle[{"f",vertexnumbers[[-1,outer[[1]]]],vertexnumbers[[-1,inner[[2]]]],vertexnumbers[[-1,inner[[1]]]]}]}];
innerprotoconchcenter=Mean[myShell[[1,inner[[1]];;inner[[2]]]]];
verts=Append[verts,StringRiffle[Prepend[innerprotoconchcenter,"v"]]];
innerprotofaces=Append[Table[StringRiffle[{"f",Length[verts],i+1,i}],{i,inner[[1]],inner[[2]]-1,1}],StringRiffle[{"f",Length[verts],inner[[1]],inner[[2]]}]];
outerprotoconchcenter=innerprotoconchcenter+{0,-1*Norm[innerprotoconchcenter-myShell[[1,1]]]/2,0};
verts=Append[verts,StringRiffle[Prepend[outerprotoconchcenter,"v"]]];
outerprotofaces=Append[Table[StringRiffle[{"f",Length[verts],i,i+1}],{i,outer[[1]],outer[[2]]-1,1}],StringRiffle[{"f",Length[verts],outer[[2]],outer[[1]]}]];
(* Writing the data *)
s=OpenWrite[path<>title<>".obj"];
WriteLine[s,"# Shell Mesh created by SaveShell[], P. David Polly "<>DateString[]];
WriteLine[s, "# Version 2.1"];
WriteLine[s,"# Parameters: W="<>ToString[W]<>" T="<>ToString[T]<>" D="<>ToString[Di]<>" turns="<>ToString[turn]];
WriteLine[s,#]&/@verts;
WriteLine[s,#]&/@outerfaces;
WriteLine[s,#]&/@innerfaces;
WriteLine[s,#]&/@endfaces;
WriteLine[s,#]&/@innerprotofaces;
WriteLine[s,#]&/@outerprotofaces;
WriteLine[s,"# End of File \n"];
Close[s]
];