(* 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\\mapsgroot-*.dat"];

filenamespos=FileNames["Data\\mapsgroot0*.dat"];

e=Flatten[{Table[ToExpression[StringTake[filenames[[i]],{-11,-7}]],{i,1,Length[filenames]}],Table[ToExpression[StringTake[filenamespos[[i]],{-10,-7}]],{i,1,Length[filenamespos]}]},1];

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

(*energy of the maps*)

e

(*prep the images*)

dataprep=ParallelTable[Table[{i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3],j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i,If[j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>34&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+21&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<-(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+275,If[j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<86&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]-82&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>-(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+172||j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<113&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]-30&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>-(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+224||j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<60&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]-31&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>-(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+118||j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<60&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]-136&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>-(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+224,0,datatot[[l,i,j]]-Min[datatot[[l]]]],0]},{i,1,Length[datatot[[l]]]},{j,1,177}],{l,1,Length[datatot]}];

ListDensityPlot[Flatten[dataprep[[8]],1]]

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,161},{l,1,177}],{i,1,Length[dataprep]}];

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

(*resulting images*)

ImageCrop[RemoveBackground[img[[8]]]]

Binarize[img[[8]],0.7]

(*Compose the background*)

image3=Table[ImageCompose[ImageResize[ImageCrop[Graphics[Rectangle[{0,0},{1,1}]]],{500,500}],ImageResize[RemoveBackground[img[[i]]],1.35ImageDimensions[img[[i]]]],Scaled[{0.49,0.46}]],{i,1,Length[datatot]}];

(*Binarize the image*)

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

Skel

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

Table[360/boxjes[[i]],{i,1,21}]

(*translate the image*)

binarytrans=Table[ParallelTable[ImageTransformation[Binary[[i]],#+{l,p} &],{l,-0.10,0.08,0.02},{p,-0.10,0.08,0.02}],{i,1,Length[Binary]}];
binaryskel=Table[ParallelTable[ImageTransformation[Skel[[i]],#+{l,p} &],{l,-0.10,0.08,0.02},{p,-0.10,0.08,0.02}],{i,1,Length[Skel]}];

(*Take the minimum of boxes that is needed to cover the set*)

datafit2skel=Table[{Table[{1/boxjes[[size]],Min[Total[Sign[(Total[#,{2,3}]&/@Map[ImageData,ParallelTable[Flatten[ImagePartition[Flatten[binaryskel[[i]]][[m]],{boxjes[[size]]}]],{m,1,100}],{2}])],{2}]]},{size,1,18,1}]},{i,1,7,1}];

datafit22skel=Table[{Table[{1/boxjes[[size]],Min[Total[Sign[(Total[#,{2,3}]&/@Map[ImageData,ParallelTable[Flatten[ImagePartition[Flatten[binaryskel[[i]]][[m]],{boxjes[[size]]}]],{m,1,100}],{2}])],{2}]]},{size,1,18,1}]},{i,8,15,1}];

datafit3trans=Table[{Table[{1/boxjes[[size]],Min[Total[Sign[(Total[#,{2,3}]&/@Map[ImageData,ParallelTable[Flatten[ImagePartition[Flatten[binarytrans[[i]]][[m]],{boxjes[[size]]}]],{m,1,100}],{2}])],{2}]]},{size,1,18,1}]},{i,1,7,1}];

datafit33trans=Table[{Table[{1/boxjes[[size]],Min[Total[Sign[(Total[#,{2,3}]&/@Map[ImageData,ParallelTable[Flatten[ImagePartition[Flatten[binarytrans[[i]]][[m]],{boxjes[[size]]}]],{m,1,100}],{2}])],{2}]]},{size,1,18,1}]},{i,8,15,1}];

datafit2totskel=Flatten[{datafit2skel,datafit22skel},1];

datafit3tottrans=Flatten[{datafit3trans,datafit33trans},1];

linfit2=Table[LinearModelFit[Log[Flatten[datafit2totskel[[i]],1][[12;;18]]],x,x],{i,1,15}];
linfitb2=Table[LinearModelFit[Log[Flatten[datafit2totskel[[i]],1][[7;;18]]],x,x],{i,1,15}];
linfitc2=Table[LinearModelFit[Log[Flatten[datafit2totskel[[i]],1][[1;;18]]],x,x],{i,1,15}];
linfitd2=Table[LinearModelFit[Log[Flatten[datafit2totskel[[i]],1][[3;;17]]],x,x],{i,1,15}];
linfite2=Table[LinearModelFit[Log[Flatten[datafit2totskel[[i]],1][[1;;11]]],x,x],{i,1,15}];
linfitf2=Table[LinearModelFit[Log[Flatten[datafit2totskel[[i]],1][[4;;15]]],x,x],{i,1,15}];
dim2=Table[linfit2[[i]][1]-linfit2[[i]][0],{i,1,15}];
dimb2=Table[linfitb2[[i]][1]-linfitb2[[i]][0],{i,1,15}];
dimc2=Table[linfitc2[[i]][1]-linfitc2[[i]][0],{i,1,15}];
dimd2=Table[linfitd2[[i]][1]-linfitd2[[i]][0],{i,1,15}];
dime2=Table[linfite2[[i]][1]-linfite2[[i]][0],{i,1,15}];
dimf2=Table[linfitf2[[i]][1]-linfitf2[[i]][0],{i,1,15}];

linfit3=Table[LinearModelFit[Log[Flatten[datafit3tottrans[[i]],1][[12;;18]]],x,x],{i,1,15}];
linfitb3=Table[LinearModelFit[Log[Flatten[datafit3tottrans[[i]],1][[7;;18]]],x,x],{i,1,15}];
linfitc3=Table[LinearModelFit[Log[Flatten[datafit3tottrans[[i]],1][[1;;18]]],x,x],{i,1,15}];
linfitd3=Table[LinearModelFit[Log[Flatten[datafit3tottrans[[i]],1][[3;;17]]],x,x],{i,1,15}];
linfite3=Table[LinearModelFit[Log[Flatten[datafit3tottrans[[i]],1][[1;;11]]],x,x],{i,1,15}];
linfitf3=Table[LinearModelFit[Log[Flatten[datafit3tottrans[[i]],1][[4;;15]]],x,x],{i,1,15}];
dim3=Table[linfit3[[i]][1]-linfit3[[i]][0],{i,1,15}];
dimb3=Table[linfitb3[[i]][1]-linfitb3[[i]][0],{i,1,15}];
dimc3=Table[linfitc3[[i]][1]-linfitc3[[i]][0],{i,1,15}];
dimd3=Table[linfitd3[[i]][1]-linfitd3[[i]][0],{i,1,15}];
dime3=Table[linfite3[[i]][1]-linfite3[[i]][0],{i,1,15}];
dimf3=Table[linfitf3[[i]][1]-linfitf3[[i]][0],{i,1,15}];

Table[{e[[i]],Skel[[i]],Plot[{linfitc2[[i]][x],linfite2[[i]][x],linfit2[[i]][x],linfitb2[[i]][x],linfitd2[[i]][x],linfitf2[[i]][x]},{x,-5,0.5},Epilog->{PointSize[Large],Point[Log[Flatten[datafit2totskel[[i]],1][[1;;18]]]]},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,10},{}},{{-6,-5,-4,-3,-2,-1,0},{}}},PlotLegends->{"= (Subscript[D, 1 - 18])"dimc2[[i]],"= (Subscript[D, 1 - 11])"dime2[[i]],"= (Subscript[D, 12 - 18])"dim2[[i]],"= (Subscript[D, 7 - 18])"dimb2[[i]],"= (Subscript[D, 3 - 17])"dimd2[[i]],"=(Subscript[D, 4 - 15])"dimf2[[i]]}]},{i,1,15}]

Multicolumn[Table[{e[[i]],Skel[[i]],Plot[{linfitc2[[i]][x],linfite2[[i]][x],linfit2[[i]][x],linfitb2[[i]][x],linfitd2[[i]][x],linfitf2[[i]][x]},{x,-5,0.5},Epilog->{PointSize[Large],Point[Log[Flatten[datafit2totskel[[i]],1][[1;;18]]]]},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 - 18])"dimc2[[i]],"= (Subscript[D, 1 - 11])"dime2[[i]],"= (Subscript[D, 12 - 18])"dim2[[i]],"= (Subscript[D, 7 - 18])"dimb2[[i]],"= (Subscript[D, 3 - 17])"dimd2[[i]],"=(Subscript[D, 4 - 15])"dimf2[[i]]}]},{i,{6,4,2,10}}],{4,2}]

Table[{e[[i]],Binary[[i]],Plot[{linfitc3[[i]][x],linfite3[[i]][x],linfit3[[i]][x],linfitb3[[i]][x],linfitd3[[i]][x],linfitf3[[i]][x]},{x,-5,0.5},Epilog->{PointSize[Large],Point[Log[Flatten[datafit3tottrans[[i]],1][[1;;18]]]]},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,10},{}},{{-6,-5,-4,-3,-2,-1,0},{}}},PlotLegends->{"= (Subscript[D, 1 - 18])"dimc3[[i]],"= (Subscript[D, 1 - 11])"dime3[[i]],"= (Subscript[D, 12 - 18])"dim3[[i]],"= (Subscript[D, 7 - 18])"dimb3[[i]],"= (Subscript[D, 3 - 17])"dimd3[[i]],"=(Subscript[D, 4 - 15])"dimf3[[i]]}]},{i,1,15}]

Multicolumn[Table[{e[[i]],Binary[[i]],Plot[{linfitc3[[i]][x],linfite3[[i]][x],linfit3[[i]][x],linfitb3[[i]][x],linfitd3[[i]][x],linfitf3[[i]][x]},{x,-5,0.5},Epilog->{PointSize[Large],Point[Log[Flatten[datafit3tottrans[[i]],1][[1;;18]]]]},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 - 18])"dimc3[[i]],"= (Subscript[D, 1 - 11])"dime3[[i]],"= (Subscript[D, 12 - 18])"dim3[[i]],"= (Subscript[D, 7 - 18])"dimb3[[i]],"= (Subscript[D, 3 - 17])"dimd3[[i]],"=(Subscript[D, 4 - 15])"dimf3[[i]]}]},{i,{6,4,2,10}}],{4,2}]