Numerické hledání Steinerova bodu pro konvoluci s gaussovským jádrem

V tomto souboru se nachází definice potřebných funkcí (definována je množina A, konvoluční jádro, funkce příslušnosti fuzzy množiny vzniklé konvolucí s daným jádrem apod.). 

Kromě definic jsou tu funkce pro výpočet bodů na alpha hladině i numerickou integraci pro výpočet Steinerova bodu. 

Na konci je rozsáhlá vizualizace. 

(Před voláním funkce doItAll je nutné inicializovat všechny funkce)
(soubory jsou ukladany / nacitany do / z domovskeho adresare - resp. do / z defaultne nastavene cesty v programu Mathematica)

Nastavení

Množina A

(* geometricky objekt, pres ktery se bude integrovat -- volano jako Element[{x,y},setA] *)

setA = Disk[{0,0},1,{0,Pi}];

Konvoluce s vybraným jádrem

Nastavení parametrů jádra

PARDSC ="sig1_2"; (* popis jadra -- pro ucely pojmenovani souboru *)
muX = 0;muY = 0; sigma = {{1/2,0},{0,1/2}};
(* kernel[x,y] je funkce jadra pouzita pri konvoluci *)
kernel[x_?NumericQ,y_?NumericQ] = PDF[MultinormalDistribution[{muX,muY},sigma],{x,y}];

Konvoluce a její maximum

(* konvoluce z definice jako integral*)
convolution[x_ ?NumericQ,y_?NumericQ]:=NIntegrate[kernel[x-u,y-v], Element[{u,v},setA]];
(* maximum -- prevracena hodnota pouzita jako normalizacni clen *)
{muStartX,muStartY} = NArgMax[convolution[x,y],{x,y}];
maximum = convolution[muStartX,muStartY];

Funkce příslušnosti fuzzy množiny vzniklé konvolucí

fuzzyMemFunction[x_?NumericQ,y_?NumericQ] := (1/maximum) * convolution[x,y];
fuzzyMemFunction::usage = "fuzzyMemFunction[x_?NumericQ,y_?NumericQ] vrati stupen prislusnosti daneho bodu (x,y) k nasi fuzzy mnozine";

Další parametry použité v průběhu výpočtu (globální proměnné):

kernelName = "Gauss"; (* nazev jadra - pro pojmenovano souboru *)
minPoints = 32; (* minimalni pocet bodu hledanych na α-hladine *)
maxPoints = 128; (* maximalni pocet bodu hledanych na α-hladine *)
PHI = 0; (* pootoceni pocatecniho smeru pri generovani smerovych vektoru, volano v generateCirclePoints *)
PRECISION = 100; (* presnost pri volani N[..., PRECISION] *)
NUMBERofPOINTS = maxPoints; (* pocet bodu hledanych na α-hladine *)
EPSILON = (10^(-10)); (* maximalni chyba pri hledani  bodu na hladinach *)
maxITER = 500; (* maximalni pocet iteraci pri hledani 1 bodu na α-hladine *)
integrationMETHOD = NrombergY;(* vyber integracni metody *)
(* NrombergX[NoScale] NrombergY[NoScale] NsimpsonY NtrapezoidX[NoScale] NtrapezoidY[NoScale] *)
integrationMETHODNS = NrombergYNoScale;(* poze NoScale metody*)

Nastavení jmen souborů a hladin, které se budou hledat

(* hledane vysky α *)
HEIGHTS = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9};

(* popisky souboru pro body na hladinach *)

heightsLABELS = Map[StringJoin["LAYER-",kernelName,"-",PARDSC,"-",ToString[NUMBERofPOINTS],
"-",ToString[#,InputForm],".dat"]&,HEIGHTS] // StringReplace[#,"/":>"_"]&;
(* popisky souboru pro Steinera *)
steinerLABELS = 
Map[StringJoin["STEINER-",kernelName,"-",PARDSC,"-",ToString[NUMBERofPOINTS],
"-",ToString[#,InputForm],".dat"]&,HEIGHTS] // StringReplace[#,"/":>"_"]&;

Základní vizualizace fuzzy množiny

Pokud byly funkce v sekci nastavení inicializovány správně, následuje obrázek fuzzy množiny a jejích vybraných hladin:

Fuzzy množina

Plot3D[fuzzyMemFunction[x,y],{x,-2,2},{y,-2,2},PlotRange->All]

Vybrané hladdiny

ContourPlot[fuzzyMemFunction[x,y],{x,-2,2},{y,-2,2},Contours->{0.2,0.4,0.6,0.8},AspectRatio->Automatic,PlotRange->All]

Výpočet

Pro spuštění výpočtu stačí po inicializaci všech funkcí zavolat následující funkci:

doItAll[HEIGHTS,heightsLABELS,steinerLABELS,NUMBERofPOINTS, EPSILON,maxITER]

printVisuals[HEIGHTS,heightsLABELS,steinerLABELS]

Funkce nutné pro výpočet

Generování směrových vektorů

generateCirclePoints[number_] := Table[{Cos[x + PHI],Sin[x + PHI]},{x,0,2Pi,2*Pi/number}][[;;-2]];
generateCirclePoints::usage="generateCirclePoints[number_] trati tabulku [number] rovnomerne rozdelenych souradnic smerovych vektoru";

Hledání bodů na α-hladině

Hledání jednoho bodu v daném směru

(dirX,dirY):    představuje směr hledání,
height:          odpovídá výšce α  hledané α-hladiny,
epsilon > 0:  je maximální rozdíl mezi stupněm příslustnosti nalezeného bodu a height,
maxIter:        je maximální počet itrací

findLayerPoint[dirX_,dirY_,height_,epsilon_,maxIter_] :=
Module[{curX,curY,nextX,nextY, currentValue, nextValue, alpha,iterations,shiftX,shiftY},
If[height ≤ 0 || height > 1, Print["Hladina musí být mezi (0,1]"];Return[{}]];

{curX,curY} = {muStartX,muStartY}; (* pocatecni bod hledani (Subscript[x, 0]) *)
currentValue = fuzzyMemFunction[curX,curY]; (* pocatecni hodnota -- bude rovna 1 *)
iterations = 0; (* pocet iteraci *)
alpha = 1.; (* skalovaci faktor -- delka kroku *)


(* hledej dokud neprekrocis maximalni pocet iteraci nebo neni dosazeno pozadovane presnosti *)
While[Abs[currentValue-height]≥ epsilon && iterations < maxIter,
(* vypocti dalsi potencialni bod (Subscript[x, i+1]) *)
{nextX,nextY} = {curX,curY} + alpha * {dirX,dirY};

(* vypocti stupen prislusntosti v (Subscript[x, i+1]) *)
nextValue = fuzzyMemFunction[nextX,nextY];

(* pokud se priblizil cili a stupen prislusnosti je mensi -- presun se do nej a zvys alpha*)
If[nextValue ≥ height, 
(* kdyz ano *)
(
currentValue = nextValue;
{curX,curY} = {nextX,nextY};
alpha *= 2;
),
(* kdyz ne, sniz alpha *)
alpha /= 2;
];

iterations += 1;
];

(*Print["Počet iterací: ",iterations];*)
If[iterations==maxIter,Print["Pozor! Dosazeno maximalniho poctu iteraci ve smeru: ",{dirX,dirY}]];
Return[{curX,curY}];
];

findLayerPoint::usage = "findLayerPoint[dirX_,dirY_,height_,epsilon_,maxIter_]  najde bod, ve kterém je stupen prislusnosti fuzzy mnoziny (hodnota fuzzyMemFunction) priblizne roven pozadovane height s maximalni odchylkou epsilon, maxIter je maximalni pocet iteraci hledani, {dirX,dirY} smer hledani";

Hledání všech bodů na dané hladině

findAllLayerPoints[number_,height_,epsilon_,maxIter_] :=
Module[{directions,layerPoints,i,shiftX,shiftY},
If[height>1 || height ≤0,Print["Hladina musí být mezi (0,1]"];Return[{}]];

layerPoints = Parallelize[Map[findLayerPoint[#[[1]],#[[2]],height,epsilon,maxIter] &,N[generateCirclePoints[number], PRECISION] ]];
Return[layerPoints];
];
findAllLayerPoints::usage = "findAllLayerPoints[number_,height_,epsilon_,maxIter_] vrati seznam bodu fuzzy mnoziny na zadane hladine";

Hodnota opěrné funkce násobená směrem x/y

findSupport[{ex_,ey_},pts_] := Max[Map[({ex,ey} . #) &, pts]]
findSupport::usage="findSupport[{ex,ey},pts] najde maximum skalarniho soucinu smeru (ex,ey) a bodu z pts";

Hodnoty integrandu při výpočtu Steinerova bodu (hodnota opěrné funkce násobená souřadnicí směru):

supportTimesX[{ex_,ey_},pts_] := ex * findSupport[{ex,ey},pts];
supportTimesY[{ex_,ey_},pts_] := ey * findSupport[{ex,ey},pts];

Numerická integrace

Lichoběžníková metoda

Pro y-ovou souřadnici:

NtrapezoidY[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesY[#,N[pts[[;;;;Length[pts]/num ]],PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi; 
(* / (Pi * (1^2)) -- tj. obsah jednotkoveho kruhu -- konstanta pri vypoctu Steinerova bodu 2D *)

NtrapezoidY::usage="NtrapezoidY[num,pts] -- vygeneruje [num] uzlovych bodu na jednotkove kruznici a pocita y-ovou souradnici Steinerova bodu pomoci lichobeznikove metody; hodnoty operne funkce hledany pres pts (nikoliv celou A), delka pts musi byt mocninou dvou a minimalne maxPoints";

NtrapezoidYNoScale[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesY[#,N[pts,PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi;

Pro x-ovou souřadnici:

NtrapezoidX[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesX[#,N[pts[[;;;;Length[pts]/num ]],PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi; 
NtrapezoidX::usage="NtrapezoidY[num,pts] -- vygeneruje [num] uzlovych bodu na jednotkove kruznici a pocita x-ovou souradnici Steinerova bodu pomoci lichobeznikove metody; hodnoty operne funkce hledany pres pts (nikoliv celou A), delka pts musi byt mocninou dvou a minimalne maxPoints";

NtrapezoidXNoScale[num_,pts_] := ((2*Pi)/num)*Total[ParallelMap[
Quiet[supportTimesX[#,N[pts,PRECISION]]] &,
N[generateCirclePoints[num],PRECISION]]] / Pi;

Simpsonova metoda

NsimpsonY[num_,pts_] := 
(((2*Pi)/num)/3)* (1/ (Pi))*
Total[Table[If[Mod[i,2]==0,2,4],{i,0,num-1}] * 
ParallelMap[
(supportTimesY[ #,N[pts,PRECISION]]) &,N[generateCirclePoints[num],PRECISION]]];

Rombergova metoda

Při volání této metody je třeba, aby délka [pts] byla mocninou 2 a byla >= maxPoints, jinak je třeba volat NrombergYNoScale

NrombergY[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidY[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];
NrombergYNoScale[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidYNoScale[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

Stejně pro x-ovou souřadnici

NrombergX[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidX[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

NrombergXNoScale[maxPoints_,pts_] :=
Module[{height,table, layerPoints,i,j},
height = Log[2,maxPoints] - Log[2,minPoints] +1;
table = Array[0 &,{height,height}];

(* vypln prvni sloupec *)
For[i=1,i≤ height,i++,table[[i,1]] = NtrapezoidXNoScale[(2^((Log[2,minPoints]+i-1))),pts] ];

(* vypln zbytek tabulky *)
For[j=2,j≤ height, j++,
For[i=j,i≤ height, i++,
table[[i,j]] = table[[i,j-1]] + (((table[[i,j-1]] - table[[i-1,j-1]]))/((4^((j-1)))-1))
]];
Return[table];
];

Odhad chyby přes “horní” a “spodní” mnohoúhelník

Opsaný mnohoúhelník L

<<ComputationalGeometry`

Funkce, která vrátí další/předchozí index (v kruhové smyčce)

nextPlus[number_,modulo_] := Mod[number,modulo]+1;
nextMinus[number_,modulo_] := If[number-1≤ 0,modulo,number-1];

Průsečík přímek daný body s indexy (i,i+1) a (i+2,i+3)

getIntersect[points_,i_]:=
Module[{aa,bb,cc,dd,AA,BB,CC,DD,mod},
mod=Length[points];
{aa,bb} = points[[i]];
{cc,dd} = points[[nextPlus[i,mod]]];
{AA,BB} = points[[ nextPlus[nextPlus[i,mod],mod] ]];
{CC,DD} = points[[ nextPlus[nextPlus[nextPlus[i,mod],mod],mod] ]];
 Return[{aa,bb}+(((-bb+BB) CC+AA (bb-DD)+aa (-BB+DD))/(-BB cc+bb CC-CC dd+AA (-bb+dd)+aa (BB-DD)+cc DD))*{aa-cc,bb-dd}];
];
getIntersect::usage="getIntersect[points,i] vrati prusecik primek dany body s indexy (i,i+1) a (i+2,i+3)";

getHull[allPoints_] :=
Module[{intersect,i,mod,hull,points},
points = allPoints[[ConvexHull[allPoints]]];
mod = Length[points];
hull = {};

For[i=1,i≤ mod,i++,
intersect = getIntersect[points,i];
AppendTo[hull,intersect]
];

Return[hull];
];
getHull::usage="getHull[points] vrati vrcholy mnohouhelnika opsaneho mnohouhelniku urceneho vrcholy [points]";

Horní a spodní část

getUpperAndLower[points_] := Module[{hullOfPoints, mapping,iMaxRight,iMaxLeft, rightUpper, leftUpper,rightLower,leftLower, iRightUpper,iLeftUpper,iRightLower,iLeftLower, upperHull,lowerHull,upperPoint,lowerPoint,xCoordRight,xCoordLeft},
hullOfPoints = getHull[points];
mapping = Map[({1,0}. # ) &, points];
(* indexy bodu s maximalnim kolmym prumetem do (1,0) a (-1,0) z [points] *)
iMaxRight = Ordering[mapping,-1][[1]];
iMaxLeft = Ordering[mapping,1][[1]];

(* body a indexy pro opsany mnohouhelnik*)
xCoordRight = points[[iMaxRight]][[1]];
xCoordLeft = points[[iMaxLeft]][[1]];

(* pro problematicky pripad konvoluce s kruhem *)
(*
(* iteruj pravy horni bod *)
While[True,
rightUpper =  getIntersect[points,nextMinus[iMaxRight,Length[points]]];
rightLower = getIntersect[points,nextMinus[nextMinus[iMaxRight,Length[points]],Length[points]]];
If[rightUpper[[1]] ≤ xCoordRight, Break[]];
iMaxRight = nextPlus[iMaxRight,Length[points]];
];

(* iteruj levy hodni bod*)
While[True,
leftUpper = getIntersect[points,nextMinus[nextMinus[iMaxLeft,Length[points]],Length[points]]];
leftLower =    getIntersect[points,nextMinus[iMaxLeft,Length[points]]];
If[leftUpper[[1]] ≥  xCoordLeft, Break[]];
iMaxLeft = nextMinus[iMaxLeft,Length[points]];
];
*)

rightUpper =  getIntersect[points,nextMinus[iMaxRight,Length[points]]];
leftUpper =    getIntersect[points,nextMinus[nextMinus[iMaxLeft,Length[points]],Length[points]]];
rightLower = getIntersect[points,nextMinus[nextMinus[iMaxRight,Length[points]],Length[points]]];
leftLower =    getIntersect[points,nextMinus[iMaxLeft,Length[points]]];

iRightUpper = Position[hullOfPoints,rightUpper][[1,1]];
iLeftUpper =   Position[hullOfPoints,   leftUpper][[1,1]];
iRightLower = Position[hullOfPoints,rightLower][[1,1]];
iLeftLower =   Position[hullOfPoints,   leftLower][[1,1]];

upperHull = If[iRightUpper > iLeftUpper, 
Join[hullOfPoints[[iRightUpper;;Length[hullOfPoints]]], hullOfPoints[[1;;iLeftUpper]]],
Join[hullOfPoints[[iRightUpper;;iLeftUpper]]]
];
	lowerHull = If[iRightUpper > iLeftUpper, 
hullOfPoints[[iLeftUpper+1;;iRightUpper]][[1;;-2]],
Join[hullOfPoints[[iLeftUpper+1;;Length[hullOfPoints]]],hullOfPoints[[1;;iRightUpper]][[1;;-2]]]
];
	upperPoint =  If[iMaxRight > iMaxLeft, 
Join[points[[iMaxRight;;Length[points]]], points[[1;;iMaxLeft]]][[2;;-2]],
Join[points[[iMaxRight;;iMaxLeft]]][[2;;-2]]
];
	lowerPoint =  If[iMaxRight > iMaxLeft, 
points[[iMaxLeft+1;;iMaxRight]][[1;;-2]],
Join[points[[iMaxLeft+1;;Length[points]]],points[[1;;iMaxRight]][[1;;-2]]]
];

Return[{Join[upperHull,lowerPoint],Join[upperPoint,lowerHull]}] (* vrat sady bodu pro {horni,spodni} odhad *)
];

Vizualizace množin pro horní a spodní odhad

visualiseHull[points_] := Module[{hullOfPoints,iMaxRight,iMaxLeft,upperSet,lowerSet,mapping},
hullOfPoints = getHull[points];
{upperSet,lowerSet} = getUpperAndLower[points];

mapping = Map[({1,0}. # ) &, points];
iMaxRight = Ordering[mapping,-1][[1]];
iMaxLeft = Ordering[mapping,1][[1]];

Print[
Graphics[{
Black,
Point[points],
Point[hullOfPoints],
Line[Append[points,points[[1]]]],
Line[Append[hullOfPoints,hullOfPoints[[1]]]] 
}]
];

Print[
Graphics[{
PointSize[Small],
Black,
Point[points],
Point[hullOfPoints],
Line[Append[points,points[[1]]]],
(*Line[Append[hullOfPoints,hullOfPoints[[1]]]] ,*)
Text[Subscript[k,"+"],points[[iMaxRight]] - {0.06,0}],
Text[Subscript[k,"-"],points[[iMaxLeft]]+{0.06,0}],
Red, 
Line[Append[upperSet,upperSet[[1]]]],
Blue,
Line[Append[lowerSet,lowerSet[[1]]]]
}]
];
]

Hledání hladin, Steinerových bodů, vizualizace ...

Hledání hladin

findAndSaveLayers[heights_,heightsLabels_,numberOfPoints_, epsilon_,maxIter_]:=
Module[{index,data},
For[index=1,index≤Length[heights],index++,
Print["Current index: ",index, ", height ", heights[[index]], ", filename: ",heightsLabels[[index]]];
data = findAllLayerPoints[numberOfPoints,heights[[index]],epsilon,maxIter];
WriteString["stdout","\n"];
Export[heightsLabels[[index]],data];
];
]
findAndSaveLayers::usage = "findAndSaveLayers[heights,heightsLabels,numberOfPoints, epsilon,maxIter] najde [numberOfPoints] bodu na kazde z [heights] hladin s presnosti [epsilon], ulozi je do soboru [heightsLabels], pri hledani kazdeho bodu nesmi prekrocit [maxIter], jinak vypise varovani";

Počítání Steinera

computeAndSaveSteiner[heightsLabels_,steinerLabels_] :=
Module[{layerData,layerSteiner,romTable,steinerUpperEst,steinerLowerEst,estTable,upperSet,lowerSet,i,j,tALower,tAUpper},
layerData = Array[0 &,Length[heightsLabels]];
layerSteiner = {};steinerUpperEst = {}; steinerLowerEst= {};
Print["Load data and compute Steiner -- start"];

Print["  -- Importing data"];
(* nacti data *)
For[i=1,i≤Length[heightsLabels],i++, layerData[[i]] = Import[heightsLabels[[i]]]];

Print["  -- Computing Steiner points"];
(* spocitej steinerovy body *)
For[i=1,i≤Length[layerData],i++, 
romTable = integrationMETHOD[maxPoints,layerData[[i]] ];
AppendTo[layerSteiner,Diagonal[romTable] ];

(* spocitej odhady  *)
tALower = {}; tAUpper={};
For[j=minPoints,j≤maxPoints,j*=2,
{upperSet,lowerSet} = getUpperAndLower[layerData[[i]][[;;;; Length[layerData[[i]]] /j ]] ];
estTable = integrationMETHODNS[maxPoints,upperSet];
AppendTo[tAUpper,estTable[[Length[estTable],Length[estTable]]] ];

estTable = integrationMETHODNS[maxPoints,lowerSet];
AppendTo[tALower,estTable[[Length[estTable],Length[estTable]]] ];
];

AppendTo[steinerUpperEst,tAUpper];
AppendTo[steinerLowerEst,tALower];
];

Print["  -- Exporting results"];
(* uloz Steinerovy body *)
For[i=1,i≤Length[layerSteiner],i++, Export[steinerLabels[[i]],{layerSteiner[[i]] , steinerUpperEst[[i]],steinerLowerEst[[i]]}]
];
];
computeAndSaveSteiner::usage = "computeAndSaveSteiner[heightsLabels,steinerLabels] spocita steinerovy body z bodu na hladinach, ktere nacte ze souboru [heightsLabels] a ulozi je do souboru [steinerLabels]";

Vizualizace

printVisuals[heights_,heightsLabels_,steinerLabels_] :=
Module[{layerData,layerSteiner,i,layerData3D,plot3Dlist,plot3Dorig, layerSteinerAll,contours,xMin,xMax,yMin,yMax, layersGraphis,layerSteinerUpper,layerSteinerLower,importBuf,upperEsts,lowerEsts,estGraphic},
layerData = Array[0 &,Length[heightsLabels]];
layerSteiner = {};layerSteinerUpper={}; layerSteinerLower={};upperEsts={}; lowerEsts={};estGraphic={};

Print["  -- Importing data"];
(* load layer data *)
For[i=1,i≤Length[layerData],i++, layerData[[i]] = Import[heightsLabels[[i]]]];
For[i=1,i≤Length[steinerLabels],i++,
importBuf = Import[steinerLabels[[i]] ];
AppendTo[layerSteiner,importBuf[[1]]];
AppendTo[layerSteinerUpper,importBuf[[2]]];
AppendTo[layerSteinerLower,importBuf[[3]]];

AppendTo[upperEsts,Flatten/@Transpose[{Array[#&,Length[layerSteinerUpper[[i]] ]],layerSteinerUpper[[i]]}] ];
AppendTo[lowerEsts,Flatten/@Transpose[{Array[#&,Length[layerSteinerLower[[i]] ]],layerSteinerLower[[i]]}] ];

];

	(* plot range *)
xMin = MinimalBy[layerData[[1]],First][[1,1]];
       xMax = MaximalBy[layerData[[1]],First][[1,1]];
       yMin = MinimalBy[layerData[[1]],Last][[1,2]];
       yMax = MaximalBy[layerData[[1]],Last][[1,2]];

Print["  -- Visualisation"];

Print["     :plot all layers"];
Print[ListPlot[layerData,AspectRatio->Automatic,PlotStyle->Map[ColorData["Rainbow"][#]&,Prepend[Array[1/(Length[layerData]-1)*#&,Length[layerData]-1],0]]]];

Print["     :plot layers in 3D"];
layerData3D = Array[0 &,Length[layerData]];
	For[i=1,i≤Length[layerData],i++,
layerData3D[[i]] = MapThread[Append,{layerData[[i]],Array[heights[[i]] &, Length[layerData[[i]] ] ]}]
];
plot3Dlist =ListPointPlot3D[Flatten[layerData3D,1],AxesLabel->{"x",Framed["y",FrameStyle->None,FrameMargins->10],Framed["z",FrameStyle->None,FrameMargins->10]}, ColorFunction->"Rainbow",PlotRange->All,AspectRatio->Automatic];
Print[plot3Dlist];

(* pro ucely exportu pro 10 hladin -- namisto bodu vykresluje cary -- setri pamet po exportu*)
(*
plot3Dlist =
Graphics3D[{
Thickness[0.01],
ColorData["Rainbow"][0],
Line[layerData3D[[1]]],
ColorData["Rainbow"][1/8],
Line[layerData3D[[2]]],
ColorData["Rainbow"][2/8],
Line[layerData3D[[3]]],
ColorData["Rainbow"][3/8],
Line[layerData3D[[4]]],
ColorData["Rainbow"][4/8],
Line[layerData3D[[5]]],
ColorData["Rainbow"][5/8],
Line[layerData3D[[6]]],
ColorData["Rainbow"][6/8],
Line[layerData3D[[7]]],
ColorData["Rainbow"][7/8],
Line[layerData3D[[8]]],
ColorData["Rainbow"][1],
Line[layerData3D[[9]]]
},PlotRange->All,Axes->True,AxesLabel->{"x",Framed["y",FrameStyle->None,FrameMargins->10],Framed["z",FrameStyle->None,FrameMargins->10]},BoxRatios->{1, 1, 1/2}];
Print[plot3Dlist];
*)

plot3Dorig =Plot3D[fuzzyMemFunction[x,y],{x,xMin,xMax},{y,yMin,yMax},PlotRange->All,AxesLabel->{"x","y",None},AspectRatio->Automatic];
Print[Show[{plot3Dlist,plot3Dorig}]];

Print[plot3Dorig];

Print["     :contours"];
contours = ContourPlot[fuzzyMemFunction[x,y],{x,xMin,xMax},{y,yMin,yMax},Contours->heights,AspectRatio->Automatic,PlotLegends->Automatic];
Print[contours];
Print[Show[{contours, ListLinePlot[layerData,AspectRatio->Automatic]}]];


Print["     :visualize Steiner points"];
layerSteinerAll = {};
For[i=1,i≤ Length[layerSteiner],i++,AppendTo[layerSteinerAll, layerSteiner[[i,-1]]];];

Print[layerSteinerAll];

Print[ListPlot[layerSteinerAll,Frame->True,GridLines->Automatic,FrameLabel->{{"(Subscript[s, y])",None},{None,"α"}}, FrameTicks->{{Automatic,None},{None,Transpose[{Array[# &, Length[heights]],Map[N[#] &,heights]}] }} ]];

Print["     :table of Steiners"];
Print[TableForm[Transpose[layerSteiner],TableHeadings->{
Array[StringJoin["Bodů: ",ToString[2^(#+Log[2,minPoints]-1)]]&,Log[2,maxPoints]-Log[2,minPoints]+1],Map[StringJoin["α = ", ToString[N[#]]] &,heights]
}]];

Print["     :error estimates"];

For[i=1,i≤Length[layerSteiner],i++,
AppendTo[estGraphic,
Show[{
ListPlot[layerSteiner[[i]],PlotStyle->{Black,PointSize[Medium]},Ticks->{Flatten/@Transpose[{Array[#&,Length[Array[ToString[2^(#+Log[2,minPoints]-1)]&,Log[2,maxPoints]-Log[2,minPoints]+1]]],Array[ToString[2^(#+Log[2,minPoints]-1)]&,Log[2,maxPoints]-Log[2,minPoints]+1]}] ,Automatic}],
Graphics[{
PointSize[Medium],
Red,
Point[upperEsts[[i]]],
Blue,
Point[lowerEsts[[i]]]
}],
ListPlot[layerSteiner[[i]],PlotStyle->{Black,PointSize[Medium]}]
},PlotRange->All,PlotLabel->ToString[heightsLabels[[i]]],AxesLabel->{N,(Subscript[s, y])}]
];
];

Print[estGraphic];

Print["     :additional graphics"];

layersGraphis = {};
For[i=1,i≤  Length[layerData],i++,
AppendTo[layersGraphis,Graphics[{
Thickness[0.01],
Line[Append[layerData[[i]],layerData[[i,1]]]],
Gray,
Point[{{0,0}}],
PointSize[Large],
Red,
Point[{0,layerSteinerAll[[i]]}],
Black
(*Inset[ StringJoin["Hladina ", ToString[N[layerHeights[[i]]]] ] ,{0,0.05}]*)
},AspectRatio->Automatic,Axes->True, PlotLabel->StringJoin["Hladina ", ToString[N[heights[[i]]]], ", sy = ", ToString[layerSteinerAll[[i]]]],PlotRange->{{xMin-0.1,xMax+0.1},{yMin-0.1,yMax+0.1}}]
];
];

Print[layersGraphis]
];

printVisuals::usage="printVisuals[heights,heightsLabels,steinerLabels] provede kompletni vizualizaci - body na hladinach a vypoctene steinerovy body nacte se souboru [heightsLabels] resp. [steinerLabels]";

Volání funkce, která najde body, spočítá steinerovy body, exportuje a vizualizuje

doItAll[heights_,heightsLabels_,steinerLabels_,numberOfPoints_, epsilon_,maxIter_] :=
Module[{},
findAndSaveLayers[heights,heightsLabels,numberOfPoints, epsilon,maxIter];
computeAndSaveSteiner[heightsLabels,steinerLabels];
printVisuals[heights,heightsLabels,steinerLabels];
];
doItAll::usage="doItAll[heights,heightsLabels,steinerLabels,numberOfPoints, epsilon,maxIter] vola funkce findAndSaveLayers, computeAndSaveSteiner, printVisuals s prislusnymi parametry";