$HistoryLength=0;
SetDirectory[NotebookDirectory[]];
ClearAll;



(********************************************************************)
(*defined FLIM descriptor-1 for best quantification of the treatment*)

DescriptorDist[FLIMbef_,FLIMaft_]:=Module[{FLIMData,descriptor1,descriptor2},

FLIMData={FLIMbef,FLIMaft};
descriptor1=Flatten[FLIMData[[2,4]]/FLIMData[[2,1]]-FLIMData[[1,4]]/FLIMData[[1,1]]];

Return[{descriptor1}];
];

(********************************************************************)
(*defined FLIM descriptor-2 for best quantification of the treatment*)

DescriptorDist2[FLIMdataPop_]:=Module[{(*FLIMdataPop,*)descriptor1,descriptor2},

descriptor1=Flatten[FLIMdataPop[[4]]/FLIMdataPop[[1]]];
descriptor2=Flatten[FLIMdataPop[[2]]/100*FLIMdataPop[[4]]-FLIMdataPop[[1]]/100*FLIMdataPop[[3]]];

Return[{descriptor1,descriptor2}];
];


(*******************************************************)
(*setting of the treatment region on color coded image*)
(*Red locator region - mark treatment region*)
(*Blue locator region - mark any non-trreatment region*)
(********************************************************)

(*Import FLIM and AF intensity images*)
colorCodedImageName=FileNames["e04_s03_t17_eye_FLIM-1PE_after-therapy3_Bin7"~~___~~".tif"];
intensityImageName=FileNames[___~~"intensity"~~___~~".tif"];

Print[{colorCodedImageName,intensityImageName}];
FLIMDataTm=Import[colorCodedImageName[[1]],"Image"];
IntensityImage=Import[intensityImageName[[#]],"Image"]&/@{1,2};
IntensityImageData=ImageData[IntensityImage[[#]]]&/@{1,2};


Manipulate[
Quiet[

crop1={Min[xyCrop1[[#]],xyCrop2[[#]]]&/@{1,2},Max[xyCrop1[[#]],xyCrop2[[#]]]&/@{1,2}};
crop2={Min[xyCrop3[[#]],xyCrop4[[#]]]&/@{1,2},Max[xyCrop3[[#]],xyCrop4[[#]]]&/@{1,2}};

Show[FLIMDataTm,Frame->True,FrameLabel->{"x [pix]","y [pix]"},ImageSize->600,Epilog->{Transparent,EdgeForm[{Thickness[0.005],Red}],Rectangle@@(crop1-{1,0}),Transparent,EdgeForm[{Thickness[0.005],Blue}],Rectangle@@(crop2-{1,0})}]
]

,{{xyCrop1,{1,1}},{1,1},{512,512},{1,1},Locator,Enabled->True}
,{{xyCrop2,{200,200}},{1,1},{512,512},{1,1},Locator,Enabled->True}
,{{xyCrop3,{300,300}},{1,1},{512,512},{1,1},Locator,Enabled->True}
,{{xyCrop4,{512,512}},{1,1},{512,512},{1,1},Locator,Enabled->True}
]

(*************************************)
(*FLIM analysis before and after treatment*)
(************************************)

(*import the datasets of the FLIM parameters a1, a2, t1 and t2 calculated in each image pixel using SPC Image software (B&H) *)
(*before and after laser treatment*)
({{fileNamesAfter=FileNames["e04_s03_t17_eye_FLIM-1PE_after-therapy3"~~__~~".asc"];
fileNamesBefore=FileNames["e04_s03_t16_eye_FLIM-1PE_before-therapy2"~~__~~".asc"];
filesLength=Length[fileNamesAfter];
Print[{fileNamesAfter,fileNamesBefore}];

(*Importing and trimming data focused on the treatment site*)
FLIMDataAfter=Import[fileNamesAfter[[#]],"Table"]&/@Range[filesLength];
FLIMDataBefore=Import[fileNamesBefore[[#]],"Table"]&/@Range[filesLength];
(*Print[FLIMDataAfter[[1,200,200]]];*)

(*Show FLIM images od individual fitted parameter - color coded with Hue*)
Print[Table[ArrayPlot[FLIMDataAfter[[i]],ColorFunction->Hue,ImageSize->150],{i,1,filesLength}]];


(*********************************************************************************************************************************************)
(*analysis of the change in FLIM descriptor t2/a1 after the treatment (1-on the site on the treatment and 2-on the surroundings, as a control*)
(*********************************************************************************************************************************************)

FLIMDataAfterTher=FLIMDataAfter[[#]][[512-crop1[[2,2]];;512-crop1[[1,2]],crop1[[1,1]];;crop1[[2,1]]]]&/@Range[filesLength];
FLIMDataBeforeTher=FLIMDataBefore[[#]][[512-crop1[[2,2]];;512-crop1[[1,2]],crop1[[1,1]];;crop1[[2,1]]]]&/@Range[filesLength];
descriptorTher=DescriptorDist[FLIMDataBeforeTher,FLIMDataAfterTher];
descriptorsN=Length[descriptorTher];
Print[descriptorsN];

(*min and max values of individual FLIM parameters in the treatment region*)
xRange=Length[FLIMDataAfterTher[[1,1]]];
yRange=Length[FLIMDataAfterTher[[1]]];
minData=Min[Cases[Flatten[FLIMDataAfterTher[[#]]],Except[0]]]&/@Range[filesLength];
maxData=Max[FLIMDataAfterTher[[#]]]&/@Range[filesLength];
Print[{"minMax",{minData,maxData}}];


bins2={0.5,20};
histListTher=HistogramList[descriptorTher[[#]],{bins2[[#]]}]&/@Range[descriptorsN];
binsN=Length[histListTher[[#,1]]]&/@Range[descriptorsN];

(*FLIM signal change in the surrouding tissue*)
FLIMDataAfterSurr=FLIMDataAfter[[#]][[512-crop2[[2,2]];;512-crop2[[1,2]],crop2[[1,1]];;crop2[[2,1]]]]&/@Range[filesLength];
FLIMDataBeforeSurr=FLIMDataBefore[[#]][[512-crop2[[2,2]];;512-crop2[[1,2]],crop2[[1,1]];;crop2[[2,1]]]]&/@Range[filesLength];
descriptorSurr=DescriptorDist[FLIMDataBeforeSurr,FLIMDataAfterSurr];
histListSurr=HistogramList[descriptorSurr[[#]],{bins2[[#]]}]&/@Range[descriptorsN];
binsN=Length[histListSurr[[#,1]]]&/@Range[descriptorsN];

plotRange={{-5,7},{-300,300}};
bins={0.2,2};
descriptorNames={Style["(Subscript[t, 2])/(Subscript[a, 1]) (aft) - (Subscript[t, 2])/(Subscript[a, 1]) (bef)",20]};
HistColor={Blue, Green};
LegendPosition={0.9,0.8};
legendLabel={"treatment area","surroundings"};
descriptorsHist=Histogram[{descriptorTher[[#]],descriptorSurr[[#]]},{bins[[#]]},"Probability",ChartStyle->{Directive[HistColor[[1]],Opacity[.5]],Directive[HistColor[[2]],Opacity[.6]]},LabelStyle->{FontSize->20},
Epilog->{Inset[Style[legendLabel[[1]],HistColor[[1]],18],Scaled[{0.2, LegendPosition[[1]]}]],
Inset[Style[legendLabel[[2]],HistColor[[2]],18],Scaled[{0.2, LegendPosition[[2]]}]]},
ImageSize->500]&/@Range[descriptorsN];
descriptorsHist=Grid[{descriptorNames,descriptorsHist}];
Print[descriptorsHist];
site="comb";
Export[StringTake[fileNamesAfter[[1]],11]<>"_FLIMdeskriptors_Bef-and-Aft_"<>site<>".tif",descriptorsHist,ImageSize->500];



(***************************************************************************)
(*analysis of the parameters on the same site, before and after the therapy*)
(***************************************************************************)

descriptorTherAft=DescriptorDist2[FLIMDataAfterTher];
descriptorTherBef=DescriptorDist2[FLIMDataBeforeTher];
histListTherAft=HistogramList[descriptorTherAft[[#]],{bins2[[#]]}]&/@Range[2];
binsN=Length[histListTherAft[[#,1]]]&/@Range[2];
slopeTherAft=Table[(histListTherAft[[#,2,i+1]]-histListTherAft[[#,2,i]]),{i,1,binsN[[#]]-2}]&/@Range[2];
slopeTherAft=slopeTherAft[[#]]/Max[slopeTherAft[[#]]]&/@Range[2];
histListTherBef=HistogramList[descriptorTherBef[[#]],{bins2[[#]]}]&/@Range[2];
binsN=Length[histListTherBef[[#,1]]]&/@Range[2];
slopeTherBef=Table[(histListTherBef[[#,2,i+1]]-histListTherBef[[#,2,i]]),{i,1,binsN[[#]]-2}]&/@Range[2];
slopeTherBef=slopeTherBef[[#]]/Max[slopeTherBef[[#]]]&/@Range[2];

plotRange={{18,28},{-300,300}};
bins={0.2,2};
descriptorNames={Style["(Subscript[t, 2])/(Subscript[a, 1])",20],Style["(Subscript[a, 2])(Subscript[t, 2]) - (Subscript[a, 1])(Subscript[t, 1])",20]};

(*remove points on the right for best fitting*)
descriptorTherAftPop=Select[descriptorTherAft[[1]],#<30&];
dataPointsN=Length[descriptorTherAftPop];

(*Fitting treatment descriptors historgram with NormalDistribution function*)
edistAft=EstimatedDistribution[descriptorTherAftPop,NormalDistribution[μ,σ],{{μ,22},{σ,1}}];
edistBef=EstimatedDistribution[descriptorTherBef[[1]],NormalDistribution[μ,σ],{{μ,22},{σ,1}}];

edistAft2=EstimatedDistribution[descriptorTherAftPop,SkewNormalDistribution[μ,σ,α](*,{{μ,22},{σ,2},{α,2}}*)];
edistBef2=EstimatedDistribution[descriptorTherBef[[1]],SkewNormalDistribution[μ,σ,α](*,{{μ,22},{σ,2},{α,2}}*)];

(*Show FLIM descriptor parameters which will be used as an imput data for model function used for a real-time laser treatment quantification*)
(*parameter t2/a1 before and after the treatment - center value, histogran width & histogram skewness*)
Print[{{"dataPointsN","edistBef","edistAft","edistBef2-skewNormDist","edistAft2-skewNormDist"},dataPointsN,{edistBef[[1]],edistBef[[2]]},{edistAft[[1]],edistAft[[2]]},edistBef2[[3]],edistAft2[[3]]}];

(*Graphical representation of data*)

legendLabel2={"After treatment","before treatment"};
descriptorsHist=Show[Histogram[{descriptorTherAft[[#]],descriptorTherBef[[#]]},{bins[[#]]},"PDF",LabelStyle->{FontSize->20},ChartStyle->{Directive[Orange,Opacity[.6]],Directive[Blue,Opacity[.6]]},
Epilog->{Inset[Style[legendLabel2[[1]],Orange,18],Scaled[{0.7, LegendPosition[[1]]}]],
Inset[Style[legendLabel2[[2]],Blue,18],Scaled[{0.7, LegendPosition[[2]]}]]},
ImageSize->500,PlotRange->{{18,34},Automatic}],
Plot[PDF[edistAft,x],{x,15,40},PlotStyle->{Thick,Orange},PlotRange->All],
Plot[PDF[edistBef,x],{x,15,40},PlotStyle->{Thick,Blue},PlotRange->All],
Plot[PDF[edistAft2,x],{x,15,40},PlotStyle->{Thick,Dashed,Orange},PlotRange->All],
Plot[PDF[edistBef2,x],{x,15,40},PlotStyle->{Thick,Dashed,Blue},PlotRange->All]
]&/@Range[1];

Print[descriptorsHist];
site="comb";
(*Save histograms*)
Export[StringTake[fileNamesAfter[[1]],11]<>"_FLIMdeskriptors_Just-after"<>site<>".tif",descriptorsHist,ImageSize->500];


(*******************************************)
(*AF intensity analysis before and after treatment*)
(*******************************************)

IntensityImageDataCrop=IntensityImageData[[#]][[512-crop1[[2,2]];;512-crop1[[1,2]],crop1[[1,1]];;crop1[[2,1]]]]&/@Range[2];

minData=Min[Cases[Flatten[IntensityImageDataCrop[[#]]],Except[0]]]&/@Range[2];
maxData=Max[Cases[Flatten[IntensityImageDataCrop[[#]]],Except[0]]]&/@Range[2];
meanData=Mean[Cases[Flatten[IntensityImageDataCrop[[#]]],Except[0]]]&/@Range[2];
meanDataStDev=StandardDeviation[Cases[Flatten[IntensityImageDataCrop[[#]]],Except[0]]]&/@Range[2];

Print[{"meanData",meanData}];

(*Save FLIM and AF intensity descriptor parameters which will be used as an imput data for model function used for a real-time laser treatment quantification*)
exportList={meanData[[2]],meanData[[1]],meanDataStDev[[2]],meanDataStDev[[1]],edistBef[[1]],edistBef[[2]],edistAft[[1]],edistAft[[2]]};

(*change . to , that can be imported in excel*)
exportList2=StringReplace[ToString[exportList[[#]]],"."->","]&/@Range[Length[exportList]];
Export[StringTake[fileNamesAfter[[1]],39]<>"_dataForTherapyFunc.xslx",exportList2,"Table"];

(*reorganize data*)
FLIMData={FLIMDataBeforeTher,FLIMDataAfterTher};
histData1={Flatten[FLIMData[[#,1]]/100],Flatten[FLIMData[[#,3]]]}&/@Range[2];
histData2={Flatten[FLIMData[[#,2]]/100],Flatten[FLIMData[[#,4]]]}&/@Range[2];
histDataJoinBefore=Join[histData1[[1]],histData2[[1]],2];
histDataJoinAfter=Join[histData1[[2]],histData2[[2]],2];

(*Histogram[WeightedData[histData1[[2,2]],histData1[[2,1]]],50,"Probability"]
Histogram[WeightedData[histData2[[2,2]],histData2[[2,1]]],{1000,3000,100},"Probability"]
Histogram[WeightedData[histDataJoinAfter[[2]],histDataJoinAfter[[1]]],{0,4000,100},"Probability",ImageSize->400]*)

(*Graphical representation of data in the treatment region before and after treatment*)
(*sum of individual weighted FLIM decays components*)
FLIMdistributions=Graphics[Show[Histogram[{WeightedData[histDataJoinBefore[[2]],histDataJoinBefore[[1]]],WeightedData[histDataJoinAfter[[2]],histDataJoinAfter[[1]]]},{0,4000,20},"Probability",
ChartStyle->{Darker[Green],Red},PlotRange-> {{100,2500},{0,All}},Frame->True,FrameLabel->{"(Subscript[τ, FLIM]) [ps]","Probability density"//DisplayForm},LabelStyle->{FontSize->18},
PlotLabel->Style["τ distribution before and after the laser treatment\n",Bold,18],ImageSize->500]]];
Print[FLIMdistributions];

site="treat";


(*Save histograms*)
Export[StringTake[fileNamesAfter[[1]],11]<>"_FLIMdist_Bef-Aft_"<>site<>".tif",FLIMdistributions,ImageSize->1000];
Export[StringTake[fileNamesAfter[[1]],11]<>"_FLIMdist_Bef-Aft_"<>site<>".jpg",FLIMdistributions,ImageSize->1000];
FLIMData={{FLIMDataBeforeTher,FLIMDataAfterTher},{FLIMDataBeforeSurr,FLIMDataAfterSurr}};
}, {\[Placeholder]}})