first

.gitignore

@@ -55,3 +55,9 @@ docs/_build/
 
 # PyBuilder
 target/
+
+*.bmp
+*.jpg
+*.dm3
+*.JPG
+*.un~

dst/readme.txt

@@ -0,0 +1 @@
+destination folder

fourier.py

@@ -0,0 +1,34 @@
+#!/bin/python2
+import numpy as np
+import matplotlib.pyplot as plt
+
+#signal = [1,2,3,2,4,2,7,2,30,4]*2
+signal = np.random.random_sample(6)
+
+fourier = np.fft.fft(signal)
+freq = np.fft.fftfreq(len(signal))
+
+print fourier, freq
+
+
+x = np.arange(0, len(signal), 0.01);
+ye = [0] * len(x)
+
+for i, c in enumerate(fourier):
+    amp = np.abs(c)/len(signal)
+    pha = np.arctan2(np.imag(c), np.real(c))
+    #print amp, pha
+    if amp != 0:
+        #y = amp*np.cos((x + pha)*2*np.pi/freq[i])
+        y = np.real(c)*np.cos(2*np.pi*x*freq[i]) - np.imag(c)*np.sin(2*np.pi*x*freq[i])
+        y = y/len(signal)
+        #y = np.cos(x*freq[i])
+        #y = np.real(c*np.exp(1j*2*np.pi*x*freq[i]))/len(signal)
+        plt.plot(x, y, ':')
+        ye = np.add(ye, y)
+
+plt.plot(x, ye)
+
+plt.plot(range(len(signal)),signal, 'o')
+
+plt.show()

mine.py

@@ -0,0 +1,223 @@
+#!/bin/python2
+import numpy as np
+import matplotlib.image as MImage
+import matplotlib.pyplot as plt
+import scipy.stats as st
+from PIL import Image, ImageFilter, ImageDraw
+
+import sys
+import os
+
+# Flush STDOUT continuously
+sys.stdout = os.fdopen(sys.stdout.fileno(), 'w', 0)
+
+orgFileName = "11_NP3_011_0.bmp"
+size = (1024, 1024)  # width and height of image in pixels
+#size = (300, 300)  # width and height of image in pixels
+
+######### SOME AUX FUNCTIONS #########
+def gkern(kernlen=21, nsig=3):
+    """Returns a 2D Gaussian kernel array."""
+
+    interval = (2*nsig+1.)/(kernlen)
+    x = np.linspace(-nsig-interval/2., nsig+interval/2., kernlen+1)
+    kern1d = np.diff(st.norm.cdf(x))
+    kernel_raw = np.sqrt(np.outer(kern1d, kern1d))
+    kernel = kernel_raw/kernel_raw.sum()
+    return kernel
+
+def saveImage(name, matrix = None, format = "BMP", image = None):
+    print " -> Saving Image " + name + " ...",
+    if image is not None:
+        pass
+    elif matrix is not None:
+        image =  Image.fromarray(matrix)
+    else:
+        raise ValueError("matrix or image must be specified")
+
+    image.save("dst/" + name + "." + format.lower(), format)
+    print "DONE"
+    return image
+
+
+#######################################
+
+print "Loading and Converting Image ...",
+orgImg = Image.open("src/" + orgFileName, mode = 'r').convert()
+
+try:
+    orgImg = orgImg.convert().crop((0,0,size[0],size[1]))
+except NameError:
+    size = orgImg.size
+
+
+orgGreyscaleMatrix = MImage.pil_to_array(orgImg.convert('L'))
+
+emptyGreyscaleMatrix = MImage.pil_to_array(Image.new('L', size, 0))
+emptyColoredMatrix = MImage.pil_to_array(Image.new('RGB', size, (0,0,0)))
+print "DONE"
+
+
+# print "Compute the 2-dimensional FFT ...",
+# #print "."
+# m = np.fft.rfft2(orgGreyscaleMatrix)
+# #print m.shape
+# #print m
+# #print np.amax(m)
+# am = np.log(np.absolute(m)+1)
+# mm = np.around(am/np.amax(am)*255, decimals=0).astype(np.uint8)
+# print "DONE"
+# #print orgGreyscaleMatrix
+# #print orgGreyscaleMatrix[1, 1]
+# #print mm
+# #print mm[1, 1]
+# saveImage("fft", mm)
+
+
+
+print "Gaussian Filter ... ",
+gaussianMatrix = emptyGreyscaleMatrix*0
+kernelSize = 7
+matrix = gkern(kernelSize)
+for i in range(size[1] - (kernelSize - 1)):
+    for j in range(size[0] - (kernelSize - 1)):
+        iS = i + (kernelSize - 1)/2
+        jS = j + (kernelSize - 1)/2
+        #print t[i:(i + kernelSize), j:(j + kernelSize)]
+        gaussianMatrix[iS, jS] = np.sum(np.multiply(orgGreyscaleMatrix[i:(i + kernelSize), j:(j + kernelSize)], matrix))
+        #s[iS, jS] = np.sum(t[i:(i + kernelSize), j:(j + kernelSize)],)/kernelSize**2
+
+print "DONE"
+
+
+
+# SAVING GAUSSIAN
+saveImage("gaussian", gaussianMatrix)
+
+
+print "Kernel Max Filter and Create Concordance ...",
+maxMatrix = emptyGreyscaleMatrix * 0
+concordanceMatrix = emptyGreyscaleMatrix * 0
+kernelSize = 9
+for i in range(size[1] - (kernelSize - 1)):
+    for j in range(size[0] - (kernelSize - 1)):
+        iS = i + (kernelSize - 1)/2
+        jS = j + (kernelSize - 1)/2
+        #print t[i:(i + kernelSize), j:(j + kernelSize)]
+        maxMatrix[iS, jS] = np.amax(gaussianMatrix[i:(i + kernelSize), j:(j + kernelSize)])
+        if maxMatrix[iS, jS] == gaussianMatrix[iS, jS]:
+            concordanceMatrix[iS, jS] = 1
+        #print iS, jS;
+        #exit();
+print "DONE"
+
+
+
+saveImage("max", maxMatrix)
+concordanceImg = saveImage("concordance", concordanceMatrix*255)
+
+
+print "Color Concordance Image ... ",
+coloredConcordanceMatrix = emptyColoredMatrix * 0
+for i in range(size[1]):
+    for j in range(size[0]):
+        if concordanceMatrix[i, j] == 1:
+            coloredConcordanceMatrix[i, j] = (255, 0, 0)
+print "DONE"
+
+
+
+
+coloredConcordanceImg = Image.fromarray(coloredConcordanceMatrix, 'RGB')
+
+compositeImg = Image.composite(coloredConcordanceImg, orgImg.convert("RGB"), concordanceImg)
+saveImage("composite", image = compositeImg)
+
+#exit()
+
+print "Finding Centers ...",
+kernelSize = 35 # ODD
+maxCenterSize = 2
+movingSize = 20
+
+centerMatrix = emptyGreyscaleMatrix * 0
+
+compImgDraw = ImageDraw.Draw(compositeImg)
+
+
+arcs = []
+dir1 = 73
+dir2 = 145
+tolerance = 20
+maxLength = 18
+
+for i in range(0, size[1] - (kernelSize - 1), movingSize):
+    for j in range(0, size[0] - (kernelSize - 1), movingSize):
+        iS = i + (kernelSize - 1)/2
+        jS = j + (kernelSize - 1)/2
+        areaCenterDict = {}
+        resetList = []
+        #print "Window centered at " + str((iS, jS)) + " size: " + str((kernelSize, kernelSize)) + " borders: " + str((iS - (kernelSize - 1)/2, jS - (kernelSize - 1)/2, iS + (kernelSize - 1)/2, jS + (kernelSize - 1)/2)) +" :"
+        for ii in range(0, kernelSize):
+            for jj in range(0, kernelSize):
+                iiS = iS - (kernelSize - 1)/2 + ii
+                jjS = jS - (kernelSize - 1)/2 + jj
+                if concordanceMatrix[iiS , jjS] == 1:
+                    concordanceMatrix[iiS, jjS] = 0
+                    resetList.append((iiS, jjS))
+                    centerInfo = { "pixels": []}
+                    centerInfo["pixels"].append((iiS, jjS))
+
+                    leftLimit = 0
+                    if ii == 0:
+                        leftLimit = -1
+
+                    topLimit = 0
+                    if jj == 0:
+                        topLimit = -1
+
+                    for ci in range(leftLimit*maxCenterSize, maxCenterSize + 1):
+                        for cj in range(topLimit*maxCenterSize, maxCenterSize + 1):
+                            if ci == 0 and cj == 0:
+                                continue
+                            if concordanceMatrix[iiS + ci, jjS + cj] == 1:
+                                concordanceMatrix[iiS + ci, jjS + cj] = 0
+                                resetList.append((iiS + ci, jjS + cj))
+                                centerInfo["pixels"].append((iiS + ci, jjS + cj))
+                    areaCenterDict[(iiS, jjS)] = centerInfo
+
+        for x, y in resetList:
+            #print "reseting " + str((x, y))
+            concordanceMatrix[x, y] = 1
+
+        #print len(areaCenterDict)
+
+        for center1 in areaCenterDict:
+            for center2 in areaCenterDict:
+                if center1 == center2:
+                    continue
+                dist1 = center1[1] - center2[1]
+                dist2 = center1[0] - center2[0]
+                vLength = np.sqrt(dist1**2 + dist2**2)
+                # Maximum Distance To Connect
+                if vLength < maxLength:
+                    arc = np.arccos(dist1/vLength)/(2*np.pi)*360
+                    if dist2 > 0:
+                        arc = 360 - arc
+                    arcs.append(arc)
+                    if abs(arc - dir1) < tolerance or abs(arc - dir2) < tolerance:
+                        compImgDraw.line([center1[1], center1[0], center2[1], center2[0]], "green")
+
+
+        #centerMatrix[iS, jS] = 255
+print "DONE"
+#a, b = np.histogram(arcs, 100)
+#for i, e in enumerate(a):
+#    print a[i], ":", b[i]
+#plt.hist(arcs, 100)
+#plt.show()
+
+#saveImage("withLines", image = compositeImg)
+
+compositeImg = Image.composite(coloredConcordanceImg, compositeImg, concordanceImg)
+saveImage("withLines", image = compositeImg)

peak.py

@@ -0,0 +1,119 @@
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.ndimage
+import matplotlib.patches
+from PIL import Image
+
+
+plt.figure(figsize=(10,10))
+ax1 = plt.subplot(221)
+ax2 = plt.subplot(222)
+ax3 = plt.subplot(223)
+ax4 = plt.subplot(224)
+
+size        = 1024 # width and height of image in pixels
+peak_height = 100  # define the height of the peaks
+num_peaks   = 20
+noise_level = 50
+threshold   = 60
+
+np.random.seed(3)
+
+#set up a simple, blank image (Z)
+x = np.linspace(0,size,size)
+y = np.linspace(0,size,size)
+
+X,Y = np.meshgrid(x,y)
+Z = X*0
+
+
+tmp = Image.open("base.jpg")
+converted = tmp.convert()
+
+
+for i in range(1024):
+    for j in range(1024):
+        Z[i,j] = converted.getpixel((i, j))[0]
+
+#exit()
+#now add some peaks
+# def gaussian(X,Y,xo,yo,amp=100,sigmax=4,sigmay=4):
+#     return amp*np.exp(-(X-xo)**2/(2*sigmax**2) - (Y-yo)**2/(2*sigmay**2))
+
+# for xo,yo in size*np.random.rand(num_peaks,2):
+#     widthx = 5 + np.random.randn(1)
+#     widthy = 5 + np.random.randn(1)
+#     Z += gaussian(X,Y,xo,yo,amp=peak_height,sigmax=widthx,sigmay=widthy)
+
+#of course, add some noise:
+# Z = Z + scipy.ndimage.gaussian_filter(0.5*noise_level*np.random.rand(size,size),sigma=5)    
+# Z = Z + scipy.ndimage.gaussian_filter(0.5*noise_level*np.random.rand(size,size),sigma=1)    
+
+t = time.time() #Start timing the peak-finding algorithm
+
+#Set everything below the threshold to zero:
+Z_thresh = np.copy(Z)
+Z_thresh[Z_thresh<threshold] = 0
+print 'Time after thresholding: %.5f seconds'%(time.time()-t)
+
+#now find the objects
+labeled_image, number_of_objects = scipy.ndimage.label(Z_thresh)
+print 'Time after labeling: %.5f seconds'%(time.time()-t)
+
+peak_slices = scipy.ndimage.find_objects(labeled_image)
+print 'Time after finding objects: %.5f seconds'%(time.time()-t)
+
+def centroid(data):
+    h,w = np.shape(data)   
+    x = np.arange(0,w)
+    y = np.arange(0,h)
+
+    X,Y = np.meshgrid(x,y)
+
+    cx = np.sum(X*data)/np.sum(data)
+    cy = np.sum(Y*data)/np.sum(data)
+
+    return cx,cy
+
+centroids = []
+
+for peak_slice in peak_slices:
+    dy,dx  = peak_slice
+    x,y = dx.start, dy.start
+    cx,cy = centroid(Z_thresh[peak_slice])
+    centroids.append((x+cx,y+cy))
+
+print 'Total time: %.5f seconds\n'%(time.time()-t)
+
+###########################################
+#Now make the plots:
+for ax in (ax1,ax2,ax3,ax4): ax.clear()
+ax1.set_title('Original image')
+ax1.imshow(Z,origin='lower')
+
+ax2.set_title('Thresholded image')
+ax2.imshow(Z_thresh,origin='lower')
+
+ax3.set_title('Labeled image')
+ax3.imshow(labeled_image,origin='lower') #display the color-coded regions
+
+for peak_slice in peak_slices:  #Draw some rectangles around the objects
+    dy,dx  = peak_slice
+    xy     = (dx.start, dy.start)
+    width  = (dx.stop - dx.start + 1)
+    height = (dy.stop - dy.start + 1)
+    rect = matplotlib.patches.Rectangle(xy,width,height,fc='none',ec='red')
+    ax3.add_patch(rect,)
+
+ax4.set_title('Centroids on original image')
+ax4.imshow(Z,origin='lower')
+
+for x,y in centroids:
+    ax4.plot(x,y,'kx',ms=10)
+
+ax4.set_xlim(0,size)
+ax4.set_ylim(0,size)
+
+plt.tight_layout
+plt.show()

src/readme.txt

@@ -0,0 +1 @@
+source folder