任意画像の輪郭を数式に変換してプロットする (Mathematica ver.8)

## gistfile1.mathematica
(* parameters *)
(* 最大次数 *)
maxOrderNum = 200;
(* 画像URL, もしくはローカルパス *)
imageURL = "http://nex.fm/wp-content/uploads/2012/08/vim-editor_logo.png";

pointListToLines[pointList_, neighbothoodSize_: 6] :=
 Module[{L = DeleteDuplicates[pointList], NF, lambda,
   lineBag, counter, seenQ, sLB, nearest,
   nearest1, nextPoint, couldReverseQ, d, n, s},
  NF = Nearest[L];
  lambda = Length[L];
  Monitor[
   (*list of segments *)
   lineBag = {};
   counter = 0;
   While[counter < lambda,
    (*new segment*)
    sLB = {RandomChoice[DeleteCases[L, _?seenQ]]};
    seenQ[sLB[[1]]] = True;
    counter++;
    couldReverseQ = True;
    (*complete segment*)
    While[
     (nearest = NF[Last[sLB], {Infinity, neighbothoodSize}];
      nearest1 =
       SortBy[DeleteCases[nearest, _?seenQ],
        1. EuclideanDistance[Last[sLB], #] &];
      nearest1 =!= {} || couldReverseQ),
     If[
      nearest1 === {},
      (*extend the other end; penalize sharp edges*)
      sLB = Reverse[sLB];
      couldReverseQ = False,

      (* prefer straight continuation *)
      nextPoint = If[Length[sLB] <= 3,
        nearest1[[1]],
        d = 1. Normalize[(sLB[[-1]] - sLB[[-2]]) +
            1/2 (sLB[[-2]] - sLB[[-3]])];
        n = {-1, 1} Reverse[d];
        s = Sort[{Sqrt[(d.(# - sLB[[-1]]))^2 +
               (*perpendicular *)2
                (n. (# - sLB[[-1]]))^2], #} & /@ nearest1];
        s[[1, 2]]];
      AppendTo[sLB, nextPoint];
      seenQ[nextPoint] = True;
      counter++]];
    AppendTo[lineBag, sLB]];
   (*return segments sorted by length*)
   Reverse[SortBy[Select[lineBag, Length[#] > 12 &],
     Length]],
   (*monitor progress*)
   Grid[
    {{Text[Style["progress point joining",
        Darker[Green, 0.66]]],
      ProgressIndicator[counter/lambda]},
     {Text[Style["number of segments",
        Darker[Green, 0.66]]],
      Length[lineBag] + 1}},
    Alignment -> Left, Dividers -> Center]]]

(* Fourier coefficients of a single curve *)
fourierComponentData[pointList_, nMax_, op_] :=
 Module[{epsilon = 10^-3, myu = 2^14, M = 10000, s, scale, delta, L,
   nds, sMax, if, xyFunction, X, Y, XFT, YFT, type},
  (* prepare curve *)
  scale =
   1. Mean[Table[
      Max[fl /@ pointList] -
       Min[fl /@ pointList], {fl, {First, Last}}]];
  delta = EuclideanDistance[First[pointList], Last[pointList]];
  L = Which[
    op === "Closed",
    type = "Closed";
    If[First[pointList] === Last[pointList], pointList,
     Append[pointList, First[pointList]]],

    op === "Open",
    type = "Open";
    pointList,

    delta == 0.,
    type = "Closed";
    pointList,

    delta/scale < op,
    type = "Closed";
    Append[pointList, First[pointList]],

    True,
    type = "Open";
    Join[pointList, Rest[Reverse[pointList]]]];
  (*re-parametrize curve by arclength *)
  xyFunction = BSplineFunction[L, SplineDegree -> 4];
  nds = NDSolve[
    {s'[t] == Sqrt[xyFunction'[t].xyFunction'[t]],
     s[0] == 0}, s, {t, 0, 1}, MaxSteps -> 10^5, PrecisionGoal -> 4];
  (* total curve length *)
  sMax = s[1] /. nds[[1]];
  if = Interpolation[
    Table[{s[rho] /. nds[[1]], rho}, {rho, 0, 1, 1/M}]];
  X[t_Real] :=
   BSplineFunction[L][Max[Min[1, if[(t + Pi)/(2 Pi) sMax]], 0]][[1]];
  Y[t_Real] :=
   BSplineFunction[L][Max[Min[1, if[(t + Pi)/(2 Pi) sMax]], 0]][[2]];
  (* extract Fourier coefficients *)
  {XFT, YFT} =
   Fourier[Table[#[N@t], {t, -Pi + epsilon,
        Pi - epsilon, (2 Pi - 2 epsilon)/myu}]] & /@ {X, Y};
  {type, 2 Pi/
     Sqrt[myu]*((Transpose[
         Table[{Re[#], Im[#]} &[Exp[I k Pi] #[[k + 1]]], {k, 0,
           nMax}]] & /@ {XFT, YFT}))}]
Options[fourierComponents] =
  {"MaxOrder" -> maxOrderNum, "OpenClose" -> 0.025};

fourierComponents[pointLists_, OptionsPattern[]] :=
 Monitor[
   Table[fourierComponentData[
     pointLists[[k]],
     If[Head[#] === List, #[[k]], #] &[OptionValue["MaxOrder"]],
     If[Head[#] === List, #[[k]], #] &[OptionValue["OpenClose"]]
     ], {k, Length[pointLists]}],
   Grid[
    {{Text[
       Style[
        "progress calculating Fourier coefficients",
        Darker[Green, 0.66]]],
      ProgressIndicator[k/Length[pointLists]]}},
    Alignment -> Left, Dividers -> Center]] /; Depth[pointLists] === 4

makeFourierSeries[
  {"Closed" | "Open", {{cax_, sax_}, {cay_, say_}}},
  t_, n_] :=
 {Sum[If[k == 0, 1/2, 1] cax[[k + 1]] Cos[k t] +
    sax[[k + 1]] Sin[k t], {k, 0, Min[n, Length[cax]]}],
  Sum[If[k == 0, 1/2, 1] cay[[k + 1]] Cos[k t] +
    say[[k + 1]] Sin[k t], {k, 0, Min[n, Length[cay]]}]}

paraplot[n_] :=
 Show[
  {ParametricPlot[
    Evaluate[makeFourierSeries[#, t, n] & /@ fCs],
    {t, -Pi, Pi}, Axes -> False](*,
   Graphics[Text[Style["n="<>ToString[n],Large,Bold],
   Scaled[{.9,.1}] ]]*)}(*,ImageSize->{500,500}*)]

(* 画像読み込み *)
img = Import[imageURL];

(* エッジ抽出 *)
edgeImage = Thinning[EdgeDetect[ColorConvert[
     ImagePad[Image[Map[Most, ImageData[img], {2}]], 20, White],
     "Grayscale"]]];
edgePoints = {#2, -#1} & @@@ Position[ImageData[edgeImage], 1, {2}];
SeedRandom[2];
hLines = pointListToLines[edgePoints, 16];

(* 線分の数を表示 *)
(* Length[hLines] *)
(* エッジ抽出\[RightArrow]線分に変換した結果を描画 *)
(* Graphics[{ColorData["DarkRainbow"][RandomReal[]],Line[#]}&/@hLines]\
 *)


(* 変換された数式を変換してプロットする *)
(* paraplotにmaxOrderNum以下の数字をわたすことで n次方程式時点での描画をする *)
fCs = fourierComponents[hLines];
paraplot[maxOrderNum]

(* 1～maxOrderNum次までの描画を全て並べてgifに出力 *)
(*
Dynamic[n]
Export["plot.gif",
 Table[(n = x; paraplot[x]), {x, 1, maxOrderNum, 1}]]
*)