moved rotProfile2D from GJEMS

regmaxsn/scripts/utils/rotProfile2D.py

@@ -0,0 +1,122 @@
+from regmaxsn.core.SWCTransforms import SWCRotate, ArgGenIterator, objFun
+import multiprocessing as mp
+import numpy as np
+from itertools import product
+import json
+from matplotlib import pyplot as plt
+import seaborn as sns
+plt.ion()
+import os
+from regmaxsn.core.matplotlibRCParams import mplPars
+
+
+sns.set(rc=mplPars)
+# ----------------------------------------------------------------------------------------------------------------------
+homeFolder = "/home/aj"
+dirPath = os.path.join(homeFolder,'DataAndResults/morphology/OriginalData/borstHSN2D/')
+#
+expNames = [
+            'HSN-fluoro01.CNG2D',
+            'HSN-fluoro01.CNG2DRandRot',
+            # 'HSN-fluoro01.CNG2DRandRot1',
+            # 'HSN-fluoro01.CNG2DRandRot2',
+            # 'HSN-fluoro01.CNG2DRandRot3',
+            # 'HSN-fluoro01.CNG2DRandRot4',
+            # 'HSN-fluoro01.CNG2DRandRot5',
+            # 'HSN-fluoro01.CNGNoiseStd42DRandRot',
+            # 'HSN-fluoro01.CNGNoiseStd42DRandRot1',
+            # 'HSN-fluoro01.CNGNoiseStd42DRandRot2',
+            # 'HSN-fluoro01.CNGNoiseStd42DRandRot3',
+            # 'HSN-fluoro01.CNGNoiseStd42DRandRot4',
+            # 'HSN-fluoro01.CNGNoiseStd42DRandRot5',
+            # 'HSN-fluoro02.CNG2D',
+            # 'HSN-fluoro03.CNG',
+            # 'HSN-fluoro04.CNG',
+            # 'HSN-fluoro05.CNG2D',
+            # 'HSN-fluoro06.CNG',
+            # 'HSN-fluoro07.CNG',
+            # 'HSN-fluoro08.CNG',
+            # 'HSN-fluoro09.CNG',
+            # 'HSN-fluoro10.CNG',
+            ]
+
+
+refInd = 0
+resDir = os.path.join(homeFolder,'DataAndResults/morphology/directPixelBased/rotProfile2D/')
+if not os.path.isdir(resDir):
+    os.mkdir(resDir)
+# ----------------------------------------------------------------------------------------------------------------------
+
+# gridSizes = [40.0, 20.0, 10.0]
+gridSizes = [40.0]
+# gridSizes = [20.0]
+# gridSizes = [10.0]
+bounds = [[0, 0], [0, 0], [-np.pi / 6, np.pi / 6]]
+
+# bounds = [[0, 0], [-np.pi / 3, np.pi / 3], [-np.pi / 3, np.pi / 3]]
+
+
+rotMinRes = np.deg2rad(1).round(4)
+nCPU = 7
+
+refSWC = os.path.join(dirPath, expNames[refInd] + '.swc')
+SWC2Align = os.path.join(dirPath, expNames[1] + '.swc')
+SWCRotCenter = np.loadtxt(SWC2Align)[:, 2:5].mean(axis=0)
+data = np.loadtxt(SWC2Align)[:, 2:5] - SWCRotCenter
+maxDist = np.linalg.norm(data, axis=1).max()
+pool = mp.Pool(processes=nCPU)
+with sns.axes_style('whitegrid'):
+    # fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(10, 8))
+    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 5.6))
+
+with open(os.path.join(dirPath, expNames[1] + '.json')) as fle:
+    pars = json.load(fle)
+    actualSol = -np.rad2deg(pars['angles'][2])
+
+
+SWCDatas = [SWCRotate(refSWC, SWC2Align, x) for x in gridSizes]
+noChangeVals = [objFun(([0, 0, 0], data)) for data in SWCDatas]
+stepSizes = [2 * x / maxDist for x in gridSizes]
+
+for gridInd, gridSize in enumerate(gridSizes):
+
+    # print('Gridsize:' + str(gridSize))
+    stepSize = stepSizes[gridInd]
+    # print('Stepsize: ' + str(np.rad2deg(stepSize)))
+    bounds = np.array(bounds)
+    boundsRoundedUp = np.sign(bounds) * np.ceil(np.abs(bounds) / stepSize) * stepSize
+    possiblePts1D = [np.round(np.arange(x[0], x[1] + 0.9 * stepSize, stepSize), 3).tolist() for x in boundsRoundedUp]
+    possiblePts3D = np.round(list(product(*possiblePts1D)), 3).tolist()
+    argGen = ArgGenIterator(possiblePts3D, SWCDatas[gridInd])
+    funcVals = pool.map(objFun, argGen)
+    outFile = os.path.join(resDir, expNames[refInd] + expNames[1] + str(int(gridSize)) + 'rot.json')
+    with open(outFile, 'w') as fle:
+        json.dump({'possiblePoints3D': possiblePts3D, 'funcVals': funcVals}, fle)
+    with sns.axes_style('whitegrid'):
+
+        # ax1 = ax[gridInd]
+        ax1 = ax
+        ax1.plot(np.rad2deg(possiblePts1D[2]), funcVals, 'b-o', lw=6, ms=12)
+        ybounds = (0.1 * np.floor(min(funcVals)*10), 0.1 * np.ceil(max(funcVals)*10))
+        ax1.set_ylim(ybounds)
+        ax1.plot([actualSol] * 2, ybounds, 'r-', lw=6)
+        ax1.set_ylabel(r'Dissimilarity')
+        # ax1.set_title(r'GridSize=' + str(int(gridSize)) + '$\mu$m')
+    if not gridInd:
+        bestSol = possiblePts3D[np.argmin(funcVals)]
+    else:
+        minimum = min(funcVals)
+        if minimum > noChangeVals[gridInd]:
+            bestSol = [0, 0, 0]
+        else:
+        # if True:
+            minimzers = [y for x, y in enumerate(possiblePts3D) if funcVals[x] == minimum]
+            prevVals = [objFun((x, SWCDatas[gridInd - 1])) for x in minimzers]
+            bestSol = minimzers[np.argmin(prevVals)]
+
+    bestVal = objFun((bestSol, SWCDatas[gridInd]))
+    print(bestSol)
+
+ax1.set_xlabel(r'rotation about Z in degrees')
+fig.tight_layout()
+fig.canvas.draw()