<commit_before># ! /usr/bin/python
# -*- coding: utf-8 -*-

"""
python src/histology_analyser.py -i ~/data/medical/data_orig/ja\
tra_mikro_data/Nejlepsi_rozliseni_nevycistene -t 6800 -cr 0 -1 100 300 100 300
"""

import sys
import os.path
path_to_script = os.path.dirname(os.path.abspath(__file__))
sys.path.append(os.path.join(path_to_script, "../extern/dicom2fem/src"))

import traceback

import logging
logger = logging.getLogger(__name__)

import argparse

from PyQt4.QtGui import QApplication

import numpy as np
import scipy.ndimage
import misc
from io3d import datareader
import csv

import skelet3d
import segmentation
import sed3 as se

import histology_analyser_gui as HA_GUI
from histology_report import HistologyReport
from skeleton_analyser import SkeletonAnalyser
import surface_measurement


class HistologyAnalyser:

    def __init__(self, data3d, metadata, threshold=-1, binaryClosing=-1,
                 binaryOpening=-1, nogui=True):
        self.data3d = data3d
        if 'voxelsize_mm' not in metadata.keys():
            metadata['voxelsize_mm'] = [1, 1, 1]
        self.metadata = metadata

        self.data3d_thr = None
        self.data3d_skel = None

        self.threshold = threshold
        
        if binaryClosing == -1:
            self.binaryClosing = 1
        else:
            self.binaryClosing = binaryClosing
        if binaryOpening == -1:
            self.binaryOpening = 1
        else:
            self.binaryOpening = binaryOpening

        self.nogui = nogui

        self.data3d_masked = None

    def get_voxelsize(self):
        return self.metadata['voxelsize_mm']

    def get_data3d(self):
        return self.data3d

    def get_data3d_segmented(self):
        return self.data3d_thr

    def get_data3d_skeleton(self):
        return self.data3d_skel

    def create_border_for_skeletonization(self, data3d, size=50):
        logger.debug('Generating border for skeletonization')
        data3d = data3d.copy()
        for i in range(size):
            logger.debug('iteration num '+str(i))

            new_shape = (data3d.shape[0]+2,
                         data3d.shape[1]+2,
                         data3d.shape[2]+2)
            logger.debug('new shape: '+str(new_shape))
            work_array = np.ones(new_shape, dtype = type(data3d[0][0][0]))

            # copy sides
            work_array[1:-1, 0, 1:-1] = data3d[:, 0, :]
            work_array[1:-1, -1, 1:-1] = data3d[:, -1, :]

            work_array[0, 1:-1, 1:-1] = data3d[0, :, :]
            work_array[-1, 1:-1, 1:-1] = data3d[-1, :, :]

            work_array[1:-1, 1:-1, 0] = data3d[:, :, 0]
            work_array[1:-1, 1:-1, -1] = data3d[:, :, -1]

            # copy corners sides
            work_array[1:-1, 0, 0] = data3d[:, 0, 0]
            work_array[1:-1, -1, 0] = data3d[:, -1, 0]
            work_array[1:-1, 0, -1] = data3d[:, 0, -1]
            work_array[1:-1, -1, -1] = data3d[:, -1, -1]

            work_array[0, 1:-1, 0] = data3d[0, :, 0]
            work_array[0, 1:-1, -1] = data3d[0, :, -1]
            work_array[-1, 1:-1, 0] = data3d[-1, :, 0]
            work_array[-1, 1:-1, -1] = data3d[-1, :, -1]

            work_array[0, 0, 1:-1] = data3d[0, 0, :]
            work_array[0, -1, 1:-1] = data3d[0, -1, :]
            work_array[-1, 0, 1:-1] = data3d[-1, 0, :]
            work_array[-1, -1, 1:-1] = data3d[-1, -1, :]

            # copy corners
            work_array[0, 0, 0] = data3d[0, 0, 0]
            work_array[0, 0, -1] = data3d[0, 0, -1]

            work_array[0, -1, 0] = data3d[0, -1, 0]
            work_array[0, -1, -1] = data3d[0, -1, -1]

            work_array[-1, 0, 0] = data3d[-1, 0, 0]
            work_array[-1, 0, -1] = data3d[-1, 0, -1]

            work_array[-1, -1, 0] = data3d[-1, -1, 0]
            work_array[-1, -1, -1] = data3d[-1, -1, -1]

            # erode
            erode_iterations = int(round(i/8))
            if erode_iterations>=1:
                work_array = scipy.ndimage.morphology.binary_erosion(work_array, border_value = 1, iterations = erode_iterations)

            # check if everything is eroded -> exit
            work_array[1:-1, 1:-1, 1:-1] = np.ones(data3d.shape, dtype = type(data3d[0][0][0]))
            eroded_sum = np.sum(np.sum(np.sum(work_array)))
            orig_sum = data3d.shape[0]*data3d.shape[1]*data3d.shape[2]

            if eroded_sum==orig_sum:
                eroded = True
            else:
                eroded = False

            # copy original data to center
            work_array[1:-1, 1:-1, 1:-1] = data3d

            data3d = work_array.copy()

            # if everything eroded -> exit
            if eroded:
                logger.debug('Everything eroded -> exiting early')
                break

        return data3d

    def data_to_binar(self):
        # ## Median filter
        filteredData = scipy.ndimage.filters.median_filter(self.data3d, size=2)

        # ## Segmentation
        data3d_thr = segmentation.vesselSegmentation(
            filteredData,  # self.data3d,
            segmentation=np.ones(filteredData.shape, dtype='int8'),
            threshold=self.threshold,
            inputSigma=0,  # 0.15,
            dilationIterations=2,
            nObj=1,
            biggestObjects=False,
            useSeedsOfCompactObjects=False,
            interactivity=not self.nogui,
            binaryClosingIterations=self.binaryClosing, # noqa 5,
            binaryOpeningIterations=self.binaryOpening # 1 # noqa
            )
        del(filteredData)

        # ## Zalepeni der
        scipy.ndimage.morphology.binary_fill_holes(data3d_thr,
                                                   output=data3d_thr)

        self.data3d_thr = data3d_thr

    def binar_to_skeleton(self):
        # create border with generated stuff for skeletonization
        #expanded_data = self.create_border_for_skeletonization(self.data3d_thr, size=50)
        expanded_data = self.data3d_thr

        expanded_skel = skelet3d.skelet3d(
            (expanded_data > 0).astype(np.int8)
        )

        # cut data shape back to original size
        border = (expanded_data.shape[0]-self.data3d_thr.shape[0])/2
        if border!=0:
            self.data3d_skel = expanded_skel[border:-border, border:-border, border:-border].copy()
        else:
            self.data3d_skel = expanded_skel

    def data_to_skeleton(self):
        self.data_to_binar()
        self.binar_to_skeleton()

    def data_to_statistics(self, guiUpdateFunction=None):
        self.stats = {}

        if (self.data3d_skel is None) or (self.data3d_thr is None):
            logger.debug('Skeleton was not generated!!! Generating now...')
            self.data_to_skeleton()

        # # add general info
        logger.debug('Computing general statistics...')
        if guiUpdateFunction is not None: guiUpdateFunction(0, 3, "Volumetric methods (Vv, Sv)")
        vs = self.metadata['voxelsize_mm']
        voxel_volume_mm3 = vs[0]*vs[1]*vs[2]
        volume_px = self.data3d.shape[0] \
            * self.data3d.shape[1] * self.data3d.shape[2]
        volume_mm3 = volume_px*voxel_volume_mm3

        # pouzij oblast zajmu
        if self.data3d_masked is None:
            self.data3d_masked = np.ones(self.data3d.shape, np.int8)
        self.data3d_masked[self.data3d_masked > 1] = 1
        used_volume_px = np.sum(np.sum(np.sum(self.data3d_masked)))
        used_volume_mm3 = used_volume_px*voxel_volume_mm3
        if guiUpdateFunction is not None: guiUpdateFunction(1, 3, "Volumetric methods (Vv, Sv)")

        vessel_volume_fraction = float(np.sum(np.sum(np.sum(
            self.data3d_thr)))) / float(used_volume_px)
        if guiUpdateFunction is not None: guiUpdateFunction(2, 3, "Volumetric methods (Vv, Sv)")
        
        sv = surface_measurement.surface_density(
            self.data3d_thr, self.metadata['voxelsize_mm'], self.data3d_masked)
        if guiUpdateFunction is not None: guiUpdateFunction(3, 3, "Volumetric methods (Vv, Sv)")

        info = {
            'voxel_size_mm': list(vs),
            'voxel_volume_mm3': float(voxel_volume_mm3),
            'shape_px': list(self.data3d.shape),
            'volume_px': float(volume_px),
            'volume_mm3': float(volume_mm3),
            'used_volume_px': float(used_volume_px),
            'used_volume_mm3': float(used_volume_mm3),
            'vessel_volume_fraction': float(vessel_volume_fraction),
            'surface_density': float(sv)
        }
        self.stats.update({'General': info})

        # # process skeleton to statistics
        logger.debug('Computing skeleton to statistics...')
        skan = SkeletonAnalyser(
            self.data3d_skel,
            volume_data=self.data3d_thr,
            voxelsize_mm=self.metadata['voxelsize_mm']
            )
        stats = skan.skeleton_analysis(guiUpdateFunction=guiUpdateFunction)
        # needed only by self.writeSkeletonToPickle()
        self.sklabel = skan.sklabel
        self.stats.update({'Graph': stats})
        
        # compute Nv statistic
        self.stats['General']['Nv'] = float(self.get_Nv())
        
    def get_Nv(self):
        logger.debug('Computing Nv...')
        
        # get lists of edges that are connected to nodes
        nodes = {}
        for key in self.stats['Graph']:
            edge = self.stats['Graph'][key]
            try:
                nodeIdA = edge['nodeIdA']
                if nodeIdA in nodes: nodes[nodeIdA] += [key] 
                else: nodes[nodeIdA] = [key]
            except Exception, e:
                logger.warning('get_Nv(): no nodeIdA')
            try:
                nodeIdB = edge['nodeIdB']
                if nodeIdB in nodes: nodes[nodeIdB] += [key] 
                else: nodes[nodeIdB] = [key]
            except Exception, e:
                logger.warning('get_Nv(): no nodeIdB')
        logger.debug('Read ' + str(len(nodes)) + ' nodes')
        
        # Get Pn (number of nodes with valence n)
        max_l = 0
        for n_key in nodes:
            max_l = max(max_l, len(nodes[n_key]))
        Pn = np.zeros(max_l+1)
        for n_key in nodes:
            Pn[len(nodes[n_key])] += 1
        logger.debug(str(Pn))
        
        # Get N(cap)
        Ncap = 0
        for n in range(len(Pn)):
            if n > 1: # dont count terminal nodes 
                Ncap += ( (n-2)/2.0 ) * Pn[n]
        
        # Get V(ref)
        Vref = self.stats['General']['used_volume_mm3']
        
        # Compute Nv
        Nv = ( Ncap / float(Vref) ) + 1
        logger.debug('Nv is ' + str(Nv))
        
        return Nv

    def showSegmentedData(self):
        skan = SkeletonAnalyser(
            self.data3d_skel,
            volume_data=self.data3d_thr,
            voxelsize_mm=self.metadata['voxelsize_mm'])
        data3d_nodes_vis = skan.sklabel.copy()
        del(skan)

        # edges
        data3d_nodes_vis[data3d_nodes_vis > 0] = 1
        # nodes and terminals
        data3d_nodes_vis[data3d_nodes_vis < 0] = 2

        if not self.nogui:
            pyed = se.sed3(
                self.data3d,
                seeds=(data3d_nodes_vis).astype(np.int8),
                contour=self.data3d_thr.astype(np.int8)
            )
            pyed.show()

    def writeStatsToYAML(self, filename='hist_stats.yaml'):<commit_msg>nv debug<commit_after># ! /usr/bin/python
# -*- coding: utf-8 -*-

"""
python src/histology_analyser.py -i ~/data/medical/data_orig/ja\
tra_mikro_data/Nejlepsi_rozliseni_nevycistene -t 6800 -cr 0 -1 100 300 100 300
"""

import sys
import os.path
path_to_script = os.path.dirname(os.path.abspath(__file__))
sys.path.append(os.path.join(path_to_script, "../extern/dicom2fem/src"))

import traceback

import logging
logger = logging.getLogger(__name__)

import argparse

from PyQt4.QtGui import QApplication

import numpy as np
import scipy.ndimage
import misc
from io3d import datareader
import csv

import skelet3d
import segmentation
import sed3 as se

import histology_analyser_gui as HA_GUI
from histology_report import HistologyReport
from skeleton_analyser import SkeletonAnalyser
import surface_measurement


class HistologyAnalyser:

    def __init__(self, data3d, metadata, threshold=-1, binaryClosing=-1,
                 binaryOpening=-1, nogui=True):
        self.data3d = data3d
        if 'voxelsize_mm' not in metadata.keys():
            metadata['voxelsize_mm'] = [1, 1, 1]
        self.metadata = metadata

        self.data3d_thr = None
        self.data3d_skel = None

        self.threshold = threshold
        
        if binaryClosing == -1:
            self.binaryClosing = 1
        else:
            self.binaryClosing = binaryClosing
        if binaryOpening == -1:
            self.binaryOpening = 1
        else:
            self.binaryOpening = binaryOpening

        self.nogui = nogui

        self.data3d_masked = None

    def get_voxelsize(self):
        return self.metadata['voxelsize_mm']

    def get_data3d(self):
        return self.data3d

    def get_data3d_segmented(self):
        return self.data3d_thr

    def get_data3d_skeleton(self):
        return self.data3d_skel

    def create_border_for_skeletonization(self, data3d, size=50):
        logger.debug('Generating border for skeletonization')
        data3d = data3d.copy()
        for i in range(size):
            logger.debug('iteration num '+str(i))

            new_shape = (data3d.shape[0]+2,
                         data3d.shape[1]+2,
                         data3d.shape[2]+2)
            logger.debug('new shape: '+str(new_shape))
            work_array = np.ones(new_shape, dtype = type(data3d[0][0][0]))

            # copy sides
            work_array[1:-1, 0, 1:-1] = data3d[:, 0, :]
            work_array[1:-1, -1, 1:-1] = data3d[:, -1, :]

            work_array[0, 1:-1, 1:-1] = data3d[0, :, :]
            work_array[-1, 1:-1, 1:-1] = data3d[-1, :, :]

            work_array[1:-1, 1:-1, 0] = data3d[:, :, 0]
            work_array[1:-1, 1:-1, -1] = data3d[:, :, -1]

            # copy corners sides
            work_array[1:-1, 0, 0] = data3d[:, 0, 0]
            work_array[1:-1, -1, 0] = data3d[:, -1, 0]
            work_array[1:-1, 0, -1] = data3d[:, 0, -1]
            work_array[1:-1, -1, -1] = data3d[:, -1, -1]

            work_array[0, 1:-1, 0] = data3d[0, :, 0]
            work_array[0, 1:-1, -1] = data3d[0, :, -1]
            work_array[-1, 1:-1, 0] = data3d[-1, :, 0]
            work_array[-1, 1:-1, -1] = data3d[-1, :, -1]

            work_array[0, 0, 1:-1] = data3d[0, 0, :]
            work_array[0, -1, 1:-1] = data3d[0, -1, :]
            work_array[-1, 0, 1:-1] = data3d[-1, 0, :]
            work_array[-1, -1, 1:-1] = data3d[-1, -1, :]

            # copy corners
            work_array[0, 0, 0] = data3d[0, 0, 0]
            work_array[0, 0, -1] = data3d[0, 0, -1]

            work_array[0, -1, 0] = data3d[0, -1, 0]
            work_array[0, -1, -1] = data3d[0, -1, -1]

            work_array[-1, 0, 0] = data3d[-1, 0, 0]
            work_array[-1, 0, -1] = data3d[-1, 0, -1]

            work_array[-1, -1, 0] = data3d[-1, -1, 0]
            work_array[-1, -1, -1] = data3d[-1, -1, -1]

            # erode
            erode_iterations = int(round(i/8))
            if erode_iterations>=1:
                work_array = scipy.ndimage.morphology.binary_erosion(work_array, border_value = 1, iterations = erode_iterations)

            # check if everything is eroded -> exit
            work_array[1:-1, 1:-1, 1:-1] = np.ones(data3d.shape, dtype = type(data3d[0][0][0]))
            eroded_sum = np.sum(np.sum(np.sum(work_array)))
            orig_sum = data3d.shape[0]*data3d.shape[1]*data3d.shape[2]

            if eroded_sum==orig_sum:
                eroded = True
            else:
                eroded = False

            # copy original data to center
            work_array[1:-1, 1:-1, 1:-1] = data3d

            data3d = work_array.copy()

            # if everything eroded -> exit
            if eroded:
                logger.debug('Everything eroded -> exiting early')
                break

        return data3d

    def data_to_binar(self):
        # ## Median filter
        filteredData = scipy.ndimage.filters.median_filter(self.data3d, size=2)

        # ## Segmentation
        data3d_thr = segmentation.vesselSegmentation(
            filteredData,  # self.data3d,
            segmentation=np.ones(filteredData.shape, dtype='int8'),
            threshold=self.threshold,
            inputSigma=0,  # 0.15,
            dilationIterations=2,
            nObj=1,
            biggestObjects=False,
            useSeedsOfCompactObjects=False,
            interactivity=not self.nogui,
            binaryClosingIterations=self.binaryClosing, # noqa 5,
            binaryOpeningIterations=self.binaryOpening # 1 # noqa
            )
        del(filteredData)

        # ## Zalepeni der
        scipy.ndimage.morphology.binary_fill_holes(data3d_thr,
                                                   output=data3d_thr)

        self.data3d_thr = data3d_thr

    def binar_to_skeleton(self):
        # create border with generated stuff for skeletonization
        #expanded_data = self.create_border_for_skeletonization(self.data3d_thr, size=50)
        expanded_data = self.data3d_thr

        expanded_skel = skelet3d.skelet3d(
            (expanded_data > 0).astype(np.int8)
        )

        # cut data shape back to original size
        border = (expanded_data.shape[0]-self.data3d_thr.shape[0])/2
        if border!=0:
            self.data3d_skel = expanded_skel[border:-border, border:-border, border:-border].copy()
        else:
            self.data3d_skel = expanded_skel

    def data_to_skeleton(self):
        self.data_to_binar()
        self.binar_to_skeleton()

    def data_to_statistics(self, guiUpdateFunction=None):
        self.stats = {}

        if (self.data3d_skel is None) or (self.data3d_thr is None):
            logger.debug('Skeleton was not generated!!! Generating now...')
            self.data_to_skeleton()

        # # add general info
        logger.debug('Computing general statistics...')
        if guiUpdateFunction is not None: guiUpdateFunction(0, 3, "Volumetric methods (Vv, Sv)")
        vs = self.metadata['voxelsize_mm']
        voxel_volume_mm3 = vs[0]*vs[1]*vs[2]
        volume_px = self.data3d.shape[0] \
            * self.data3d.shape[1] * self.data3d.shape[2]
        volume_mm3 = volume_px*voxel_volume_mm3

        # pouzij oblast zajmu
        if self.data3d_masked is None:
            self.data3d_masked = np.ones(self.data3d.shape, np.int8)
        self.data3d_masked[self.data3d_masked > 1] = 1
        used_volume_px = np.sum(np.sum(np.sum(self.data3d_masked)))
        used_volume_mm3 = used_volume_px*voxel_volume_mm3
        if guiUpdateFunction is not None: guiUpdateFunction(1, 3, "Volumetric methods (Vv, Sv)")

        vessel_volume_fraction = float(np.sum(np.sum(np.sum(
            self.data3d_thr)))) / float(used_volume_px)
        if guiUpdateFunction is not None: guiUpdateFunction(2, 3, "Volumetric methods (Vv, Sv)")
        
        sv = surface_measurement.surface_density(
            self.data3d_thr, self.metadata['voxelsize_mm'], self.data3d_masked)
        if guiUpdateFunction is not None: guiUpdateFunction(3, 3, "Volumetric methods (Vv, Sv)")

        info = {
            'voxel_size_mm': list(vs),
            'voxel_volume_mm3': float(voxel_volume_mm3),
            'shape_px': list(self.data3d.shape),
            'volume_px': float(volume_px),
            'volume_mm3': float(volume_mm3),
            'used_volume_px': float(used_volume_px),
            'used_volume_mm3': float(used_volume_mm3),
            'vessel_volume_fraction': float(vessel_volume_fraction),
            'surface_density': float(sv)
        }
        self.stats.update({'General': info})

        # # process skeleton to statistics
        logger.debug('Computing skeleton to statistics...')
        skan = SkeletonAnalyser(
            self.data3d_skel,
            volume_data=self.data3d_thr,
            voxelsize_mm=self.metadata['voxelsize_mm']
            )
        stats = skan.skeleton_analysis(guiUpdateFunction=guiUpdateFunction)
        # needed only by self.writeSkeletonToPickle()
        self.sklabel = skan.sklabel
        self.stats.update({'Graph': stats})
        
        # compute Nv statistic
        self.stats['General']['Nv'] = float(self.get_Nv())
        
    def get_Nv(self):
        logger.debug('Computing Nv...')
        
        # get lists of edges that are connected to nodes
        nodes = {}
        for key in self.stats['Graph']:
            edge = self.stats['Graph'][key]
            try:
                nodeIdA = edge['nodeIdA']
                if nodeIdA in nodes: nodes[nodeIdA] += [key] 
                else: nodes[nodeIdA] = [key]
            except Exception, e:
                logger.warning('get_Nv(): no nodeIdA')
            try:
                nodeIdB = edge['nodeIdB']
                if nodeIdB in nodes: nodes[nodeIdB] += [key] 
                else: nodes[nodeIdB] = [key]
            except Exception, e:
                logger.warning('get_Nv(): no nodeIdB')
        logger.debug('Read ' + str(len(nodes)) + ' nodes')
        
        # Get Pn (number of nodes with valence n)
        max_l = 0
        for n_key in nodes:
            max_l = max(max_l, len(nodes[n_key]))
        Pn = np.zeros(max_l+1)
        for n_key in nodes:
            Pn[len(nodes[n_key])] += 1
        logger.debug(str(Pn))
        
        # Get N(cap)
        Ncap = 0
        for n in range(len(Pn)):
            if n > 1: # dont count terminal nodes 
                Ncap += ( (n-2)/2.0 ) * Pn[n]
        
        # Get V(ref)
        Vref = self.stats['General']['used_volume_mm3']
        
        # Compute Nv
        Nv = ( Ncap / float(Vref) ) + 1
        import ipdb; ipdb.set_trace() #  noqa BREAKPOINT

        logger.debug('Nv is ' + str(Nv))
        
        return Nv

    def showSegmentedData(self):
        skan = SkeletonAnalyser(
            self.data3d_skel,
            volume_data=self.data3d_thr,
            voxelsize_mm=self.metadata['voxelsize_mm'])
        data3d_nodes_vis = skan.sklabel.copy()
        del(skan)

        # edges
        data3d_nodes_vis[data3d_nodes_vis > 0] = 1
        # nodes and terminals
        data3d_nodes_vis[data3d_nodes_vis < 0] = 2

        if not self.nogui:
            pyed = se.sed3(
                self.data3d,
                seeds=(data3d_nodes_vis).astype(np.int8),
                contour=self.data3d_thr.astype(np.int8)
            )
            pyed.show()

    def writeStatsToYAML(self, filename='hist_stats.yaml'):
{'fraction_non_alphanumeric': 0.12216197152052198, 'fraction_numerical': 0.033417096503989774, 'mean_word_length': 1.9625179574245788, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1750770', 'n_tokens_mistral': 8100, 'n_tokens_neox': 7729, 'n_words': 1509}