regmaxs/regmaxsn/scripts/utils/rotProfile2D.py

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()