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

filenamespos=FileNames["Data\\mapsmiddel0*.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>27.5&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+11&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<-(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+169,If[j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<53&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]-42&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i>-(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+117,0,datatot[[l,i,j]]-Min[datatot[[l]]]],0]},{i,1,Length[datatot[[l]]]},{j,1,101}],{l,1,Length[datatot]}];

dataprepnomask=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>25&&j*Cos[0.11496092050070644`]*0.5Sqrt[3]+Sin[0.11496092050070644`]*i<(i*Cos[0.11496092050070644`]-Sin[0.11496092050070644`]j*0.5Sqrt[3])*Sqrt[3]+15&&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,datatot[[l,i,j]]-Min[datatot[[l]]],0]},{i,1,Length[datatot[[l]]]},{j,1,101}],{l,1,Length[datatot]}];

ListDensityPlot[Flatten[dataprep[[1]],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

Length[dataprep[[1,1]]]

(*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.75*max[[i]],1,0]},{j,1,101},{l,1,101}],{i,1,Length[dataprep]}];

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

(*resulting images*)

img[[8]]

Binarize[img[[8]],0.7]

(*Compose the background*)

image3=Table[ImageCompose[ImageResize[ImageCrop[Graphics[Rectangle[{0,0},{1,1}]]],{360,360}],ImageResize[RemoveBackground[img[[i]]],1.5 ImageDimensions[img[[i]]]],Scaled[{0.5,0.325}]],{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]}]};

Binary

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}]

datafit=Table[{ParallelTable[{1/boxjes[[size]],Total[Sign/@(Total[#,2]&/@(ImageData/@Flatten[ImagePartition[Binary[[i]],boxjes[[size]]]]))]},{size,1,Length[boxjes],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]}];

Table[{e[[i]],linfit[[i]]},{i,1,Length[linfit]}];

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,10},{}},{{-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,{30,20,10,56}}],{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]}]