# Sets up vars needed
    meanSum, meanTotal, count = 0.0, 0.0, 0
    hc, hd, m, sdc, sdd = 4, 4, 20, 3, 3
    hci, hdi = 1.0/hc**2, 1.0/hd**2

    # Sets current pixel and compares against rest
    x = img[i,j,0]
    for k in range(img.shape[0]):
        for l in range(img.shape[1]):
            # Gets next pixel intensity
            xi = img[k,l,0]

            # Calculates spatial and intensity distances
            magHc = abs(x-xi)
            magHd = math.sqrt(math.pow(i-k, 2) + math.pow(j-l, 2))

            # If within distances
            if magHc <= sdc*hc and magHd <= sdd*hd:
                count += 1

                # Calculates
                xxia, xxib, xxic = (x-xi)**2, (i-k)**2, (j-l)**2
                xxi = xxia*hci + xxib*hdi + xxic*hdi
                exp = math.exp(-0.5 * xxi)
                # magxi = math.sqrt(math.pow(xi, 2) + math.pow(k, 2) + math.pow(l, 2))

                # Adds to sums
                meanSum += xi * exp
                meanTotal += exp

    # Clustering
    if m < count:
        newImg[i,j,0] = meanSum / meanTotal