２次元点群（point cloud）からパーシステント図を計算する

方針：
- 点群の要素数 N に対して，可能であればO(N)のアルゴリズムを目指す．
- 点群からドロネー図の構成はMathematicaの関数を使う．
- lineは{i,j}, faceは{i,j,k}のように頂点のラベルで記録されている．これをソートしておく．
- 境界作用素{line}->{vertex}はlineのラベル{i,j}そのもの．
- faceは{i,j,k} -> {i,j},{i,k},{j,k} として境界のlineが取り出せる．平坦化(Flatten)をかけて再びソートすると，境界作用素{face}->{line}が得られる．この際，Orderingを使えばソート前後の置換が得られる．
- セルの出現順は以下：
	- vertexは時刻0，
	- lineは，lineの長さ/2, 
	- faceは，face（三角形）の外接円の半径
- 出現順にセルをソートすると，境界作用素: {cell}->{cell}が得られる．
- ここまではO(N)で実行される．N=10^5くらいであれば数秒で計算が終わる．

- パーシステント図を得るアルゴリズムとして，制限された掃き出し法を使う．アルゴリズムは以下

大林一平，「位相的データ解析の現在」
https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2057-05.pdf
より引用．
- 制限掃き出し法の計算量は，セルのラベルjを走り，さらにwhileの条件が満たされなくなるまで，iのラベルがi<jの範囲内で走るため，最悪の場合でO(N^3)程度となるらしい[Nina Otter, Mason A. Porter, Ulrike Tillmann, Peter Grindrod, Heather A. Harrington, “A roadmap for the computation of persistent homology”, arXiv:1506.08903]．
- よって制限掃き出し法の計算がボトルネックとなる．この点の改善を調べる必要がある．
- ラップトップだとN=500程度であれば数秒で計算が実行される．

(*画像の準備．画像からエッジ（色勾配が大きな点）を抽出する．以下のweb pageを参考にした: https://aomoriringo.hateblo.jp/entry/2013/11/30/074758 *)
imageURL(*任意の画像URL*)="https://scholar.googleusercontent.com/citations?view_op=medium_photo&user=mbTEhB4AAAAJ&citpid=9"(*これは自分のgoogle scholarプロフィール画像*);
(*imageURL(*任意の画像URL*)="https://upload.wikimedia.org/wikipedia/commons/thumb/a/a7/Olympic_flag.svg/1280px-Olympic_flag.svg.png"(*これはオリンピックのロゴ*);*)
img=Import[imageURL]
edgeImage=Thinning[EdgeDetect[ColorConvert[ImagePad[Image[Map[Most,ImageData[img],{2}]],20,White],"Grayscale"]]]
edgeimg=ImageData[edgeImage];
pts(*画像のエッジ上の点群*)=Position[ImageData[edgeImage],1];
Nsample(*点群のサンプル数*)=500;
Spts(*Nsampleだけサンプルする．排他的サンプリング*)=RandomSample[pts,Nsample];
matimg=SparseArray[Table[x->1,{x,Spts}],Dimensions[edgeimg],0]//Normal;
Image[matimg]

tSpts=Spts.{{0,-1},{1,0}}(*もとの画像データは(-y,x)点の点群データであるので，90度回転してxy座標をあわせる*);
(*ランダムに点を選ぶ場合*)
(*Nsample(*点群の要素数*)=500;
tSpts=RandomReal[{0,1},{Nsample,2}];*)
DT=DelaunayMesh[tSpts]
Print["Numbers of cells =",MeshCellCount[DT]]
vertices=MeshCoordinates[DT];
(*-------境界作用素の計算-------*)
lines(*lineの要素をソート*)=Sort[Sort/@MeshCells[DT,1][[All,1]]];
faces(*faceの要素をソート*)=Sort/@MeshCells[DT,2][[All,1]];
d1(*境界作用素:{line}->{vertex}*)=SparseArray[Reap[Do[Sow[{j,i}->1],{i,Length[lines]},{j,lines[[i]]}]][[2,1]]];
facebdy(*faceを境界に分割し,全てを順番に並べる*)=Reap[Do[
Sow[{f[[1]],f[[2]]}];Sow[{f[[1]],f[[3]]}];Sow[{f[[2]],f[[3]]}],{f,faces}];
][[2,1]];
permd2(*faceの境界のlineをソートする置換*)=Ordering[facebdy];
lis(*ソートされたlineの出現数∈{1,2}を数える*)=Tally[Sort[facebdy]][[All,2]];
d2(*境界作用素:{face}->{line}*)=SparseArray[Reap[n=0;
Do[n++;Sow[{i,Ceiling[permd2[[n]]/3]}->1];,{i,Length[lis]},{j,lis[[i]]}];
][[2,1]]];
bdy(*全体の境界作用素*)=ArrayFlatten[{{SparseArray[{},{Length[vertices],Length[vertices]}],d1,0},{0,0,d2},{0,0,SparseArray[{},{Length[faces],Length[faces]}]}}];
(*-------セルの出現順の計算-------*)
(*0セルは時刻0で出現．1セルは長さ/2で出現．2セル（三角形）は外接円の半径で出現*)
(*線の集合の長さ/2の配列*)
dist=(Norm[Differences[vertices[[#]]]]&/@lines)/2;
(*三角形の外接円の半径の配列. Circumsphereを使った方が早い*)
radius=(Circumsphere[vertices[[#]]]//Last)&/@faces;
(*{{線の長さ/2},{外接円の半径}}をソートする置換行列*)
perm=Ordering[Join[dist,radius]];
P(*置換行列．元の基底 B に対して，B'=BP が出現順の基底*)=ArrayFlatten[{{SparseArray[{{i_,i_}->1},{Length[vertices],Length[vertices]}],0},{0,SparseArray[Table[{perm[[i]],i}->1,{i,Length[perm]}]]}}];
bdy(*出現順で基底変換された境界作用素*)=P\[Transpose].bdy.P;
(*ソート後の各セルの次元*)
vdim=Join[Table[0,Length[vertices]],Table[1,Length[lines]],Table[2,Length[faces]]].P;

(*出現順にセルを表示*)
cells=Join[Point/@vertices,Line/@(vertices[[#]]&/@lines),Polygon/@(vertices[[#]]&/@faces)];
cells=cells.P;
Manipulate[Graphics[cells[[1;;time]]],{time,Length[vertices],Length[cells],1}]

(*-------パーシステント図の計算-------*)
L[sa_]:=Module[{rule},rule=ArrayRules[sa];If[Length[rule]==1,Return[0],Return[First[First[rule[[-2]]]]]]];
m(*Mathematicaで扱いやすいように転置を取る*)=bdy\[Transpose];
Nvertices=Length[vertices];
Do[
l=L[m[[j]]];
If[l>0,
i=Nvertices+1(*頂点はスキップ*);
While[i<j,
If[L[m[[i]]]==l,
m[[j]]=SparseArray[DeleteCases[ArrayRules[Mod[m[[j]]+m[[i]],2]],_->0],Dimensions[m[[j]]]];
l=L[m[[j]]](*lを再定義*);
If[l==0,Break[];];
i=Nvertices(*iをリセット*);
];
i++;
];
];
If[Mod[j,100]==0,Print[j]];
,{j,Nvertices+1,Length[m]}]
Print["Done"];
(*有限時間で生まれ死ぬセル*)
finite=Reap[Do[If[L[m[[j]]]>0,Sow[{L[m[[j]]],j}]],{j,Length[m]}]][[2,1]];
(*∞まで生き残るセル*)
inf=Complement[Position[L/@m,0]//Flatten,Select[L/@m,#>0&]//Sort];
inf=Table[{i,∞},{i,inf}];
PHD(*パーシステント図*)=Join[finite,inf];

tPHD(*セルの次元に分ける.2セルは消えないので不要*)=Table[Select[PHD,vdim[[#[[1]]]]==d&],{d,{0,1}}];
times=Sort[Join[Table[0,Length[vertices]],dist,radius]];
tPHD=Map[If[#≠∞,times[[#]],∞]&,tPHD,{3}];

(*パーシステント図*)
Show[
ListPlot[tPHD[[1]],PlotStyle->{Red,PointSize[0.02]},PlotMarkers->"o",PlotLegends->{"vertex"}],
ListPlot[tPHD[[2]],PlotStyle->{Black,PointSize[0.02]},PlotMarkers->"x",PlotLegends->{"loop"}],PlotRange->All,AxesLabel->{"birth","death"}]