Python scipy.ndimage.measurements.center_of_mass() Examples

def test_kdakernel1(testdata):
    """
    COMs of KDA kernel maps should match the foci fed in (assuming focus isn't
    masked out and spheres don't overlap).
    Test on 1mm template.
    """
    id_ = 'pain_01.nidm-1'
    kern = kernel.KDAKernel(r=4, value=1)
    ma_maps = kern.transform(testdata['dset'].coordinates, testdata['dset'].masker)

    ijk = testdata['dset'].coordinates.loc[testdata['dset'].coordinates['id'] == id_,
                                           ['i', 'j', 'k']]
    ijk = np.squeeze(ijk.values.astype(int))
    kern_data = ma_maps[0].get_fdata()
    com = np.array(center_of_mass(kern_data)).astype(int).T
    com = np.squeeze(com)
    assert np.array_equal(ijk, com) 

def find_centroid(region):
    """
    Finds an approximate centroid for a region that is within the region.
    
    Parameters
    ----------
    region : np.ndarray(shape=(m, n), dtype='bool')
        mask of the region.

    Returns
    -------
    i, j : tuple(int, int)
        2d index within the region nearest the center of mass.
    """

    x, y = center_of_mass(region)
    w = np.argwhere(region)
    i, j = w[np.argmin(np.linalg.norm(w - (x, y), axis=1))]
    return i, j 

def test_kdakernel2(testdata):
    """
    COMs of KDA kernel maps should match the foci fed in (assuming focus isn't
    masked out and spheres don't overlap).
    Test on 2mm template.
    """
    id_ = 'pain_01.nidm-1'
    kern = kernel.KDAKernel(r=4, value=1)
    ma_maps = kern.transform(testdata['dset'].coordinates, testdata['dset'].masker)

    ijk = testdata['dset'].coordinates.loc[testdata['dset'].coordinates['id'] == id_,
                                           ['i', 'j', 'k']]
    ijk = np.squeeze(ijk.values.astype(int))
    kern_data = ma_maps[0].get_fdata()
    com = np.array(center_of_mass(kern_data)).astype(int).T
    com = np.squeeze(com)
    assert np.array_equal(ijk, com) 

def test_mkdakernel2(testdata):
    """
    COMs of MKDA kernel maps should match the foci fed in (assuming focus isn't
    masked out and spheres don't overlap).
    Test on 2mm template.
    """
    id_ = 'pain_01.nidm-1'
    kern = kernel.MKDAKernel(r=4, value=1)
    ma_maps = kern.transform(testdata['dset'].coordinates, testdata['dset'].masker)

    ijk = testdata['dset'].coordinates.loc[testdata['dset'].coordinates['id'] == id_,
                                           ['i', 'j', 'k']]
    ijk = np.squeeze(ijk.values.astype(int))
    kern_data = ma_maps[0].get_fdata()
    com = np.array(center_of_mass(kern_data)).astype(int).T
    com = np.squeeze(com)
    assert np.array_equal(ijk, com) 

def test_mkdakernel1(testdata):
    """
    COMs of MKDA kernel maps should match the foci fed in (assuming focus isn't
    masked out and spheres don't overlap).
    Test on 1mm template.
    """
    id_ = 'pain_01.nidm-1'
    kern = kernel.MKDAKernel(r=4, value=1)
    ma_maps = kern.transform(testdata['dset'].coordinates, testdata['dset'].masker)

    ijk = testdata['dset'].coordinates.loc[testdata['dset'].coordinates['id'] == id_,
                                           ['i', 'j', 'k']]
    ijk = np.squeeze(ijk.values.astype(int))
    kern_data = ma_maps[0].get_fdata()
    com = np.array(center_of_mass(kern_data)).astype(int).T
    com = np.squeeze(com)
    assert np.array_equal(ijk, com) 

def test_integrity(param_test):
    """
    Test integrity of function
    """

    # open ground truth
    im_seg_manual = Image(param_test.fname_gt).change_orientation("RPI")

    # Compute center of mass of the SC seg on each axial slice.
    center_of_mass_x_y_z_lst = [[int(center_of_mass(im_seg_manual.data[:, :, zz])[0]),
                                 int(center_of_mass(im_seg_manual.data[:, :, zz])[1]),
                                 zz] for zz in range(im_seg_manual.dim[2])]

    im_ctr_manual = msct_image.zeros_like(im_seg_manual)
    for x_y_z in center_of_mass_x_y_z_lst:
        im_ctr_manual.data[x_y_z[0], x_y_z[1], x_y_z[2]] = 1

    # open output segmentation
    im_ctr = Image(param_test.file_ctr).change_orientation("RPI")

    # compute MSE between generated ctr and ctr from database
    mse_detection = compute_mse(im_ctr, im_ctr_manual)

    param_test.output += 'Computed MSE: ' + str(mse_detection)
    param_test.output += 'MSE threshold (if computed MSE higher: fail): ' + str(param_test.mse_threshold)

    if mse_detection > param_test.mse_threshold:
        param_test.status = 99
        param_test.output += '--> FAILED'
    else:
        param_test.output += '--> PASSED'

    # update Panda structure
    param_test.results['mse_detection'] = mse_detection

    return param_test 

def find_local_maxima(scmap, radius, threshold):
    grid = peak_local_max(
        scmap,
        min_distance=radius,
        threshold_abs=threshold,
        exclude_border=False,
        indices=False,
    )
    labels = measurements.label(grid)[0]
    xy = measurements.center_of_mass(grid, labels, range(1, np.max(labels) + 1))
    return np.asarray(xy, dtype=np.int).reshape((-1, 2)) 

def extract_sagital_slice(self):
        """Extract the sagital slice where the detection is done.

        If the segmentation is provided,
            the 2D sagital slice is choosen accoding to the segmentation.

        If the segmentation is not provided,
            the 2D sagital slice is choosen as the mid-sagital slice of the input image.
        """
        # TODO: get the mean across multiple sagittal slices to reduce noise

        if self.fname_seg is not None:
            img_seg = Image(self.fname_seg)

            z_mid_slice = img_seg.data[:, int(img_seg.dim[1] / 2), :]
            if 1 in z_mid_slice:  # if SC segmentation available at this slice
                self.rl_coord = int(center_of_mass(z_mid_slice)[1])  # Right_left coordinate
            else:
                self.rl_coord = int(img_seg.dim[2] / 2)
            del img_seg

        else:
            img = Image(self.fname_im)
            self.rl_coord = int(img.dim[2] / 2)  # Right_left coordinate
            del img

        sct.run(['sct_crop_image', '-i', self.fname_im, '-zmin', str(self.rl_coord), '-zmax', str(self.rl_coord + 1), '-o', self.slice2D_im]) 

def scan_slice(z_slice, model, mean_train, std_train, coord_lst, patch_shape, z_out_dim):
    """Scan the entire axial slice to detect the centerline."""
    z_slice_out = np.zeros(z_out_dim)
    sum_lst = []
    # loop across all the non-overlapping blocks of a cross-sectional slice
    for idx, coord in enumerate(coord_lst):
        block = z_slice[coord[0]:coord[2], coord[1]:coord[3]]
        block_nn = np.expand_dims(np.expand_dims(block, 0), -1)
        block_nn_norm = _normalize_data(block_nn, mean_train, std_train)
        block_pred = model.predict(block_nn_norm, batch_size=BATCH_SIZE)

        if coord[2] > z_out_dim[0]:
            x_end = patch_shape[0] - (coord[2] - z_out_dim[0])
        else:
            x_end = patch_shape[0]
        if coord[3] > z_out_dim[1]:
            y_end = patch_shape[1] - (coord[3] - z_out_dim[1])
        else:
            y_end = patch_shape[1]

        z_slice_out[coord[0]:coord[2], coord[1]:coord[3]] = block_pred[0, :x_end, :y_end, 0]
        sum_lst.append(np.sum(block_pred[0, :x_end, :y_end, 0]))

    # Put first the coord of the patch were the centerline is likely located so that the search could be faster for the
    # next axial slices
    coord_lst.insert(0, coord_lst.pop(sum_lst.index(max(sum_lst))))

    # computation of the new center of mass
    if np.max(z_slice_out) > 0.5:
        z_slice_out_bin = z_slice_out > 0.5
        labeled_mask, numpatches = label(z_slice_out_bin)
        largest_cc_mask = (labeled_mask == (np.bincount(labeled_mask.flat)[1:].argmax() + 1))
        x_CoM, y_CoM = center_of_mass(largest_cc_mask)
        x_CoM, y_CoM = int(x_CoM), int(y_CoM)
    else:
        x_CoM, y_CoM = None, None

    return z_slice_out, x_CoM, y_CoM, coord_lst 

def crop_image_around_centerline(im_in, ctr_in, crop_size):
    """Crop the input image around the input centerline file."""
    data_ctr = ctr_in.data
    data_ctr = data_ctr if len(data_ctr.shape) >= 3 else np.expand_dims(data_ctr, 2)
    data_in = im_in.data.astype(np.float32)
    im_new = empty_like(im_in)  # but in fact we're going to crop it

    x_lst, y_lst, z_lst = [], [], []
    data_im_new = np.zeros((crop_size, crop_size, im_in.dim[2]))
    for zz in range(im_in.dim[2]):
        if np.any(np.array(data_ctr[:, :, zz])):
            x_ctr, y_ctr = center_of_mass(np.array(data_ctr[:, :, zz]))

            x_start, x_end = _find_crop_start_end(x_ctr, crop_size, im_in.dim[0])
            y_start, y_end = _find_crop_start_end(y_ctr, crop_size, im_in.dim[1])

            crop_im = np.zeros((crop_size, crop_size))
            x_shape, y_shape = data_in[x_start:x_end, y_start:y_end, zz].shape
            crop_im[:x_shape, :y_shape] = data_in[x_start:x_end, y_start:y_end, zz]

            data_im_new[:, :, zz] = crop_im

            x_lst.append(str(x_start))
            y_lst.append(str(y_start))
            z_lst.append(zz)

    im_new.data = data_im_new
    return x_lst, y_lst, z_lst, im_new 

def centroids(A, window=None):
    A = A.load().dropna('unit_id', how='all')
    if not A.size > 0:
        return pd.DataFrame()
    if window is None:
        window = A.isnull().sum('unit_id') == 0
    try:
        A = A.where(window, drop=True)
    except:
        set_trace()
    A = A.fillna(0)
    meta_dims = set(A.coords.keys()) - set(A.dims)
    meta_dict = {dim: A.coords[dim].values for dim in meta_dims}
    cur_meta = pd.Series(meta_dict)
    cts_list = []
    for uid, cur_uA in A.groupby('unit_id'):
        cur_A = cur_uA.values
        if not (cur_A > 0).any():
            continue
        cur_idxs = cur_uA.dims
        cur_cts = center_of_mass(cur_A)
        cur_cts = pd.Series(cur_cts, index=cur_idxs)
        cur_cts = cur_cts.append(pd.Series(dict(unit_id=uid)))
        cur_cts = cur_cts.append(cur_meta)
        cts_list.append(cur_cts)
    try:
        cts_df = pd.concat(cts_list, axis=1, ignore_index=True).T
    except ValueError:
        cts_df = pd.DataFrame()
    return cts_df 

def calculate_centroids(a):
    print("calculating centroids for " + a.name)
    centroids = np.zeros((a.shape[2], 2))
    for idu, u in enumerate(centroids):
        centroids[idu, :] = center_of_mass(a.values[:, :, idu])
    centroids = xr.DataArray(
        centroids.T,
        coords={'centloc': ['cy', 'cx'],
                'unitid': range(a.shape[2])},
        dims=('centloc', 'unitid'),
        name=a.name)
    return centroids 

def centroid(A, verbose=False):
    def rel_cent(im):
        im_nan = np.isnan(im)
        if im_nan.all():
            return np.array([np.nan, np.nan])
        if im_nan.any():
            im = np.nan_to_num(im)
        cent = np.array(center_of_mass(im))
        return cent / im.shape
    gu_rel_cent = da.gufunc(
        rel_cent,
        signature='(h,w)->(d)',
        output_dtypes=float,
        output_sizes=dict(d=2),
        vectorize=True
    )
    cents = (xr.apply_ufunc(
        gu_rel_cent, A.chunk(dict(height=-1, width=-1)),
        input_core_dims=[['height', 'width']],
        output_core_dims=[['dim']],
        dask='allowed')
             .assign_coords(dim=['height', 'width']))
    if verbose:
        print("computing centroids")
        with ProgressBar():
            cents=cents.compute()
    cents_df = (cents.rename('cents').to_series().dropna()
                .unstack('dim').rename_axis(None, axis='columns')
                .reset_index())
    h_rg = (A.coords['height'].min().values, A.coords['height'].max().values)
    w_rg = (A.coords['width'].min().values, A.coords['width'].max().values)
    cents_df['height'] = cents_df['height'] * (h_rg[1] - h_rg[0]) + h_rg[0]
    cents_df['width'] = cents_df['width'] * (w_rg[1] - w_rg[0]) + w_rg[0]
    return cents_df 

def _calculate_contours_centroids(self):
        cnts_df_list = []
        cts_df_list = []
        A = self.cnmf['A'].load()
        for uid in range(self._u):
            cur_A = A.sel(unit_id=uid)
            cur_idxs = cur_A.squeeze().dims
            cur_thres = dask.delayed(cur_A.max)()
            cur_thres = dask.delayed(float)(cur_thres * .3)
            cur_cnts = dask.delayed(find_contours)(cur_A, cur_thres)
            cur_cnts = dask.delayed(np.concatenate)(cur_cnts)
            cur_cnts = dask.delayed(pd.DataFrame)(cur_cnts, columns=cur_idxs)
            cur_cnts = cur_cnts.assign(unit_id=uid)
            cur_cts = dask.delayed(center_of_mass)(cur_A.values)
            cur_cts = dask.delayed(pd.Series)(cur_cts, index=cur_idxs)
            cur_cts = cur_cts.append(pd.Series(dict(unit_id=uid)))
            cnts_df_list.append(cur_cnts)
            cts_df_list.append(cur_cts)
        cnts_df_list = dask.compute(*cnts_df_list)
        cts_df_list = dask.compute(*cts_df_list)
        cnts_df = pd.concat(cnts_df_list)
        cts_df = pd.concat(cts_df_list, axis=1).T
        for dim in cur_idxs:
            cnts_df[dim].update(cnts_df[dim] / A.sizes[dim] * self._dims[dim])
            cts_df[dim].update(cts_df[dim] / A.sizes[dim] * self._dims[dim])
        return cnts_df, cts_df 

def calculate_centroids_old(*args):
    ndims = np.array(np.ndim(a) for a in args)
    if np.any(ndims < 3):
        raise AssertionError("not a spatial matrix. reshape first!")
    nunits = tuple(a.shape[-1] for a in args)
    centroids = list()
    for ida, cur_a in enumerate(args):
        print ("calculating centroids for matrix " + str(ida))
        cur_centroid = np.zeros((nunits[ida], 2))
        for idu, u in enumerate(cur_centroid):
            cur_centroid[idu, :] = center_of_mass(cur_a[:, :, idu])
        centroids.append(cur_centroid)
    return centroids 

def unet_candidates():
    cands = glob.glob("../data/predictions_epoch9_23_all/*.png")
    #df = pd.DataFrame(columns=['seriesuid','coordX','coordY','coordZ','class'])
    data = []
    imname = ""
    origin = []
    spacing = []
    nrimages = 0
    for name in tqdm(cands):

        #image = imread(name)
        image_t = imread(name)
        image_t = image_t.transpose()
        #Thresholding
        image_t[image_t<THRESHOLD] = 0
        image_t[image_t>0] = 1
        #erosion
        selem = morphology.disk(1)
        image_eroded = image_t
        image_eroded = morphology.binary_erosion(image_t,selem=selem)

        label_im, nb_labels = ndimage.label(image_eroded)
        imname3 = os.path.split(name)[1].replace('.png','')

        splitted = imname3.split("slice")
        slice = splitted[1]
        imname2 = splitted[0][:-1]
        centers = []
        for i in xrange(1,nb_labels+1):
            blob_i = np.where(label_im==i,1,0)
            mass = center_of_mass(blob_i)
            centers.append([mass[1],mass[0]])


        if imname2 != imname:
            if os.path.isfile("../data/1_1_1mm_512_x_512_annotation_masks/spacings/{0}.pickle".format(imname2)):
                with open("../data/1_1_1mm_512_x_512_annotation_masks/spacings/{0}.pickle".format(imname2), 'rb') as handle:
                    dic = pickle.load(handle)
                    origin = dic["origin"]
                    spacing = dic["spacing"]

            imname = imname2
            nrimages +=1

        for center in centers:
            coords = voxel_2_world([int(slice),center[1]+(512-324)*0.5,center[0]+(512-324)*0.5],origin,spacing)
            data.append([imname2,coords[2],coords[1],coords[0],'?'])

        #if nrimages == 5:
        #    break

    df = pd.DataFrame(data,columns=CANDIDATES_COLUMNS)
    save_candidates("../data/candidates_unet_final_23.csv",df) 

def _peak_position(img, nwindow, plot=False):
    """
    function to determine the peak position on the detector using the center of
    mass (COM)

    Parameters
    ----------
    img :       array-like
        detector image data as 2D array
    nwindow :   int
        to avoid influence of hot pixels far away from the peak position the
        center of mass approach is repeated with a window around the COM of the
        full image. COM of the size (nwindow, nwindow) is returned
    plot :      bool, optional
        the result of the of the determination can be saved as a plot
    """
    nw = nwindow // 2
    [cen1r, cen2r] = center_of_mass(img)
    for i in range(11):  # refine center of mass multiple times
        [cen1, cen2] = center_of_mass(
            img[max(int(cen1r) - nw, 0):
                min(int(cen1r) + nw, img.shape[0]),
                max(int(cen2r) - nw, 0):
                min(int(cen2r) + nw, img.shape[1])])
        cen1 += max(int(cen1r) - nw, 0)
        cen2 += max(int(cen2r) - nw, 0)
        if numpy.linalg.norm((cen1 - cen1r, cen2 - cen2r)) > 3:
            cen1r, cen2r = (cen1, cen2)
        else:
            break
    if i == 10 and config.VERBOSITY >= config.INFO_LOW:
        print("XU.analysis._peak_position: Warning: peak position "
              "determination not converged, consider debug mode!")
    if plot:
        plot, plt = utilities.import_matplotlib_pyplot('XU.analysis._peak_'
                                                       'position')
    if plot:
        plt.figure("_ccd")
        plt.imshow(utilities.maplog(img), origin='low')
        plt.plot(cen2, cen1, 'ow', mfc='none')
        plt.axis([cen2 - nw, cen2 + nw, cen1 - nw, cen1 + nw])
        plt.colorbar()
        fnr = len(glob.glob('xu_calib_ccd_img*.png'))
        plt.savefig("xu_calib_ccd_img%d.png" % (fnr + 1))
        plt.close("_ccd")

    return cen1, cen2