Debug

FF_HEDM/GetPackages.sh

@@ -19,15 +19,15 @@ if [ ! -d ${dirThis}/NLOPT ]; then # NLOPT INSTALL
 	make -j8 install
 fi
 
-if [ ! -d ${dirThis}/NLOPTShared ]; then # NLOPT SHARED INSTALL
-	cd $dirThis
-	echo $(pwd)
-	rm -rf nlopt-2.4.2
-	tar -xzf nlopt.tar.gz
-	cd nlopt-2.4.2
-	./configure --prefix=${dirThis}/NLOPTShared --enable-shared
-	make -j8 install
-fi
+# if [ ! -d ${dirThis}/NLOPTShared ]; then # NLOPT SHARED INSTALL
+# 	cd $dirThis
+# 	echo $(pwd)
+# 	rm -rf nlopt-2.4.2
+# 	tar -xzf nlopt.tar.gz
+# 	cd nlopt-2.4.2
+# 	./configure --prefix=${dirThis}/NLOPTShared --enable-shared
+# 	make -j8 install
+# fi
 
 if [ ! -d ${dirThis}/swift ]; then # SWIFT
 	cd $dirThis

FF_HEDM/Makefile

@@ -31,8 +31,7 @@ endif
 SRCDIR=src/
 BINDIR=bin/
 
-all: help bindircheck calibrant imagemax \
-	fittiltbclsdsample fitposorstrains peaksfitting \
+all: help bindircheck calibrant  imagemax fitposorstrains peaksfitting fittiltbclsdsample \
 	mergeoverlaps calcradius findsaturatedpx genmediandark fitgrain tiff2ge \
 	mergerings fittiltx fitwedge hkls indexer bindata processgrains graintracking\
 	mapmultdetectors matchgrains detectormapper mergemultiplescans \

FF_HEDM/src/CalibrantOMP.c

@@ -405,11 +405,13 @@ void CalcFittedMean(int nIndices, int *NrEachIndexBin, int **Indices, double *Av
 			Rmi = Rs[j] - Rstep/2;
 			Rma = Rs[j] + Rstep/2;
 			CalcPeakProfileParallel(IndicesThis,NrIndicesThis,idxThis,Average,Rmi,Rma,EtaMi,EtaMa,ybc,zbc,px,NrPixels, &RetVal);
+			// printf("%lf ",RetVal);
 			PeakShape[j] = RetVal;
 			if (RetVal != 0){
 				AllZero = 0;
 			}
 		}
+		// printf("\n");
 		for (j=0;j<NrPtsForFit;j++){
 			Etas[j]=EtaMean[idxThis];
 		}
@@ -426,6 +428,7 @@ void CalcFittedMean(int nIndices, int *NrEachIndexBin, int **Indices, double *Av
 		}else{
 			Rfit = 0;
 		}
+		// printf("\t\t%lf %lf %lf %lf %lf\n",Rmin,Rmax,Rfit,ytr,ztr);
 		RMean[idxThis] = Rfit;
 		free(NrPts);
 		free(Rm);
@@ -597,6 +600,7 @@ void CorrectTiltSpatialDistortion(int nIndices, double MaxRad, double *YMean, do
 		Diffs[i] = Diff;
 		MeanDiff += Diff;
 		RadOuts[i] = Rcorr;
+		// printf("%lf %lf %lf %lf %lf %lf %lf\n",Rad,Lsd,XYZ[0],XYZ[1],XYZ[2],YMean[i],ZMean[i]);
 	}
 	MeanDiff /= nIndices;
 	double StdDiff2=0;
@@ -755,10 +759,11 @@ int fileReader (FILE *f,char fn[], int dType, int NrPixels, double *returnArr, c
 		int j,k;
 	    hid_t file;
 	    herr_t status, status_n;                             
-	    hid_t dataset;  
+	    hid_t dataset;
 		hid_t dataspace;
 	    hsize_t dims[3];
 	    int ndims;
+		printf("%s\n",fn);
 	    file = H5Fopen(fn,H5F_ACC_RDONLY, H5P_DEFAULT);
 	    dataset = H5Dopen(file, DATASETNAME,H5P_DEFAULT);
 	    dataspace = H5Dget_space(dataset);
@@ -838,12 +843,27 @@ int main(int argc, char *argv[])
 	int HeadSize = 8192;
 	int dType = 1;
 	char GapFN[4096], BadPxFN[4096];
+	char darkDatasetName[4096], dataDatasetName[4096];
+	sprintf(darkDatasetName,"exchange/dark");
+	sprintf(dataDatasetName,"exchange/data");
     while (fgets(aline,1000,fileParam)!=NULL){
 		str = "FileStem ";
         LowNr = strncmp(aline,str,strlen(str));
         if (LowNr==0){
             sscanf(aline,"%s %s", dummy, fn);
             continue;
+        }
+		str = "darkDataset ";
+        LowNr = strncmp(aline,str,strlen(str));
+        if (LowNr==0){
+            sscanf(aline,"%s %s", dummy, darkDatasetName);
+            continue;
+        }
+		str = "dataDataset ";
+        LowNr = strncmp(aline,str,strlen(str));
+        if (LowNr==0){
+            sscanf(aline,"%s %s", dummy, dataDatasetName);
+            continue;
         }
 		str = "Folder ";
         LowNr = strncmp(aline,str,strlen(str));
@@ -1338,6 +1358,11 @@ int main(int argc, char *argv[])
 			TotFrames+=nFrames;
 			fclose(fp);
 		} else{
+			printf("Reading HDF5.\n");
+			printf("%s\n",FileName);
+			// sprintf(dname,"%s",darkDatasetName);
+			// printf("%s\n",dname);
+			// return 1;
 			dname = "exchange/dark";
 			rc = fileReader(fd,FileName,dType,NrPixelsY*NrPixelsZ,DarkFile,dname);
 			MakeSquare(NrPixels,NrPixelsY,NrPixelsZ,DarkFile,DarkFile2);
@@ -1360,6 +1385,7 @@ int main(int argc, char *argv[])
 			}
 			printf("Dark file read.\n");
 			for (j=0;j<(NrPixels*NrPixels);j++)AverageDark[j]=DarkFile2[j];
+			// sprintf(dname,"%s",dataDatasetName);
 			dname = "exchange/data";
 			rc = fileReader(fd,FileName,dType,NrPixelsY*NrPixelsZ,Image,dname);
 			MakeSquare(NrPixels,NrPixelsY,NrPixelsZ,Image,Image2);
@@ -1427,6 +1453,7 @@ int main(int argc, char *argv[])
 				RMean2[countr] = RMean[i];
 				EtaMean2[countr] = EtaMean[i];
 				IdealTtheta2[countr] = IdealTtheta[i];
+				// printf("%lf %lf %lf %lf\n",RMean[i],IdealR[i],EtaMean[i],IdealTtheta[i]);
 				countr++;
 			}
 		}

FF_HEDM/src/FitPosOrStrainsOMP.c

@@ -39,6 +39,7 @@
 #include <sys/mman.h>
 #include <omp.h>
 #include <libgen.h>
+#include <unistd.h>
 
 #define deg2rad 0.0174532925199433
 #define rad2deg 57.2957795130823
@@ -51,6 +52,18 @@
 #define TestBit(A,k)  (A[(k/32)] &   (1 << (k%32)))
 #define MAXNOMEGARANGES 2000
 
+int
+CalcDiffractionSpots(double Distance,
+	double ExcludePoleAngle,
+	double OmegaRanges[MAXNOMEGARANGES][2],
+	int NoOfOmegaRanges,
+	double **hkls,
+	int n_hkls,
+	double BoxSizes[MAXNOMEGARANGES][4],
+	int *nTspots,
+	double OrientMatr[3][3],
+	double **TheorSpots);
+
 // For detector mapping!
 extern int BigDetSize;
 extern int *BigDetector;
@@ -669,7 +682,7 @@ double FitErrorsPosT(double x[12],int nSpotsComp,double spotsYZOIn[nSpotsComp][9
 }
 
 static inline
-double FitErrorsOrientStrains(double x[9],int nSpotsComp,double spotsYZO[nSpotsComp][12],int nhkls,double hklsIn[nhkls][7],
+double FitErrorsOrientStrains(double x[9],int nSpotsComp,double spotsYZO[nSpotsComp][9],int nhkls,double hklsIn[nhkls][7],
 					 double Lsd,double Wavelength,int nOmeRanges,double OmegaRanges[nOmeRanges][2],
 					 double BoxSizes[nOmeRanges][4],double MinEta,double wedge,double chi, double Pos[3])
 {
@@ -1488,7 +1501,7 @@ int main(int argc, char *argv[])
 	FILE *hklf = fopen(hklfn,"r");
 	if (hklf == NULL){
 		printf("Could not read the hkl file. Exiting.\n");
-		return;
+		return 1;
 	}
 	fgets(aline,1000,hklf);
 	int h,kt,l,Rnr, nhkls=0;
@@ -1511,7 +1524,7 @@ int main(int argc, char *argv[])
 		}
 	}
 	fclose(hklf);
-	if (nOmeRanges != nBoxSizes){printf("Number of omega ranges and number of box sizes don't match. Exiting!\n");return;}
+	if (nOmeRanges != nBoxSizes){printf("Number of omega ranges and number of box sizes don't match. Exiting!\n");return 1;}
 	double MargOme=0.01,MargPos=Rsample,MargPos2=Rsample/2,MargOme2=2,chi=0;
 	int thisRowNr;
 	# pragma omp parallel for num_threads(numProcs) private(thisRowNr) schedule(dynamic)

FF_HEDM/src/ForwardSimulationCompressed.c

@@ -1363,14 +1363,13 @@ main(int argc, char *argv[])
 			// }
 			// printf("\n");
 			OrientMat2Euler(OM,EulerThis);
-			Euler2OrientMat(EulerThis,OM);
-			// for (i=0;i<3;i++){
-			// 	for (j=0;j<3;j++){
-			// 		printf("%lf ",OM[i][j]);
-			// 	}
-			// 	printf("\n");
-			// }
-			// printf("\n");
+			double OMT[9];
+			Euler2OrientMat(EulerThis,OMT);
+			for (i=0;i<3;i++){
+				for (j=0;j<3;j++){
+					OM[i][j] = OMT[i*3+j];
+				}
+			}
 			// Calculate the spots now.
 			CalcDiffrSpots_Furnace(hklsOut,OM,Lsd,Wavelength,TheorSpots,&nTspots);
 			// For each spot, calculate displacement, calculate tilt and wedge effect.

FF_HEDM/src/GetMisorientation.c

@@ -24,6 +24,9 @@
 #define deg2rad 0.0174532925199433
 #define rad2deg 57.2957795130823
 #define EPS 1e-9
+static inline double sin_cos_to_angle (double s, double c){return (s >= 0.0) ? acos(c) : 2.0 * M_PI - acos(c);}
+static inline double cosd(double x){return cos(deg2rad*x);}
+static inline double sind(double x){return sin(deg2rad*x);}
 
 static inline void normalizeQuat(double quat[4]){
 	double norm = sqrt(quat[0]*quat[0]+quat[1]*quat[1]+quat[2]*quat[2]+quat[3]*quat[3]);

FF_HEDM/src/IndexerOMP.c

@@ -24,6 +24,7 @@
 #include <sys/types.h>
 #include <omp.h>
 #include <libgen.h>
+#include <unistd.h>
 
 static void
 check (int test, const char * message, ...)
@@ -62,6 +63,17 @@ check (int test, const char * message, ...)
 RealType *ObsSpotsLab;
 int n_spots = 0;
 
+int
+CalcDiffractionSpots(double Distance,
+	double ExcludePoleAngle,
+	double OmegaRanges[MAX_N_OMEGARANGES][2],
+	int NoOfOmegaRanges,
+	double **hkls,
+	int n_hkls,
+	double BoxSizes[MAX_N_OMEGARANGES][4],
+	int *nTspots,
+	double OrientMatr[3][3],
+	double **TheorSpots);
 
 // hkls to use
 double hkls[MAX_N_HKLS][7];

FF_HEDM/src/PeaksFittingOMPZarr.c

@@ -551,11 +551,11 @@ check (int test, const char * message, ...)
     }
 }
 
-void main(int argc, char *argv[]){
+int main(int argc, char *argv[]){
     double start_time = omp_get_wtime();
 	if (argc < 5){
 		printf("Usage: %s DataFile blockNr nBlocks numProcs (optional)ResultFolder\n",argv[0]);
-		return;
+		return 1;
 	}
 	char *DataFN = argv[1];
 	int blockNr = atoi(argv[2]);
@@ -594,6 +594,7 @@ void main(int argc, char *argv[]){
 	double BadPxIntensity = 0;
     char *resultFolder=NULL, dummy[2048];
 	int maskLoc = -1;
+	int doPeakFit = 1;
     while ((zip_stat_index(arch, count, 0, finfo)) == 0) {
         if (strstr(finfo->name,"exchange/data/.zarray")!=NULL){
             s = calloc(finfo->size + 1, sizeof(char));
@@ -714,6 +715,18 @@ void main(int argc, char *argv[]){
             free(arr);
             free(data);
         }
+        if (strstr(finfo->name,"measurement/process/scan_parameters/doPeakFit/0")!=NULL){
+            arr = calloc(finfo->size + 1, sizeof(char)); 
+            fd = zip_fopen_index(arch, count, 0);
+            zip_fread(fd, arr, finfo->size);
+            dsize = sizeof(int);
+            data = (char*)malloc((size_t)dsize);
+            dsize = blosc1_decompress(arr,data,dsize);
+            doPeakFit = *(int *)&data[0];
+			printf("doPeakFit: %d\n",doPeakFit);
+            free(arr);
+            free(data);
+        }
         if (strstr(finfo->name,"analysis/process/analysis_parameters/ResultFolder/0")!=NULL){
             arr = calloc(finfo->size + 1, sizeof(char)); 
             fd = zip_fopen_index(arch, count, 0);
@@ -1533,12 +1546,38 @@ void main(int argc, char *argv[]){
 				}
 			}
 			printf("nPeaks %d\n",nPeaks);
-			double retVal;
-			int rc = Fit2DPeaks(nPeaks,NrPixelsThisRegion,z,UsefulPixels,MaximaValues,MaximaPositions,IntegratedIntensity,IMAX,YCEN,ZCEN,Rads,Etass,Ycen,Zcen,Thresh,NrPx,OtherInfo,NrPixels,&retVal);
+			double retVal=0;
+			int rc = 0;
+			// If we don't want to fit, we can just compute weighted center of mass, but first put nPeaks =1;
+			// We need {IntegratedIntensity}, {IMAX}, {YCEN}, {ZCEN}, {Rads}, {Etass}, sigmaR (0), sigmaEta (0), {NrPx}, returnCode (0), retVal (0)
+			// \tBG\tSigmaGR\tSigmaLR\tSigmaGEta\tSigmaLEta\tMU\n" All of these will be 0
+			// OtherInfo is set already, will be 0 and will populate the other values in the line above.
+			if (doPeakFit == 0){
+				nPeaks = 1;
+				IMAX[0] = MaximaValues[0];
+				NrPx[0] = NrPixelsThisRegion;
+				YCEN[0] = 0;
+				ZCEN[0] = 0;
+				IntegratedIntensity[0] = 0;
+				for (i=0;i<NrPixelsThisRegion;i++){
+
+				}
+				IntegratedIntensity[0];
+				YCEN[0];
+				ZCEN[0];
+				Rads[0] = CalcNorm2(YCEN[0],ZCEN[0]);
+				Etass[0] = CalcEtaAngle(YCEN[0],ZCEN[0]);
+
+			} else{
+				rc = Fit2DPeaks(nPeaks,NrPixelsThisRegion,z,UsefulPixels,MaximaValues,MaximaPositions,IntegratedIntensity,IMAX,YCEN,ZCEN,Rads,Etass,Ycen,Zcen,Thresh,NrPx,OtherInfo,NrPixels,&retVal);
+			}
 			for (i=0;i<nPeaks;i++){
-				fprintf(outfilewrite,"%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t",(SpotIDStart+i),IntegratedIntensity[i],Omega,-YCEN[i]+Ycen,ZCEN[i]+Zcen,IMAX[i],Rads[i],Etass[i]);
+				fprintf(outfilewrite,"%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t",(SpotIDStart+i),IntegratedIntensity[i],
+					Omega,-YCEN[i]+Ycen,ZCEN[i]+Zcen,IMAX[i],Rads[i],Etass[i]);
 				for (j=0;j<2;j++) fprintf(outfilewrite, "%f\t",OtherInfo[8*i+6+j]);
-				fprintf(outfilewrite,"%d\t%d\t%d\t%d\t%d\t%f\t%f\t%f\t%d\t%lf",NrPx[i],NrPixelsThisRegion,nPeaks,MaximaPositions[i*2+0],MaximaPositions[i*2+1],(double)MaximaPositions[i*2+0]+YCEN[i]-Ycen,(double)MaximaPositions[i*2+1]-ZCEN[i]-Zcen,MaximaValues[i],rc,retVal);
+				fprintf(outfilewrite,"%d\t%d\t%d\t%d\t%d\t%f\t%f\t%f\t%d\t%lf",NrPx[i],NrPixelsThisRegion,
+					nPeaks,MaximaPositions[i*2+0],MaximaPositions[i*2+1],(double)MaximaPositions[i*2+0]+YCEN[i]-Ycen,
+					(double)MaximaPositions[i*2+1]-ZCEN[i]-Zcen,MaximaValues[i],rc,retVal);
 				for (j=0;j<6;j++) fprintf(outfilewrite, "\t%f",OtherInfo[8*i+j]);
 				fprintf(outfilewrite,"\n");
 			}

FF_HEDM/src/ProcessGrainsZarr.c

@@ -42,6 +42,9 @@
 #define rad2deg 57.2957795130823
 static inline double sin_cos_to_angle (double s, double c){return (s >= 0.0) ? acos(c) : 2.0 * M_PI - acos(c);}
 
+inline
+double GetMisOrientation(double quat1[4], double quat2[4], double axis[3], double *Angle,int SGNr);
+
 static inline
 void OrientMat2Euler(double m[3][3],double Euler[3])
 {
@@ -125,6 +128,9 @@ FreeMemMatrix(double **mat,int nrows)
     free(mat);
 }
 
+inline 
+void OrientMat2Quat(double OrientMat[9], double Quat[4]);
+
 static inline
 int
 FindInternalAnglesTwins(int nrIDs, int *IDs, int *IDsPerGrain,
@@ -241,11 +247,22 @@ QuatToOrientMat(
     OrientMat[8] = 1 - 2*(Q1_2+Q2_2);
 }
 
+inline void
+CalcStrainTensorFableBeaudoin(double LatCin[6],double LatticeParameterFit[6],
+	double Orient[3][3], double StrainTensorSample[3][3]);
+
+inline int
+StrainTensorKenesei(int nspots,double **SpotsInfo, double Distance, double wavelength,
+		double StrainTensorSample[3][3], int **IDHash,
+		double *dspacings, int nRings, int startSpotMatrix, double **SpotMatrix, double *RetVal,
+		double StrainTensorInput[3][3]);
+
+
 int main(int argc, char *argv[])
 {
 	if (argc != 2){
 		printf("Usage: ProcessGrains ZarrZip\n");
-		return;
+		return 1;
 	}
 	clock_t start, end;
     double diftotal;
@@ -767,7 +784,7 @@ int main(int argc, char *argv[])
 		}
 		if (retval == 0){
 			printf("Did not read correct hash table for IDs. Exiting\n");
-			return;
+			return 1;
 		}
 		FinalMatrix[nGrains][0] = GrainIDThis;
 		for (j=0;j<21;j++){

gui/ff_asym.py

@@ -1040,7 +1040,7 @@ def replot():
 # Main function
 root = Tk.Tk()
 root.wm_title("FF display v0.2 Dt. 2024/02/10 [email protected]")
-figur = Figure(figsize=(17,8),dpi=100)
+figur = Figure(figsize=(15,6),dpi=100)
 canvas = FigureCanvasTkAgg(figur,master=root)
 a = figur.add_subplot(121,aspect='equal')
 b = figur.add_subplot(122,aspect='equal')