pattern_disctionary.py 11.3 KB
Newer Older
Jiri Borovec committed
1 2 3 4 5 6 7 8 9
"""


Copyright (C) 2015-2016 Jiri Borovec <jiri.borovec@fel.cvut.cz>
"""

import logging

from scipy import ndimage
10
from skimage import morphology, feature, filters
11
from scipy import ndimage as ndi
Jiri Borovec committed
12
import numpy as np
Jiri Borovec committed
13 14

import dataset_utils as data
Jiri Borovec committed
15
import pattern_weights as ptn_weight
Jiri Borovec committed
16

17
REINIT_PATTERN_COMPACT = True
Jiri Borovec committed
18

19 20 21 22 23 24

# TODO: init: Otsu threshold on sum over all input images -> WaterShade on distance
# TODO: init: sum over all input images and use it negative as distance for WaterShade


def initialise_atlas_random(im_size, label_max):
Jiri Borovec committed
25 26 27
    """ initialise atlas with random labels

    :param im_size: (w, h) size og image
28
    :param label_max: int, number of labels
Jiri Borovec committed
29 30
    :return: np.array<w, h>
    """
Jiri Borovec committed
31
    logging.debug('initialise atlas %s as random labeling', repr(im_size))
32
    nb_labels = label_max + 1
Jiri Borovec committed
33
    np.random.seed()  # reinit seed to have random samples even in the same time
34 35
    img_init = np.random.randint(1, nb_labels, im_size)
    return np.array(img_init, dtype=np.int)
Jiri Borovec committed
36 37


38
def initialise_atlas_mosaic(im_size, nb_labels, coef=1.):
39
    """ generate grids texture and into each rectangle plase a label,
Jiri Borovec committed
40 41 42 43 44 45
    each row contains all labels (permutation)

    :param im_size: (w, h) size og image
    :param max_lb: int, number of labels
    :return: np.array<w, h>
    """
Jiri Borovec committed
46
    logging.debug('initialise atlas %s as grid labeling', repr(im_size))
47
    nb_labels = int(nb_labels * coef)
Jiri Borovec committed
48
    np.random.seed()  # reinit seed to have random samples even in the same time
49 50 51
    block_size = np.ceil(np.array(im_size) / float(nb_labels))
    block = np.ones(block_size.astype(np.int))
    vec = range(1, nb_labels + 1) * int(np.ceil(coef))
Jiri Borovec committed
52
    logging.debug('block size is %s', repr(block.shape))
53 54 55
    for label in range(nb_labels):
        idx = np.random.permutation(vec)[:nb_labels]
        for k in range(nb_labels):
Jiri Borovec committed
56 57 58 59 60
            b = block.copy() * idx[k]
            if k == 0:
                row = b
            else:
                row = np.hstack((row, b))
61
        if label == 0: mosaic = row
Jiri Borovec committed
62
        else: mosaic = np.vstack((mosaic, row))
Jiri Borovec committed
63
    logging.debug('generated mosaic %s with labeling %s',
Jiri Borovec committed
64
                 repr(mosaic.shape), repr(np.unique(mosaic).tolist()))
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
    img_init = mosaic[:im_size[0], :im_size[1]]
    img_init = np.remainder(img_init, nb_labels)
    return np.array(img_init, dtype=np.int)


def initialise_atlas_otsu_watershed_2d(imgs, nb_labels, bg='none'):
    """ do some simple operations to get better initialisation
    1] sum over all images, 2] Otsu thresholding, 3] watershed

    :param imgs: [np.array<w, h>]
    :param nb_labels: int
    :param bg: str, set weather the Otsu backround sould be filled randomly
    :return: np.array<w, h>
    """
    logging.debug('initialise atlas for %i labels from %i images of shape %s '
                  'with Otsu-Watershed', nb_labels, len(imgs), repr(imgs[0].shape))
    img_sum = np.sum(np.asarray(imgs), axis=0) / float(len(imgs))
    img_gauss = filters.gaussian_filter(img_sum, 1)
    # http://scikit-image.org/docs/dev/auto_examples/plot_otsu.html
    thresh = filters.threshold_otsu(img_gauss)
    img_otsu = (img_gauss >= thresh)
    # http://scikit-image.org/docs/dev/auto_examples/plot_watershed.html
    img_dist = ndi.distance_transform_edt(img_otsu)
    local_maxi = feature.peak_local_max(img_dist, labels=img_otsu,
                                        footprint=np.ones((2, 2)))
    seeds = np.zeros_like(img_sum)
    seeds[local_maxi[:,0], local_maxi[:,1]] = range(1, len(local_maxi) + 1)
    labels = morphology.watershed(-img_dist, seeds)
    img_init = np.remainder(labels, nb_labels)
    if bg == 'rand':
        # add random labels on the potential backgound
        img_rand = np.random.randint(1, nb_labels, img_sum.shape)
        img_init[img_otsu == 0] = img_rand[img_otsu == 0]
    return img_init.astype(np.int)


def initialise_atlas_gauss_watershed_2d(imgs, nb_labels):
    """ do some simple operations to get better initialisation
    1] sum over all images, 2]watershed

    :param imgs: [np.array<w, h>]
    :param nb_labels: int
    :return: np.array<w, h>
    """
    logging.debug('initialise atlas for %i labels from %i images of shape %s '
                  'with Gauss-Watershed', nb_labels, len(imgs), repr(imgs[0].shape))
    img_sum = np.sum(np.asarray(imgs), axis=0) / float(len(imgs))
    img_gauss = filters.gaussian_filter(img_sum, 1)
    local_maxi = feature.peak_local_max(img_gauss, footprint=np.ones((2, 2)))
    seeds = np.zeros_like(img_sum)
    seeds[local_maxi[:,0], local_maxi[:,1]] = range(1, len(local_maxi) + 1)
    # http://scikit-image.org/docs/dev/auto_examples/plot_watershed.html
    labels = morphology.watershed(-img_gauss, seeds) # , mask=im_diff
    img_init = np.remainder(labels, nb_labels)
    return img_init.astype(np.int)
Jiri Borovec committed
120 121 122 123 124 125 126 127 128


def initialise_atlas_deform_original(atlas):
    """take the orginal atlas and use geometrical deformation
    to generate new deformed atlas

    :param atlas: np.array<w, h>
    :return: np.array<w, h>
    """
Jiri Borovec committed
129
    logging.debug('initialise atlas by deforming original one')
Jiri Borovec committed
130 131 132 133
    res = data.image_deform_elastic(atlas)
    return np.array(res, dtype=np.int)


Jiri Borovec committed
134
def reconstruct_samples(atlas, w_bins):
Jiri Borovec committed
135 136 137
    """ create reconstruction of binary images according given atlas and weights

    :param atlas: np.array<w, h> input atlas
Jiri Borovec committed
138
    :param w_bins: np.array<nb_imgs, nb_lbs>
Jiri Borovec committed
139 140
    :return: [np.array<w, h>] * nb_imgs
    """
Jiri Borovec committed
141 142 143
    # w_bins = np.array(weights)
    w_bin_ext = np.append(np.zeros((w_bins.shape[0], 1)), w_bins, axis=1)
    imgs = [None] * w_bins.shape[0]
Jiri Borovec committed
144 145
    for i, w in enumerate(w_bin_ext):
        imgs[i] = np.asarray(w)[np.asarray(atlas)]
Jiri Borovec committed
146
        assert atlas.shape == imgs[i].shape
Jiri Borovec committed
147 148 149
    return imgs


150 151
def prototype_new_pattern(imgs, imgs_reconst, diffs, atlas,
                          ptn_compact=REINIT_PATTERN_COMPACT):
Jiri Borovec committed
152 153 154 155
    """ estimate new pattern that occurs in input images and is not cover
    by any label in actual atlas, remove collision with actual atlas

    :param imgs: [np.array<w, h>] list of input images
156
    :param imgs_reconst: [np.array<w, h>] list of image reconstructions
Jiri Borovec committed
157 158 159 160
    :param diffs: [int] list of differences among input and reconstruct images
    :return: np.array<w, h> binary single pattern
    """
    id_max = np.argmax(diffs)
161
    # im_diff = np.logical_and(imgs[id_max] == True, imgs_reconst[id_max] == False)
Jiri Borovec committed
162
    # take just positive differences
163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
    im_diff = (imgs[id_max] - imgs_reconst[id_max]) > 0
    if ptn_compact:  # WaterShade
        logging.debug('.. reinit pattern using WaterShade')
        # im_diff = morphology.opening(im_diff, morphology.disk(3))
        # http://scikit-image.org/docs/dev/auto_examples/plot_watershed.html
        dist = ndi.distance_transform_edt(im_diff)
        local_maxi = feature.peak_local_max(dist, indices=False, labels=im_diff,
                                            footprint=np.ones((3, 3)))
        labels = morphology.watershed(-dist, ndi.label(local_maxi)[0], mask=im_diff)
    else:
        logging.debug('.. reinit pattern as major component')
        im_diff = morphology.closing(im_diff, morphology.disk(1))
        labels = None
    # find largest connected component
    img_ptn = data.extract_image_largest_element(im_diff, labels)
    # ptn_size = np.sum(ptn) / float(np.product(ptn.shape))
    # if ptn_size < 0.01:
    #     logging.debug('new patterns was too small %f', ptn_size)
    #     ptn = data.extract_image_largest_element(im_diff)
182
    img_ptn = (img_ptn == True)
183
    # img_ptn = np.logical_and(img_ptn == True, atlas == 0)
184 185 186 187 188
    return img_ptn


def insert_new_pattern(imgs, imgs_reconst, atlas, label,
                       ptn_compact=REINIT_PATTERN_COMPACT):
Jiri Borovec committed
189 190 191
    """ with respect to atlas empty spots inset new patterns

    :param imgs: [np.array<w, h>] list of input images
192
    :param imgs_reconst: [np.array<w, h>]
Jiri Borovec committed
193
    :param atlas: np.array<w, h>
194
    :param label: int
Jiri Borovec committed
195 196 197
    :return: np.array<w, h> updated atlas
    """
    # count just positive difference
198 199
    diffs = [np.sum((im - im_rc) > 0) for im, im_rc in zip(imgs, imgs_reconst)]
    im_ptn = prototype_new_pattern(imgs, imgs_reconst, diffs, atlas, ptn_compact)
Jiri Borovec committed
200
    # logging.debug('new im_ptn: {}'.format(np.sum(im_ptn) / np.prod(im_ptn.shape)))
Jiri Borovec committed
201
    # plt.imshow(im_ptn), plt.title('im_ptn'), plt.show()
202
    atlas[im_ptn == True] = label
203
    logging.debug('area of new pattern %i is %i', label, np.sum(atlas == label))
Jiri Borovec committed
204 205 206
    return atlas


207 208
def reinit_atlas_likely_patterns(imgs, w_bins, atlas, label_max=None,
                                 ptn_compact=REINIT_PATTERN_COMPACT):
Jiri Borovec committed
209 210
    """ walk and find all all free labels and try to reinit them by new patterns

Jiri Borovec committed
211
    :param label_max:
Jiri Borovec committed
212 213 214 215 216 217
    :param imgs: [np.array<w, h>] list of input images
    :param w_bins: np.array<nb_imgs, nb_lbs>
    :param atlas: np.array<w, h>
    :return: np.array<w, h>, np.array<nb_imgs, nb_lbs>
    """
    # find empty patterns
Jiri Borovec committed
218 219
    if label_max is None:
        label_max = max(np.max(atlas), w_bins.shape[1])
Jiri Borovec committed
220
    else:
Jiri Borovec committed
221
        logging.debug('compare w_bin %s to max %i', repr(w_bins.shape), label_max)
Jiri Borovec committed
222
        for i in range(w_bins.shape[1], label_max):
Jiri Borovec committed
223
            logging.debug('adding disappeared weigh column %i', i)
Jiri Borovec committed
224
            w_bins = np.append(w_bins, np.zeros((w_bins.shape[0], 1)), axis=1)
225 226
    # w_bin_ext = np.append(np.zeros((w_bins.shape[0], 1)), w_bins, axis=1)
    # logging.debug('IN > sum over weights: %s', repr(np.sum(w_bin_ext, axis=0)))
Jiri Borovec committed
227
    # add one while indexes does not cover 0 - bg
Jiri Borovec committed
228
    logging.debug('total nb labels: %i', label_max)
Jiri Borovec committed
229
    atlas_new = atlas.copy()
230 231
    for label in range(1, label_max + 1):
        w_index = label - 1
Jiri Borovec committed
232
        w_sum = np.sum(w_bins[:, w_index])
233
        logging.debug('reinit. label: %i with weight sum %i', label, w_sum)
Jiri Borovec committed
234 235
        if w_sum > 0:
            continue
236 237 238
        imgs_reconst = reconstruct_samples(atlas_new, w_bins)
        atlas_new = insert_new_pattern(imgs, imgs_reconst, atlas_new, label,
                                       ptn_compact)
239
        # logging.debug('w_bins before: %i', np.sum(w_bins[:, w_index]))
Jiri Borovec committed
240
        lim_repopulate = 100. / np.prod(atlas_new.shape)
241 242
        w_bins[:, w_index] = ptn_weight.weights_label_atlas_overlap_threshold(
                                        imgs, atlas_new, label, lim_repopulate)
Jiri Borovec committed
243
        logging.debug('w_bins after: %i', np.sum(w_bins[:, w_index]))
Jiri Borovec committed
244
    return atlas_new, w_bins
Jiri Borovec committed
245 246


247
def atlas_split_indep_ptn(atlas, label_max):
Jiri Borovec committed
248 249 250
    """ split  independent patterns labeled equally

    :param atlas: np.array<w, h>
251
    :param label_max: int
Jiri Borovec committed
252 253
    :return:
    """
Jiri Borovec committed
254
    patterns = []
255 256 257
    for label in np.unique(atlas):
        labeled, nb_objects = ndimage.label(atlas == label)
        logging.debug('for label %i detected #%i', label, nb_objects)
Jiri Borovec committed
258 259
        ptn = [(labeled == j) for j in np.unique(labeled)]
        # skip the largest one assuming to be background
Jiri Borovec committed
260 261
        patterns += sorted(ptn, key=lambda x: np.sum(x), reverse=True)[1:]
    patterns = sorted(patterns, key=lambda x: np.sum(x), reverse=True)
Jiri Borovec committed
262
    logging.debug('list of all areas %s', repr([np.sum(p) for p in patterns]))
Jiri Borovec committed
263
    atlas_new = np.zeros(atlas.shape, dtype=np.int)
264 265 266
    # take just label_max largest elements
    for i, ptn in enumerate(patterns[:label_max]):
        label = i + 1
267
        # logging.debug('pattern #%i area %i', lb, np.sum(ptn))
268
        atlas_new[ptn] = label
Jiri Borovec committed
269 270 271 272 273

    # plt.figure()
    # plt.subplot(121), plt.imshow(atlas), plt.colorbar()
    # plt.subplot(122), plt.imshow(atlas_new), plt.colorbar()
    # plt.show()
Jiri Borovec committed
274
    logging.debug('atlas unique %s', repr(np.unique(atlas_new)))
Jiri Borovec committed
275
    return atlas_new