(*  The following are examples of curvature analysis  on  individual tubercles in seeds of Silene species*)






	(*******************************  S. ciliata 1 to 6 *)

Needs["Splines`"]

ciliata[1]={{0.5000,9.0000},{2.5000,8.0000},{4.5000,7.5000},{6.5000,6.5000},{8.5000,7.0000},{10.5000,8.0000},{12.5000,10.0000},{14.5000,13.0000},{16.5000,16.0000},{18.5000,18.5000},{20.5000,21.5000},{22.5000,25.0000},{24.5000,29.0000},{26.5000,32.5000},{28.5000,36.5000},{30.5000,40.0000},{32.5000,43.5000},{34.5000,46.0000},{36.5000,48.0000},{38.5000,50.0000},{40.5000,51.5000},{42.5000,52.5000},{44.5000,53.5000},{46.5000,54.5000},{48.5000,54.5000},{50.5000,54.5000},{52.5000,54.5000},{54.5000,54.5000},{56.5000,54.0000},{58.5000,53.0000},{60.5000,51.5000},{62.5000,49.0000},{64.5000,46.5000},{66.5000,44.0000},{68.5000,41.0000},{70.5000,38.0000},{72.5000,34.5000},{74.5000,31.0000},{76.5000,27.0000},{78.5000,22.5000},{80.5000,19.5000},{82.5000,16.0000},{84.5000,14.0000},{86.5000,13.0000},{88.5000,11.0000},{90.5000,10.0000},{92.5000,8.0000},{94.5000,6.5000}}
ciliata[2]={{4.5000,1.0000},{6.5000,3.0000},{8.5000,5.0000},{10.5000,9.0000},{12.5000,15.0000},{14.5000,20.0000},{16.5000,24.5000},{18.5000,28.5000},{20.5000,32.0000},{22.5000,35.0000},{24.5000,36.0000},{26.5000,38.0000},{28.5000,38.5000},{30.5000,39.5000},{32.5000,39.5000},{34.5000,40.5000},{36.5000,40.5000},{38.5000,39.5000},{40.5000,39.0000},{42.5000,38.5000},{44.5000,37.5000},{46.5000,36.0000},{48.5000,35.0000},{50.5000,33.0000},{52.5000,30.0000},{54.5000,26.0000},{56.5000,20.0000},{58.5000,13.5000},{60.5000,9.5000},{62.5000,7.0000},{64.5000,5.5000},{66.5000,4.5000},{68.5000,4.0000}}
ciliata[3]={{0.5000,7.0000},{2.5000,7.5000},{4.5000,9.0000},{6.5000,10.0000},{8.5000,11.5000},{10.5000,12.5000},{12.5000,15.0000},{14.5000,18.0000},{16.5000,19.5000},{18.5000,22.0000},{20.5000,24.0000},{22.5000,26.0000},{24.5000,28.0000},{26.5000,30.0000},{28.5000,32.0000},{30.5000,34.0000},{32.5000,35.0000},{34.5000,35.5000},{36.5000,36.5000},{38.5000,37.5000},{40.5000,38.0000},{42.5000,38.5000},{44.5000,39.0000},{46.5000,38.5000},{48.5000,38.5000},{50.5000,38.5000},{52.5000,38.0000},{54.5000,37.5000},{56.5000,36.5000},{58.5000,35.5000},{60.5000,34.0000},{62.5000,33.0000},{64.5000,31.0000},{66.5000,29.5000},{68.5000,28.0000},{70.5000,25.0000},{72.5000,22.5000},{74.5000,20.0000},{76.5000,17.0000},{78.5000,13.0000},{80.5000,9.0000},{82.5000,6.0000},{84.5000,2.5000},{86.5000,0.5000}}
ciliata[4]={{6.5000,1.0000},{8.5000,4.0000},{10.5000,8.5000},{12.5000,13.5000},{14.5000,18.5000},{16.5000,21.5000},{18.5000,25.0000},{20.5000,28.0000},{22.5000,30.0000},{24.5000,32.5000},{26.5000,34.5000},{28.5000,37.0000},{30.5000,39.0000},{32.5000,40.5000},{34.5000,42.0000},{36.5000,43.0000},{38.5000,44.5000},{40.5000,45.5000},{42.5000,45.5000},{44.5000,45.5000},{46.5000,45.5000},{48.5000,44.5000},{50.5000,44.0000},{52.5000,43.0000},{54.5000,42.5000},{56.5000,40.5000},{58.5000,39.0000},{60.5000,37.0000},{62.5000,34.5000},{64.5000,32.5000},{66.5000,28.5000},{68.5000,24.0000},{70.5000,19.5000},{72.5000,15.5000},{74.5000,11.5000},{76.5000,9.0000},{78.5000,7.0000},{80.5000,6.0000},{82.5000,5.5000},{84.5000,4.5000},{86.5000,3.5000}}
ciliata[5]={{0.500,24.5000},{2.500,24.5000},{4.500,25.5000},{6.500,25.5000},{8.500,26.5000},{10.500,27.0000},{12.500,29.0000},{14.500,31.0000},{16.500,33.0000},{18.500,35.5000},{20.500,39.5000},{22.500,43.5000},{24.500,48.5000},{26.500,53.5000},{28.500,57.5000},{30.500,60.5000},{32.500,64.0000},{34.500,67.0000},{36.500,68.0000},{38.500,70.0000},{40.500,71.0000},{42.500,72.0000},{44.500,73.0000},{46.500,73.5000},{48.500,75.0000},{50.500,75.5000},{52.500,76.0000},{54.500,76.5000},{56.500,76.5000},{58.500,76.5000},{60.500,75.5000},{62.500,75.5000},{64.500,74.5000},{66.500,73.5000},{68.500,73.0000},{70.500,72.0000},{72.500,71.5000},{74.500,70.5000},{76.500,69.5000},{78.500,68.5000},{80.500,67.0000},{82.500,64.5000},{84.500,62.0000},{86.500,58.5000},{88.500,55.0000},{90.500,51.0000},{92.500,49.0000},{94.500,46.0000},{96.500,43.5000},{98.500,41.0000},{100.500,39.0000},{102.500,38.0000},{104.500,36.0000},{106.500,35.0000},{108.500,33.5000},{110.500,32.5000},{112.500,31.5000}}
ciliata[6]={{0.5000,0.5000},{2.5000,2.0000},{4.5000,4.0000},{6.5000,6.0000},{8.5000,7.0000},{10.5000,9.5000},{12.5000,13.5000},{14.5000,19.5000},{16.5000,24.5000},{18.5000,28.5000},{20.5000,33.0000},{22.5000,37.5000},{24.5000,41.5000},{26.5000,45.0000},{28.5000,47.0000},{30.5000,49.0000},{32.5000,50.5000},{34.5000,52.0000},{36.5000,53.0000},{38.5000,54.0000},{40.5000,55.0000},{42.5000,56.0000},{44.5000,56.5000},{46.5000,56.5000},{48.5000,56.5000},{50.5000,57.5000},{52.5000,57.5000},{54.5000,56.5000},{56.5000,54.0000},{58.5000,53.0000},{60.5000,50.0000},{62.5000,47.5000},{64.5000,42.5000},{66.5000,36.0000},{68.5000,31.5000},{70.5000,29.0000},{72.5000,26.0000},{74.5000,21.5000},{76.5000,17.0000},{78.5000,14.0000},{80.5000,10.5000},{82.5000,7.0000},{84.5000,5.5000},{86.5000,3.5000},{88.5000,2.5000},{90.5000,1.0000},{92.5000,0.5000},{94.5000,0.5000}}

(*Poner aqui la familia*)
var=ciliata;
(* el número de curvas estudiadas/semillas de la familia*)
m=6;
For[j=1,j≤m,j++,
puntos=var[j];
(* pru es el spline cubico
d es el vector que contiene la información de pru 
n es el número de trozos del spline, es decir, el número de puntos menos 1
*)
pru=SplineFit[puntos,Cubic];
d=InputForm[pru];
n=IntegerPart[d[[1,2,2]]];
spli=ParametricPlot[pru[t],{t,0,n},PlotRange->All, DisplayFunction->Identity];
(* Con los coeficientes se calculan los polinomios cúbicos que dan el spline. Habrá n polinomios para la primera coordenada (X) y otros tantos para la segunda
*)
For[i=0,i<n,i++,
SX[i,t_]=d[[1,4,i+1,1,1]]+d[[1,4,i+1,1,2]]*(t-i)+d[[1,4,i+1,1,3]]*(t-i)^2+d[[1,4,i+1,1,4]]*(t-i)^3];
For[i=0,i<n,i++,SY[i,t_]=d[[1,4,i+1,2,1]]+d[[1,4,i+1,2,2]]*(t-i)+d[[1,4,i+1,2,3]]*(t-i)^2+d[[1,4,i+1,2,4]]*(t-i)^3];
(* Derivada primera y segunda de cada trozo del spline
*)
For[i=0,i<n,i++,DX[i,t_]=d[[1,4,i+1,1,2]]+2*d[[1,4,i+1,1,3]]*(t-i)+3*d[[1,4,i+1,1,4]]*(t-i)^2];
For[i=0,i<n,i++,DY[i,t_]=d[[1,4,i+1,2,2]]+2*d[[1,4,i+1,2,3]]*(t-i)+3*d[[1,4,i+1,2,4]]*(t-i)^2];
For[i=0,i<n,i++,DDX[i,t_]=2*d[[1,4,i+1,1,3]]+6*d[[1,4,i+1,1,4]]*(t-i)];
For[i=0,i<n,i++,DDY[i,t_]=2*d[[1,4,i+1,2,3]]+6*d[[1,4,i+1,2,4]]*(t-i)];
(*Cálculo de la curvatura a trozos
*)
For[i=0,i<n,i++,cur[i,t_]=(DX[i,t]DDY[i,t]-DDX[i,t]DY[i,t])/((DX[i,t])^2+(DY[i,t])^2)^(3/2)];
(* Dibujos de los trozos de curvatura 
*)
For[i=0,i<n,i++,dib[i]=Plot[cur[i,t],{t,i,i+1},DisplayFunction->Identity]];
pp=dib[0];
For[i=1,i<n,i++,pp=Show[pp,dib[i],PlotRange->All]];
(* bez es el spline de Bezier
b es el vector que contiene la información de bez
nb es el número de trozos del spline, es decir, el número de puntos menos 1 (Como nb=n, no se define nuevo)
*)
bez=SplineFit[puntos,Bezier];
b=InputForm[bez];
(*Calculo de las funciones de Bezier para la coordenada X y la Y; la curva es la imagen del intervalo[0,1]*)
For[i=1,i≤n+1,i++,BX[t_]=bez[[4,1,1]]/.#1->t];
For[i=1,i≤n+1,i++,BY[t_]=bez[[4,1,2]]/.#1->t];
dibez[j]=ParametricPlot[{BX[t],BY[t]},{t,0,1},PlotRange->All,DisplayFunction->Identity];
curbe[j,t_]=(BX'[t]BY''[t]-BX''[t]BY'[t])/((BX'[t])^2+(BY'[t])^2)^(3/2);
cur[j]=Abs[NMinValue[{curbe[j,t],0.30<t<0.7},t]];
vmed[j]=Abs[NIntegrate[curbe[j,t],{t,0.3,0.7}]/0.4];
dibcurbe[j]=Plot[{curbe[j,t],-vmed[j]},{t,0,1},PlotRange->{-0.2,0.1},DisplayFunction->Identity,PlotStyle->{RGBColor[0,1,0],{Dotted,Dashed}}];
total=GraphicsGrid[{
(*{spli,pp},*)
{dibez[j],dibcurbe[j]}},Spacings->0.3];
Print[var,j, ":", "Máximo valor de la curvatura (en valor absoluto)",cur[j]];
Print[total]
]
Print["Curvaturas de cada ", var ,":  ", Table[cur[j],{j,1,6}]];
Print["Valores medios de curvatura de cada ", var ," en [0.3,0.7]:  ", Table[vmed[j],{j,1,6}]];

Show[Table[dibcurbe[j],{j,1,6}]]

Table[NMinValue[{curbe[j,t],0.25<t<0.75},t],{j,1,6}]

dibcurbe[j],{j,1,6}

Table[dibcurbe[j],{j,1,6}]

curbe[j,t_]
dibcurbe[j]