Added DL2FF.py, experimental workflow to use DL peaksearch output and…

… run analysis.

FF_HEDM/src/CalcRadius.c

@@ -87,6 +87,10 @@ int main(int argc, char *argv[]){
     int RingNr = atoi(argv[2]);
     char aline[1000], *str, dummy[1000];
     fileParam = fopen(ParamFN,"r");
+    if (fileParam == NULL){
+		printf("Could not read file %s\n",ParamFN);
+		return 1;
+	}
     int LowNr = 1;
     char Folder[1024], FileStem[1024], fs[1024];
     int StartNr, EndNr, LayerNr;
@@ -205,6 +209,10 @@ int main(int argc, char *argv[]){
     sprintf(InputFile,"%s/PeakSearch/%s/Result_StartNr_%d_EndNr_%d_RingNr_%d.csv",Folder,FileStem,StartNr,EndNr,RingNr);
 	FILE *Infile;
 	Infile = fopen(InputFile,"r");
+    if (Infile == NULL){
+		printf("Could not read file %s\n",InputFile);
+		return 1;
+	}
 	fgets(aline,1000,Infile);
 	int counter = 0;
 	double **SpotsMat;

FF_HEDM/src/MergeOverlappingPeaks.c

@@ -163,6 +163,7 @@ static inline int ReadSortFiles (char OutFolderName[1024], char FileStem[1024],
 	sprintf(InFile,"%s/%s_%0*d_%d_PS.csv",OutFolderName,FileStem,Padding,FileNr,RingNr);
     FILE *infileread;
     infileread = fopen(InFile,"r");
+    if (infileread == NULL) printf("Could not read the input file %s\n",InFile);
     struct InputData *MyData;
     MyData = malloc(nOverlapsMaxPerImage*sizeof(*MyData));
     int counter = 0;
@@ -211,6 +212,10 @@ int main(int argc, char *argv[]){
     fflush(stdout);
     char aline[1000], *str, dummy[1000];
     fileParam = fopen(ParamFN,"r");
+    if (fileParam == NULL){
+		printf("Could not read file %s\n",ParamFN);
+		return 1;
+	}
     int LowNr = 1;
     char Folder[1024], FileStem[1024],*TmpFolder,fs[1024];
     int LayerNr;

utils/DL2FF.py

@@ -0,0 +1,123 @@
+import numpy as np
+from subprocess import check_call
+from glob import glob
+from os.path import expanduser
+import os
+from math import sqrt, degrees, acos
+
+peaksFN = 'center-of-mass-dl.csv'
+pFile = 'psAu.txt'
+nCPUs = 2
+
+def getValueFromParamFile(paramfn,searchStr,nLines=1,wordNr=1,nWords=1):
+	ret_list = []
+	nrLines = 0
+	f = open(paramfn,'r')
+	PSContents = f.readlines()
+	for line in PSContents:
+		if line.startswith(searchStr+' '):
+			line = line.replace('\t',' ')
+			line = line.replace('\n',' ')
+			words = line.split(' ')
+			words = [_f for _f in words if _f]
+			ret_list.append(words[wordNr:wordNr+nWords])
+			nrLines += 1
+			if (nrLines == nLines):
+				return ret_list
+	return ret_list
+
+def CalcEtaAngle(y,z):
+	alpha = -degrees(acos(z/sqrt(y*y+z*z))) if y > 0 else degrees(acos(z/sqrt(y*y+z*z)))
+	return alpha
+
+thisFolder = os.getcwd() + '/'
+paramFile = thisFolder + pFile
+binFolder = expanduser('~')+ '/opt/MIDAS/FF_HEDM/bin/'
+fileStem = getValueFromParamFile(paramFile,'FileStem')[0][0]
+px = float(getValueFromParamFile(paramFile,'px')[0][0])
+padding = int(getValueFromParamFile(paramFile,'Padding')[0][0])
+BC = [float(s) for s in getValueFromParamFile(paramFile,'BC',nWords=2)[0]]
+Rings = [int(s[0]) for s in getValueFromParamFile(paramFile,'RingThresh',nLines=100)]
+omegaStep = float(getValueFromParamFile(paramFile,'OmegaStep')[0][0])
+omegaStart = float(getValueFromParamFile(paramFile,'OmegaFirstFile')[0][0])
+
+width=float(getValueFromParamFile(paramFile,'Width')[0][0])
+folder = fileStem + '_Layer1_Analysis_Time_2020_06_11_12_55_10/'
+check_call('mkdir -p '+thisFolder+folder,shell=True)
+check_call(binFolder+'/GetHKLList ' + paramFile,shell=True)
+hkls = np.genfromtxt('hkls.csv',skip_header=1)
+check_call('cp '+thisFolder+'hkls.csv '+thisFolder+folder,shell=True)
+ringNrs=(set(list(hkls[:,4].flat)))
+ringNrs = [int(s) for s in ringNrs]
+ringSzs = np.zeros(len(ringNrs))
+for ctr,ringNr in enumerate(ringNrs):
+	for row in hkls:
+		if int(row[4]) == ringNr:
+			ringSzs[ctr] = float(row[10])
+peaksInfo = np.genfromtxt(peaksFN,delimiter=',',skip_header=1)
+nPeaks,ncols = peaksInfo.shape
+
+firstFrameNr = int(peaksInfo[0,1]) + 1
+lastFrameNr = int(peaksInfo[-1,1]) + 1
+head = 'SpotID IntegratedIntensity Omega(degrees) YCen(px) ZCen(px) IMax Radius(px) Eta(degrees) SigmaR SigmaEta\n'
+for ringNr in ringNrs:
+	folN = thisFolder+folder+'Ring'+str(ringNr)+'/Temp/'
+	check_call('mkdir -p '+folN,shell=True)
+	for fNr in range(firstFrameNr,lastFrameNr+1):
+		fn = folN + fileStem + '_1_' + str(fNr).zfill(padding) + '_' + str(ringNr) + '_PS.csv'
+		f = open(fn,'w')
+		f.write(head)
+		f.close()
+
+NrSpots = 0
+for peakNr in range(nPeaks):
+	thisY = peaksInfo[peakNr,2] - BC[0]
+	thisZ = peaksInfo[peakNr,3] - BC[1]
+	thisRad = sqrt(thisY**2+thisZ**2) * px
+	thisR = thisRad / px
+	thisEta = CalcEtaAngle(thisY,thisZ)
+	thisFrame = int(peaksInfo[peakNr,1]) + 1
+	thisOmega = omegaStart + (thisFrame-1)*omegaStep
+	for ctr,ringNr in enumerate(ringNrs):
+		tmpRingSz = ringSzs[ctr]
+		if thisRad < tmpRingSz + width and thisRad > tmpRingSz - width:
+			# now we write
+			folN = thisFolder+folder+'Ring'+str(ringNr)+'/Temp/'
+			fn = folN + fileStem + '_1_' + str(thisFrame).zfill(padding) + '_' + str(ringNr) + '_PS.csv'
+			f = open(fn,'r')
+			SpotID = len(f.readlines()) - 1
+			f.close()
+			outstr = str(SpotID) + ' 100 ' + str(thisOmega) + ' ' + str(thisY+BC[0]) + ' ' + str(thisZ+BC[1]) + ' 100 ' + str(thisR) + ' ' + str(thisEta) + ' 1 1\n'
+			f = open(fn,'a')
+			f.write(outstr)
+			f.close()
+			NrSpots += 1
+
+# Now run MergeOverlappingPeaks, then CalcRadius, then run a new FF analysis without peaksearch. A bit convoluted, but easiest way for now
+# First MergeOverlappingPeaks
+os.chdir(thisFolder+folder)
+for ringNr in ringNrs:
+	os.chdir(thisFolder+folder)
+	folN = thisFolder+folder+'Ring'+str(ringNr)
+	thisParamFile = thisFolder + folder + 'Layer1_' + pFile
+	check_call('cp '+paramFile+' '+thisParamFile,shell=True)
+	pf = open(thisParamFile,'a')
+	pf.write('Folder '+folN+'\n')
+	pf.write('RingToIndex '+str(ringNr)+'\n')
+	pf.write('RingNumbers '+str(ringNr)+'\n')
+	pf.write('LayerNr 1\n')
+	pf.close()
+	check_call(binFolder+'MergeOverlappingPeaks '+thisParamFile+' '+str(ringNr),shell=True)
+	check_call(binFolder+'CalcRadius '+thisParamFile+' '+str(ringNr),shell=True)
+	resFile=folN + '/PeakSearch/' + fileStem + '_1/Radius_StartNr_' + str(firstFrameNr) + '_EndNr_' + str(lastFrameNr) + '_RingNr_' + str(ringNr) + '.csv'
+	check_call('cp ' + resFile + ' ' + thisFolder + folder,shell=True)
+
+os.chdir(thisFolder)
+newParamFile = paramFile + '.New.txt'
+check_call('cp '+paramFile+' '+newParamFile,shell=True)
+f = open(newParamFile,'a')
+f.write('FolderName '+folder[:-1]+'\n')
+f.write('SeedFolder '+thisFolder+'\n')
+f.close()
+
+check_call(expanduser('~')+'/.MIDAS/MIDAS_V5_FarField_Layers '+pFile+'.New.txt 1 1 0 '+ str(nCPUs) + ' local sd',shell=True)