# -*- coding: utf-8 -*-
# **************************************************************************
# *
# * Authors: Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
# *
# * BCU, Centro Nacional de Biotecnologia, CSIC
# *
# *
# * This program is free software; you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation; either version 2 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program; if not, write to the Free Software
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
# * 02111-1307 USA
# *
# * All comments concerning this program package may be sent to the
# * e-mail address 'scipion@cnb.csic.es'
# *
# **************************************************************************
from os.path import basename
import numpy as np
from pwem.protocols import EMProtocol
from pyworkflow import BETA
from pyworkflow.protocol.params import PointerParam, FloatParam
from tomo.protocols import ProtTomoBase
import tomo.constants as const
[docs]class XmippProtConnectedComponents(EMProtocol, ProtTomoBase):
""" This protocol takes a set of coordinates and identifies connected
components among the picked particles."""
_label = 'connected components'
_devStatus = BETA
def __init__(self, **args):
EMProtocol.__init__(self, **args)
# --------------------------- DEFINE param functions ------------------------
def _defineParams(self, form):
form.addSection(label='Coordinates')
form.addParam('inputCoordinates', PointerParam, label="Coordinates",
pointerClass='SetOfCoordinates3D', help='Select the SetOfCoordinates3D.')
form.addParam('distance', FloatParam, label='Distance', help='Maximum radial distance (in voxels) between '
'particles to consider that they are in the same '
'connected component. Wizard returns three times '
'the box size of the input coordinates.')
# --------------------------- INSERT steps functions ----------------------
def _insertAllSteps(self):
self._insertFunctionStep('computeConnectedComponentsStep')
self._insertFunctionStep('createOutputStep')
# --------------------------- STEPS functions -------------------------------
[docs] def computeConnectedComponentsStep(self):
inputCoors = self.inputCoordinates.get()
coorSetList = []
tomoList = []
# Create a separate setOfCoordinates for each tomogram
for coor in inputCoors.iterCoordinates():
tomoName = coor.getVolName()
if tomoName not in tomoList:
tomoList.append(tomoName)
tomocoorset = self._createSetOfCoordinates3D(inputCoors.getPrecedents(), '_' + basename(tomoName))
coorSetList.append(tomocoorset)
idx = tomoList.index(tomoName)
coorSetList[idx].append(coor)
# For each tomogram coordinates, perform "connected components" logic
outputsetIndex = 0
self.outputSet = self._createSetOfCoordinates3D(inputCoors.getPrecedents())
self.outputSet.copyInfo(inputCoors)
self.outputSet.setBoxSize(inputCoors.getBoxSize())
self.outputSet.setSamplingRate(inputCoors.getSamplingRate())
for coorSet in coorSetList:
coorSet.write()
minDist = self.distance.get()
coorlist = []
for i, coor in enumerate(coorSet.iterCoordinates()):
coorlist.append([coor.getX(const.BOTTOM_LEFT_CORNER),
coor.getY(const.BOTTOM_LEFT_CORNER),
coor.getZ(const.BOTTOM_LEFT_CORNER)])
A = np.zeros([len(coorlist), len(coorlist)])
for j, coor1 in enumerate(coorlist):
for k, _ in enumerate(coorlist, start=j+1):
if k == len(coorlist):
break
else:
coor2 = coorlist[k]
if abs(coor1[0]-coor2[0]) <= minDist and abs(coor1[1]-coor2[1]) <= minDist \
and abs(coor1[2]-coor2[2]) <= minDist:
A[j, k] = 1
A[k, j] = 1
tomoNameshort = basename(coorSet.getFirstItem().getVolName())
np.savetxt(self._getExtraPath('adjacency_matrix_%s' % tomoNameshort), A)
D = np.diag(A.sum(axis=1))
np.savetxt(self._getExtraPath('degree_matrix_%s' % tomoNameshort), D)
L = D - A
np.savetxt(self._getExtraPath('laplacian_matrix_%s' % tomoNameshort), L)
vals, vecs = np.linalg.eig(L)
vals = vals.real
vecs = vecs.real
np.savetxt(self._getExtraPath('eigenvecs_matrix_%s' % tomoNameshort), vecs)
np.savetxt(self._getExtraPath('eigenvalues_matrix_%s' % tomoNameshort), vals)
vals0list = [i for i, x in enumerate(vals.real) if abs(x) < (1/(np.sqrt(len(coorSet)))*1e-3)]
listOfSets = [[] for x in range(len(vals0list))]
for k in range(len(coorSet)):
row = []
for j in vals0list:
row.append(vecs[k, j])
row = [abs(number)for number in row]
ixMax = np.argmax(row)
listOfSets[ixMax].append(k)
for ix, coorInd in enumerate(listOfSets):
if len(coorInd) != 0:
outputsetIndex += 1
for id, coor3D in enumerate(coorSet.iterItems()):
if id in coorInd:
coor3D.setGroupId(outputsetIndex)
self.outputSet.append(coor3D)
print('Connected components found: %d' % outputsetIndex)
[docs] def createOutputStep(self):
self._defineOutputs(outputSetOfCoordinates3D=self.outputSet)
self._defineSourceRelation(self.inputCoordinates, self.outputSet)
# --------------------------- INFO functions --------------------------------
def _validate(self):
validateMsgs = []
return validateMsgs
def _summary(self):
summary = []
summary.append("Maximum radial distance between particles in the same connected component: %d voxels"
% (self.distance.get()))
return summary
def _methods(self):
methods = []
methods.append("Coordinates grouped in connected component with a maximum radial distance of %d voxels."
% self.distance.get())
return methods