(* In the following code, we calculate the box counting dimension of an image. We do this by loading the files, set a treshold for the maximum value that should be taken into account, binarize the files and partition the image on a rectangular grid. Next, we count how many partitions of the image contain a white element, such that we obtain a number of boxes on a square grid that is part of the fractal set. We plot the logs of the size and the amount of boxes and calculate the slope. The counting of the boxes is based on https://mathematica.stackexchange.com/questions/13125/measuring-fractal-dimension-of-natural-objects-from-digital-images
Sander Kempkes 02-11-2020 *)

(*Load the files*)

SetDirectory[NotebookDirectory[]];

filenames=FileNames["Data\\Sierpinski3_2-20170413-102354--259_*-bwd.txt"];

datatot=Table[Import[filenames[[i]],"Table"],{i,1,Length[filenames]}];

ListDensityPlot[datatot[[1]],PlotRange->All]

e=Flatten[{Table[-0.2+0.025 i,{i,0,20}],Table[-0.325+0.025i,{i,0,4}]}];

(*energy of the maps*)

e

(*prep the images*)

dataprep=ParallelTable[Table[{j,-i,If[i<157&&i>(j)*Sqrt[3]-134&&i>-(j)*Sqrt[3]+209,If[i>101&&i<(j)*Sqrt[3]-23&&i<-(j)*Sqrt[3]+317||i>75&&i<(j)*Sqrt[3]-79&&i<-(j)*Sqrt[3]+262||i>129&&i<(j)*Sqrt[3]+31&&i<-(j)*Sqrt[3]+262||i>129&&i<(j)*Sqrt[3]-79&&i<-(j)*Sqrt[3]+370,0,datatot[[l,i,j]]-Min[datatot[[l]]]],0]},{i,1,Length[datatot[[l]]]},{j,1,200}],{l,1,Length[datatot]}];

datanomaskprepnomask=ParallelTable[Table[{j,-i,datatot[[l,i,j]]-Min[datatot[[l]]]},{i,1,Length[datatot[[l]]]},{j,1,200}],{l,1,Length[datatot]}];

ListDensityPlot[Flatten[dataprep[[19]],1],PlotRange->All]

ListDensityPlot[Flatten[datanomaskprepnomask[[20]],1],PlotRange->All]

max=Table[Max[Flatten[dataprep[[i]],1][[All,3]]],{i,1,Length[datatot]}];

max

min=Table[Min[Flatten[dataprep[[i]],1][[All,3]]],{i,1,Length[datatot]}];

min

(*select a certain threshold of the data files*)

dat=Table[Table[{dataprep[[i,j,l,1]],dataprep[[i,j,l,2]],If[dataprep[[i,j,l,3]]>0.65*max[[i]],1,0]},{j,1,200},{l,1,200}],{i,1,Length[dataprep]}];

img=Table[ListDensityPlot[Flatten[dat[[i]],1],Frame->False,AspectRatio->Automatic],{i,1,Length[dataprep]}];

(*resulting images*)

(*Compose the background*)

image3=Table[ImageCompose[ImageResize[ImageCrop[Graphics[Rectangle[{0,0},{1,1}]]],{360,360}],ImageResize[img[[i]],1.45*ImageDimensions[img[[i]]]],Scaled[{0.515,0.4}]],{i,1,Length[datatot]}];

image3

(*Binarize the image*)

{Table[Binarize[image3[[i]],0.75],{i,1,Length[image3]}],Binary=Table[Binarize[image3[[i]],0.5],{i,1,Length[image3]}]};

Binary

(*Count the number of white elements for the given size of the radius of the circle. For each regime there are different step sizes. Small sizes take a longer time, since the partitioning into smaller pieces creates more elements in the list*)

boxjes=Reverse[{1,2,3,4,5,6,8,9,10,12,15,18,20,24,30,36,40,45,60,72,90}];

datafit=Table[{ParallelTable[{1/boxjes[[size]],Total[Sign/@(Total[#,2]&/@(ImageData/@Flatten[ImagePartition[Binary[[i]],boxjes[[size]]]]))]},{size,1,21,1}]},{i,1,Length[Binary]}];

(*Linearfit to the log of the data, slope gives the fractal dimension*)

linfit=Table[LinearModelFit[Log[Flatten[datafit[[i]],1][[12;;21]]],x,x],{i,1,Length[datafit]}];
linfitb=Table[LinearModelFit[Log[Flatten[datafit[[i]],1][[7;;21]]],x,x],{i,1,Length[datafit]}];
linfitc=Table[LinearModelFit[Log[Flatten[datafit[[i]],1][[1;;21]]],x,x],{i,1,Length[datafit]}];
linfitd=Table[LinearModelFit[Log[Flatten[datafit[[i]],1][[3;;17]]],x,x],{i,1,Length[datafit]}];
linfite=Table[LinearModelFit[Log[Flatten[datafit[[i]],1][[1;;11]]],x,x],{i,1,Length[datafit]}];
linfitf=Table[LinearModelFit[Log[Flatten[datafit[[i]],1][[4;;15]]],x,x],{i,1,Length[datafit]}];
dim=Table[linfit[[i]][1]-linfit[[i]][0],{i,1,Length[datafit]}];
dimb=Table[linfitb[[i]][1]-linfitb[[i]][0],{i,1,Length[datafit]}];
dimc=Table[linfitc[[i]][1]-linfitc[[i]][0],{i,1,Length[datafit]}];
dimd=Table[linfitd[[i]][1]-linfitd[[i]][0],{i,1,Length[datafit]}];
dime=Table[linfite[[i]][1]-linfite[[i]][0],{i,1,Length[datafit]}];
dimf=Table[linfitf[[i]][1]-linfitf[[i]][0],{i,1,Length[datafit]}];

Multicolumn[Table[{e[[i]],Binary[[i]],Plot[{linfitc[[i]][x],linfite[[i]][x],linfit[[i]][x],linfitb[[i]][x],linfitd[[i]][x],linfitf[[i]][x]},{x,-5,0.5},Epilog->{PointSize[Large],Point[Log[Flatten[datafit[[i]],1][[1;;21]]]]},PlotStyle->{Green,Red,Blue,Gray,Orange,Cyan},Frame->True,Axes->False,FrameLabel->{"Log((r^(-1)))","Log(N(r))"},FrameStyle->{{Large,Black},{Large,Black},{Large,Black},{Large,Black}},ImageSize->Medium,FrameTicks->{{{0,2,4,6,8},{}},{{-6,-5,-4,-3,-2,-1,0},{}}},PlotLegends->{"=(Subscript[D, 1 - 21])"NumberForm[dimc[[i]],3],"=(Subscript[D, 1 - 11])"NumberForm[dime[[i]],3],"=(Subscript[D, 12 - 21])"NumberForm[dim[[i]],3],"=(Subscript[D, 7 - 21])"NumberForm[dimb[[i]],3],"=(Subscript[D, 3 - 17])"NumberForm[dimd[[i]],3],"=(Subscript[D, 4 - 15])"NumberForm[dimf[[i]],3]}]},{i,{22}}],{1,2}]

Multicolumn[Table[{e[[i]],Binary[[i]],Plot[{linfitc[[i]][x],linfite[[i]][x],linfit[[i]][x],linfitb[[i]][x],linfitd[[i]][x],linfitf[[i]][x]},{x,-5,0.5},Epilog->{PointSize[Large],Point[Log[Flatten[datafit[[i]],1][[1;;21]]]]},PlotStyle->{Green,Red,Blue,Gray,Orange,Cyan},Frame->True,Axes->False,FrameLabel->{"Log((r^(-1)))","Log(N(r))"},FrameStyle->{{Large,Black},{Large,Black},{Large,Black},{Large,Black}},ImageSize->Medium,FrameTicks->{{{0,2,4,6,8},{}},{{-6,-5,-4,-3,-2,-1,0},{}}},PlotLegends->{"=(Subscript[D, 1 - 21])"NumberForm[dimc[[i]],3],"=(Subscript[D, 1 - 11])"NumberForm[dime[[i]],3],"=(Subscript[D, 12 - 21])"NumberForm[dim[[i]],3],"=(Subscript[D, 7 - 21])"NumberForm[dimb[[i]],3],"=(Subscript[D, 3 - 17])"NumberForm[dimd[[i]],3],"=(Subscript[D, 4 - 15])"NumberForm[dimf[[i]],3]}]},{i,{23,1,5,13}}],{4,2}]

(*Check the fit*)

Table[{e[[i]],Binary[[i]],Plot[{linfitc[[i]][x],linfite[[i]][x],linfit[[i]][x],linfitb[[i]][x],linfitd[[i]][x],linfitf[[i]][x]},{x,-5,0.5},Epilog->{PointSize[Large],Point[Log[Flatten[datafit[[i]],1][[1;;21]]]]},PlotStyle->{Green,Red,Blue,Gray,Orange,Cyan},Frame->True,Axes->False,FrameLabel->{"Log((r^(-1)))","Log(N(r))"},FrameStyle->{{Large,Black},{Large,Black},{Large,Black},{Large,Black}},ImageSize->Medium,FrameTicks->{{{0,2,4,6,8},{}},{{-6,-5,-4,-3,-2,-1,0},{}}},PlotLegends->{"=(Subscript[D, 1 - 21])"dimc[[i]],"=(Subscript[D, 1 - 11])"dime[[i]],"=(Subscript[D, 12 - 21])"dim[[i]],"=(Subscript[D, 7 - 21])"dimb[[i]],"=(Subscript[D, 3 - 17])"dimd[[i]],"=(Subscript[D, 4 - 15])"dimf[[i]]}]},{i,1,Length[datafit]}]