任意画像の輪郭を数式に変換してプロットする (Mathematica ver.8)
(* 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}]] | |
*) |
This comment has been minimized.
This comment has been minimized.
|
This comment has been minimized.
This comment has been minimized.
お借りさせていただきました。 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This comment has been minimized.
RaspberryPi+Mathematicaでも、何十分かの処理でできましたです