python-siftgpu(OpenCL版本)--并且改为ASIFT-GPU

https://github.com/silx-kit/silx

http://www.silx.org/doc/silx/dev/install.html

http://www.silx.org/doc/silx/

http://www.silx.org/doc/silx/latest/

http://www.silx.org/doc/silx/dev/sample_code/index.html

 

https://github.com/pierrepaleo/sift_pyocl/blob/master/test/demo_match.py

sudo pip install silx 

from silx.image import sift
import cv2
img = cv2.imread("./11_IMG_Texture_8Bit.png", cv2.IMREAD_GRAYSCALE)
siftp = sift.SiftPlan(img.shape, img.dtype, devicetype="GPU")
kp = siftp.keypoints(img)

1原始SIFT-GPU (OpenCL版本)

# Python 2/3 compatibility
from __future__ import print_function
import sys
import numpy as np
import cv2 as cv

# local modules
from common import Timer
from descriptor import init_feature, filter_matches, explore_match

from silx.image import sift



import numpy
import pylab
from matplotlib.patches import ConnectionPatch
import cv2


exp = 0
if exp == 0:
    fn1 = "./save_ply/OtherSampleFrame_IMG_Texture_8Bit_48.png"
    fn2 = "./save_ply/OtherSampleFrame_IMG_Texture_8Bit_52.png"

    depfn1 = "./save_ply/OtherSampleFrame_IMG_DepthMap_48.tif"
    depfn2 = "./save_ply/OtherSampleFrame_IMG_DepthMap_52.tif"
else:
    fn1 = "./save_ply/11_IMG_Texture_8Bit.png"
    fn2 = "./save_ply/22_IMG_Texture_8Bit.png"

    depfn1 = "./save_ply/11_IMG_DepthMap.tif"
    depfn2 = "./save_ply/22_IMG_DepthMap.tif"

img1 = cv.imread(fn1, cv.IMREAD_GRAYSCALE)
img2 = cv.imread(fn2, cv.IMREAD_GRAYSCALE)

depImage1 = cv.imread(depfn1, -1)
depImage2 = cv.imread(depfn2, -1)

if img1 is None:
    print('Failed to load fn1:', fn1)
    sys.exit(1)

if img2 is None:
    print('Failed to load fn2:', fn2)
    sys.exit(1)


siftp = sift.SiftPlan(img1.shape, img1.dtype, devicetype="GPU")
keypoints1 = siftp.keypoints(img1)
keypoints2 = siftp.keypoints(img2)

print("After keypoints1 %i" % keypoints1.shape[0])
print("After keypoints2 %i" % keypoints2.shape[0])

fig = pylab.figure()
sp1 = fig.add_subplot(122)
sp2 = fig.add_subplot(121)
im1 = sp1.imshow(img1)
im2 = sp2.imshow(img2)

siftmatch = sift.MatchPlan(devicetype="GPU")
matching = siftmatch.match(keypoints1, keypoints2)

h1, w1 = img1.shape[:2]
h2, w2 = img2.shape[:2]
vis = np.zeros((max(h1, h2), w1 + w2), np.uint8)
vis[:h1, :w1] = img1
vis[:h2, w1:w1 + w2] = img2
vis = cv2.cvtColor(vis, cv2.COLOR_GRAY2BGR)

green = (0, 255, 0)
for left, right in zip(matching[:, 0], matching[:, 1]):
    cv2.line(vis, (int(left.x), int(left.y)), (int(right.x+w1), int(right.y)), green, thickness=4)
cv2.imwrite("./save_ply/No2_FeaturePointElimination.png", vis)

 

2改为ASIFT-GPU(OpenCL版本)

 

# Python 2/3 compatibility
from __future__ import print_function
import sys
import numpy as np
import cv2 as cv

# local modules
from common import Timer
from descriptor import init_feature, filter_matches, explore_match

from silx.image import sift



import numpy
import pylab
from matplotlib.patches import ConnectionPatch
import cv2

def affine_skew(tilt, phi, img, mask=None):
    '''
    affine_skew(tilt, phi, img, mask=None) -> skew_img, skew_mask, Ai
    Ai - is an affine transform matrix from skew_img to img
    '''
    h, w = img.shape[:2]
    if mask is None:
        mask = np.zeros((h, w), np.uint8)
        mask[:] = 255
    A = np.float32([[1, 0, 0], [0, 1, 0]])
    if phi != 0.0:
        phi = np.deg2rad(phi)
        s, c = np.sin(phi), np.cos(phi)
        A = np.float32([[c,-s], [ s, c]])
        corners = [[0, 0], [w, 0], [w, h], [0, h]]
        tcorners = np.int32( np.dot(corners, A.T) )
        x, y, w, h = cv.boundingRect(tcorners.reshape(1,-1,2))
        A = np.hstack([A, [[-x], [-y]]])
        img = cv.warpAffine(img, A, (w, h), flags=cv.INTER_LINEAR, borderMode=cv.BORDER_REPLICATE)
    if tilt != 1.0:
        s = 0.8*np.sqrt(tilt*tilt-1)
        img = cv.GaussianBlur(img, (0, 0), sigmaX=s, sigmaY=0.01)
        img = cv.resize(img, (0, 0), fx=1.0/tilt, fy=1.0, interpolation=cv.INTER_NEAREST)
        A[0] /= tilt
    if phi != 0.0 or tilt != 1.0:
        h, w = img.shape[:2]
        mask = cv.warpAffine(mask, A, (w, h), flags=cv.INTER_NEAREST)
    Ai = cv.invertAffineTransform(A)
    return img, mask, Ai

def affine_detect(img, mask=None):
    '''
    affine_detect(detector, img, mask=None, pool=None) -> keypoints, descrs
    Apply a set of affine transformations to the image, detect keypoints and
    reproject them into initial image coordinates.
    See http://www.ipol.im/pub/algo/my_affine_sift/ for the details.
    ThreadPool object may be passed to speedup the computation.
    '''

    params = [(1.0, 0.0)]
    for t in 2**(0.5*np.arange(1,6)):
        for phi in np.arange(0, 180, 72.0 / t):
            params.append((t, phi))

    hh,ww = img.shape[:2]

    keypointa_all = []

    total_size = 0

    dtype_kp = numpy.dtype([('x', numpy.float32),
                            ('y', numpy.float32),
                            ('scale', numpy.float32),
                            ('angle', numpy.float32),
                            ('desc', (numpy.uint8, 128))
                            ])

    for i, (k, d) in enumerate(params):
        t, phi = k, d
        timg, tmask, Ai = affine_skew(t, phi, img)
        img_disp = cv.bitwise_and(timg, timg, mask=tmask);

        siftp = sift.SiftPlan(img_disp.shape, img_disp.dtype, devicetype="GPU")
        keypoints = siftp.keypoints(img_disp)

        for kp in keypoints:
            x = kp.x
            y = kp.y
            pt_x,pt_y = tuple(np.dot(Ai, (x, y, 1)))
            # Out of bounds judgment
            if ((pt_x<0) or (pt_y<0) or (pt_x > ww-1) or (pt_y > hh-1)):
                if (pt_x < 0):
                    kp.x, kp.y = (0, pt_y)
                if (pt_y < 0):
                    kp.x, kp.y = (pt_x, 0)
                if (pt_x > ww-1):
                    kp.x, kp.y = (ww-1, pt_y)
                if (pt_y > hh-1):
                    kp.x, kp.y = (pt_x, hh-1)
        keypointa_all.append(keypoints)
        total_size += len(keypoints)

    output = numpy.recarray(shape=(total_size,), dtype=dtype_kp)
    last = 0
    for ds in keypointa_all:
        l = ds.shape[0]
        if l > 0:
            output[last:last + l].x = ds[0:l].x
            output[last:last + l].y = ds[0:l].y
            output[last:last + l].scale = ds[0:l].scale
            output[last:last + l].angle = ds[0:l].angle
            output[last:last + l].desc = ds[0:l].desc
            last += l
    return output


def main():
    feature_name = 'sift'
    exp = 1
    if exp == 0:
        fn1 = "./save_ply/OtherSampleFrame_IMG_Texture_8Bit_48.png"
        fn2 = "./save_ply/OtherSampleFrame_IMG_Texture_8Bit_52.png"

        depfn1 = "./save_ply/OtherSampleFrame_IMG_DepthMap_48.tif"
        depfn2 = "./save_ply/OtherSampleFrame_IMG_DepthMap_52.tif"
    else:
        fn1 = "./save_ply/P1000965.JPG"
        fn2 = "./save_ply/P1000966.JPG"

        depfn1 = "./save_ply/11_IMG_DepthMap.tif"
        depfn2 = "./save_ply/22_IMG_DepthMap.tif"

    img1 = cv.imread(fn1, cv.IMREAD_GRAYSCALE)
    img2 = cv.imread(fn2, cv.IMREAD_GRAYSCALE)

    depImage1 = cv.imread(depfn1, -1)
    depImage2 = cv.imread(depfn2, -1)

    if img1 is None:
        print('Failed to load fn1:', fn1)
        sys.exit(1)

    if img2 is None:
        print('Failed to load fn2:', fn2)
        sys.exit(1)

    with Timer('kp1'):
        keypoints1 = affine_detect(img1)
    with Timer('kp2'):
        keypoints2 = affine_detect(img2)

    print("keypoints1 %i" % len(keypoints1))
    print("keypoints2 %i" % len(keypoints2))

    siftmatch = sift.MatchPlan(devicetype="GPU")
    matching = siftmatch.match(keypoints1, keypoints2)

    h1, w1 = img1.shape[:2]
    h2, w2 = img2.shape[:2]
    vis = np.zeros((max(h1, h2), w1 + w2), np.uint8)
    vis[:h1, :w1] = img1
    vis[:h2, w1:w1 + w2] = img2
    vis = cv2.cvtColor(vis, cv2.COLOR_GRAY2BGR)

    green = (0, 255, 0)
    for left, right in zip(matching[:, 0], matching[:, 1]):
        cv2.line(vis, (int(left.x), int(left.y)), (int(right.x+w1), int(right.y)), green, thickness=4)
    cv2.imwrite("./save_ply/No2_FeaturePointElimination.png", vis)

if __name__ == '__main__':
    with Timer('all time:'):
        main()