Posts Tagged ‘Gaussian Curvature’

This code snippet shown below shows the calculation of Gaussian Curvature on a parametric Bezier Patch. The program flow goes like this –

a) Forming Coon point control grid .
b) Constructing Bezier patch corresponding to the Coon Point Control grid points using De Casteljau’s algorithm.
c) Finding partial derivatives for a given parametric point (u,v) along u direction, v direction, uv direction using De casteljau’s algorithm.
d) Finding normal and double partial derivatives along u and v direction using De Casteljau’s algorithm.
e) Finding Gaussian curvature using the above calculated derivatives and normal.

The following code is  written in Mathematica. The code can be optimized a lot . This will only serve to get a basic understanding of computing Gaussian curvature using De casteljau’s algorithm on Bezier patches.

Output:


(*  Change the values of U and V as  you wish between 0 and 1*)

u = 0.25;
v = 0.25;

(* Boundary Polygons *)
coonPoints = {{{0, 0, 0}, {1, 0, 1}, {2, 0, 1}, {3, 0, 1}, {4, 0, 1}},
   		    {{0, 1, 0}, {}, {}, {}, {4, 1, 3}},
   		    {{0, 2, 0}, {}, {}, {}, {4, 2, 3}},
   		    {{0, 3, 0}, {1, 3, 1}, {2, 3, 1}, {3, 3, 1}, {4, 3, 1}}};

m = 3;
n = 4;
For[i = 1, i < m, i++,
  For[j = 1, j < n, j++,
   temp = ((({{(1 - (i/m))*(coonPoints[[1, j + 1]])}}) + ({{(i/
              m)*(coonPoints[[m + 1, 
              j + 1]])}}) + ({{(1 - (j/n))*(coonPoints[[i + 1, 
              1]])}}) + ({{(j/n)*(coonPoints[[i + 1, n + 1]])}})) - 
      Transpose[{{1 - (j/n)}, {j/n}}].Transpose[{{1 - (i/m), 
           i/m}}.{{coonPoints[[1, 1]], 
           coonPoints[[1, n + 1]]}, {coonPoints[[m + 1, 1]], 
           coonPoints[[m + 1, n + 1]]}}]);
   
   coonPoints[[i + 1, j + 1]] = temp[[1, 1]];]];
       coonPointsbackup = coonPoints;

(* De Cateljau's algorithm function definition *)
dca[c_, r_, i_, t_] := 
  If [r == 0, 
   c[[i + 1]], (1 - t)*dca[c, r - 1, i, t] + 
    t*dca[c, r - 1, i + 1, t]];

dcaPart1[u_, points_] := dca[points, m, 0, u];
dcaPart2[u_, v_, points_] := dca[dcaPart1[u, points], n, 0, v];

dcaPart1upartial[u_, points_] := dca[points, m - 1, 0, u];
dcaPart2upartial[u_, v_, points_] := 
  dca[dcaPart1upartial[u, points], n, 0, v]; 

dcaPart1uupartial[u_, points_] := dca[points, m - 2, 0, u];
dcaPart2uupartial[u_, v_, points_] := 
  dca[dcaPart1uupartial[u, points], n, 0, v]; 

dcaPart1vpartial[u_, points_] := dca[points, m, 0, u];
dcaPart2vpartial[u_, v_, points_] := 
  dca[dcaPart1vpartial[u, points], n - 1, 0, v]; 

dcaPart1vvpartial[u_, points_] := dca[points, m, 0, u];
dcaPart2vvpartial[u_, v_, points_] := 
  dca[dcaPart1vvpartial[u, points], n - 2, 0, v]; 

dcaPart1uvpartial[u_, points_] := dca[points, m - 1, 0, u];
dcaPart2uvpartial[u_, v_, points_] := 
  dca[dcaPart1uvpartial[u, points], n - 1, 0, v]; 

MatrixForm[coonPoints];

(* U partial*)
For[i = 1, i <= m , i++, 
  coonPoints[[i]] = coonPoints[[i + 1]] - coonPoints[[i]]];
a = coonPoints[[1]];
b = coonPoints[[2]];
c = coonPoints[[3]];
upartialmatrix = {a, b, c};
upartialbackup = upartialmatrix;

(* uu partial *)
For[i = 1, i <= m - 1, i++, 
  upartialbackup [[i]] = 
   upartialbackup [[i + 1]] - upartialbackup [[i]]];
aa = upartialbackup [[1]];
bb = upartialbackup [[2]];
uupartialmatrix = {aa, bb};


(* coonPoints again *)

     coonPoints = coonPointsbackup ;

(* V partial*)
For[j = 1, j <= n , j++, 
  coonPoints[[All, j]] = 
   coonPoints[[All, j + 1]] - coonPoints[[All, j]]];
a = coonPoints[[All, 1]];
b = coonPoints[[All, 2]];
c = coonPoints[[All, 3]];
d = coonPoints[[All, 4]];
vpartialmatrix = Transpose[{a, b, c, d}];
vpartialbackup = vpartialmatrix;

(* vv partial *)
For[j = 1, j <= n - 1, j++, 
  vpartialbackup [[All, j]] = 
   vpartialbackup [[All, j + 1]] - vpartialbackup [[All, j]]];
a1 = vpartialbackup [[All, 1]];
a2 = vpartialbackup [[All, 2]];
a3 = vpartialbackup [[All, 3]];
vvpartialmatrix = Transpose[{a1, a2, a3}];


(* uv partial calculation *)

  coonPoints = coonPointsbackup ;

For[i = 1, i <= m , i++, 
  coonPoints[[i]] = coonPoints[[i + 1]] - coonPoints[[i]]];
a = coonPoints[[1]];
b = coonPoints[[2]];
c = coonPoints[[3]];
upartialfinal = {a, b, c};
MatrixForm[upartialfinal];
For[j = 1, j <= n , j++, 
  upartialfinal[[All, j]] = 
   upartialfinal[[All, j + 1]] - upartialfinal[[All, j]]];
MatrixForm[upartialfinal];
a = upartialfinal [[All, 1]];
b = upartialfinal[[All, 2]];
c = upartialfinal [[All, 3]];
d = upartialfinal [[All, 4]];
uvpartialmatrix = Transpose[{a, b, c, d}];


(* Important functions to calculate partials *)
upartialfunction[u_, v_] := m*dcaPart2upartial[u, v, upartialmatrix ];
vpartialfunction[u_, v_] := n*dcaPart2vpartial[u, v, vpartialmatrix ];
uupartialfunction[u_, v_] := 
  m*(m - 1)*dcaPart2uupartial[u, v, uupartialmatrix ];
vvpartialfunction[u_, v_] := 
  n*(n - 1)*dcaPart2vvpartial[u, v, vvpartialmatrix ];
uvpartialfunction[u_, v_] := 
  m*n*dcaPart2uvpartial[u, v, uvpartialmatrix ];
normalfunction[u_, v_] := 
  Cross[upartialfunction[u, v], vpartialfunction[u, v]]/
   Norm[Cross[upartialfunction[u, v], vpartialfunction[u, v]]];


(* Calculating Gaussian Curvature Functions *)
f[upar_, uvpar_, vpar_] := 
  Det[{ {upar.upar, upar.vpar} , {upar.vpar, vpar.vpar }}];
s[norm_, uupar_, uvpar_, vvpar_] := 
  Det[{ {norm.uupar, norm.uvpar} , {norm.uvpar, norm.vvpar}}];
gauss[upar_, vpar_, norm_, uupar_, vvpar_, uvpar_] := 
  s[norm, uupar, uvpar, vvpar]/f[upar, uvpar, vpar];
gaussfinal[u_, v_] := 
  gauss[upartialfunction[u, v], vpartialfunction[u, v], 
   normalfunction[u, v], uupartialfunction[u, v], 
   vvpartialfunction[u, v], uvpartialfunction[u, v]];

(* calculating  gaussian curvature *)
gaussiancurvature = gaussfinal[u, v];
Print["The value of gaussian curvature at " , "(", u, " ," , v, ")", 
  "is", "   ", gaussiancurvature];

coonPoints = coonPointsbackup;
controlGrid = 
  Graphics3D[{PointSize[Medium], Red, Map[Point, coonPoints], Gray, 
    Line[coonPoints], Line[Transpose[coonPoints]]}];

(* Create BezierSurface using DCA *)
bezierSurface = 
  ParametricPlot3D[dcaPart2[u, v, coonPoints], {u, 0, 1}, {v, 0, 1}];




(* calculating actual partial values values here *)

upartial = upartialfunction[u, v];
vpartial = vpartialfunction[u, v];
normal = normalfunction[u, v];
uvpartial = uvpartialfunction[u, v];
uupartial = uupartialfunction[u, v];
vvpartial = vvpartialfunction[u, v];

p = dcaPart2[u, v, coonPoints];
Print[" upartial at " , "(", u, " ," , v, ")", "is", "   ", upartial ];
Print[" vpartial at " , "(", u, " ," , v, ")", "is", "   ", vpartial ];
Print[" normal at " , "(", u, " ," , v, ")", "is", "   ", normal];
Print[" uvpartial at " , "(", u, " ," , v, ")", "is", "   ", 
  uvpartial ];
Print[" uupartial at " , "(", u, " ," , v, ")", "is", "   ", 
  uupartial ];
Print[" vvpartial at " , "(", u, " ," , v, ")", "is", "   ", 
  vvpartial ];
(* Calculating to plot *)
n = p - 0.5*normal;
up = p - 0.25*upartial;
vp = p - 0.25*vpartial;
uvp = p - 0.25*uvpartial;
plotpoints1 = 
  Graphics3D[{RGBColor[1, 0, 0], Cylinder[{p, n}, .03], 
    RGBColor[0, 1, 0], Cylinder[{p, up}, .03], RGBColor[0, 0, 1], 
    Cylinder[{p, vp}, .03], RGBColor[1, 1, 0], 
    Cylinder[{p, uvp}, .03]}];

(* Just Graphics to plot partials and their colors *)
kk = Graphics[{
    
    Text[Style["Normal Vector", Medium, "Label", Darker[Brown]], {1, 
      38}, {-1, 0}],
    Red, Thick, Line[{{1, 36}, {10, 36}}],
    Text[Style["Vpartial Vector", Medium, "Label", Darker[Brown]], {1,
       33}, {-1, 0}],
    Blue, Thick, Line[{{1, 31}, {10, 31}}],
    Text[Style["Upartial Vector", Medium, "Label", Darker[Brown]], {1,
       28}, {-1, 0}],
    Green, Thick, Line[{{1, 26}, {10, 26}}],
    Text[Style["UVpartial Vector", Medium, "Label", 
      Darker[Brown]], {1, 23}, {-1, 0}],
    Yellow, Thick, Line[{{1, 21}, {10, 21}}],
    }, ImageSize -> {150, 300}, PlotRange -> All];


Show[kk, PlotRange -> All]
Show[plotpoints1, controlGrid, bezierSurface, 
 PlotLabel -> 
  Style["upartial, vpartial, normal, uvpartial vectors plotted", 16],
 PlotRange -> All, Axes -> True]