sift without normalization.

display.lua

@@ -65,18 +65,22 @@ function fex.imshow(im,params)
     require 'qtwidget'
     require 'qttorch'
     params = params or {}
-    local xx = fex.imToDisplay(im,params)
+    local title = params.title or "Image Display"
+    local xx = fex.imToDisplay(im:clone(),params)
     local w,h = xx:size(xx:dim()),xx:size(xx:dim()-1)
-    local ww = params.ww or qtwidget.newwindow(w,h,'Image Display')
+    local ww = params.win or qtwidget.newwindow(w,h,title)
     local xi = params.x or 0
     local yi = params.y or 0
     local qim = qt.QImage.fromTensor(xx)
-    local wr,hr = ww:currentsize()
-    local ss = math.min(wr/w,hr/h)
-    local wi,hi = w*ss,h*ss
-    ww:resize(wi,hi)
+    if params.win == nil then
+        local wr,hr = ww:currentsize()
+        local ss = math.min(wr/w,hr/h)
+        local wi,hi = w*ss,h*ss
+        ww:resize(wi,hi)
+        ww:onResize(function(w,h) ww:image(xi,yi,w,h,qim) end)
+    end
+    ww.widget.windowTitle = title
     ww:image(xi,yi,wi,hi,qim)
-    ww:onResize(function(w,h) ww:image(xi,yi,w,h,qim) end)
     return ww
 end
 

sift.lua

@@ -1,3 +1,33 @@
+
+
+local function csize(x,k,s)
+    local o = math.floor((x-k)/s+1)
+    local xi = (o-1)*s+k
+    return xi,o
+end
+-- patch_size is the region over which the SIFT feature is computed
+-- grid_spacing is the step size between each patch
+-- num_bins is the 
+function sift_sampler(io,patch_size,grid_spacing,num_bins)
+    local num_angles = io:size(1)
+    local binstep = math.floor(patch_size/num_bins)
+    local nr,nor = csize(io:size(2),patch_size,grid_spacing)
+    local nc,noc = csize(io:size(3),patch_size,grid_spacing)
+    local shiftr = 1+math.floor((io:size(2)-nr)/2)
+    local shiftc = 1+math.floor((io:size(3)-nc)/2)
+    local nio = io:narrow(2,shiftr,nr):narrow(3,shiftc,nc)
+    local uio = nio:unfold(2,patch_size,grid_spacing):unfold(3,patch_size,grid_spacing)
+    local sift = torch.Tensor(num_bins*num_bins*num_angles,nor,noc)
+    local b = 1
+    for i=1,num_bins do
+        for j=1,num_bins do
+            sift:narrow(1,b,num_angles):copy(uio:select(5,(i-1)*binstep+1):select(4,(j-1)*binstep+1))
+            b = b + num_angles
+        end
+    end
+    return sift
+end
+
 -- SIFT feature extractor as implemented in Svetlana Lazebnik's
 -- Pyramid code.
 function fex.dsift(im, params)
@@ -6,7 +36,7 @@ function fex.dsift(im, params)
 
     params = params or {}
     local patch_size = params.patch_size or 16
-    local grid_spacing = params.grid_spacing or 1
+    local grid_spacing = params.grid_spacing or 8
 
     local num_angles = 8
     local num_bins = 4
@@ -20,8 +50,10 @@ function fex.dsift(im, params)
 
     local I
     if im:dim() == 3 then
-        I = torch.sum(im,1)/3
-        I = I[1]
+        I = torch.Tensor(im:size(2),im:size(3))
+        torch.add(I,im[1],im[2])
+        I:add(im[3])
+        I:div(3)
     elseif im:dim() == 2 then
         I = im
     else
@@ -41,7 +73,7 @@ function fex.dsift(im, params)
 
     I_X = I_X:narrow(1,3,hgt):narrow(2,3,wid):clone()
     I_Y = I_Y:narrow(1,3,hgt):narrow(2,3,wid):clone()
-    print('1',torch.toc(t))
+    if fex.verbose then print('1',torch.toc(t)) end
 
     local I_theta = torch.atan2(I_Y,I_X)
     I_theta[torch.ne(I_theta,I_theta)]=0
@@ -54,41 +86,144 @@ function fex.dsift(im, params)
     local grid_x = torch.range(patch_size/2,wid-patch_size/2+1,grid_spacing)
     local grid_y = torch.range(patch_size/2,hgt-patch_size/2+1,grid_spacing)
 
-    print('2',torch.toc(t))
+    if fex.verbose then print('2',torch.toc(t)) end
     local I_orientation = torch.Tensor(num_angles, hgt, wid)
 
     local cosI = torch.cos(I_theta)
     local sinI = torch.sin(I_theta)
 
-    print('3',torch.toc(t))
+    if fex.verbose then print('3',torch.toc(t)) end
 
+    local tmp = torch.Tensor(num_angles,hgt,wid)
     for a=1,num_angles do
-        local tmp = cosI*math.cos(angles[a]) + sinI*math.sin(angles[a])
-        tmp:pow(alpha)
-        tmp[torch.le(tmp,0)] = 0
-        torch.cmul(I_orientation[a], tmp, I_mag)
+        torch.mul(tmp[a],cosI,math.cos(angles[a]))
+        tmp[a]:add(math.sin(angles[a]),sinI)
     end
-    print('4',torch.toc(t))
+    tmp:pow(alpha)
+    tmp[torch.le(tmp,0)] = 0
+    local tt = torch.Tensor(I_mag:storage(),I_mag:storageOffset(),num_angles,0,hgt,I_mag:stride(1),wid,I_mag:stride(2))
+    torch.cmul(I_orientation,tmp,tt)
+
+    if fex.verbose then print('4',torch.toc(t)) end
 
     local weight_kernel = torch.zeros(patch_size,patch_size)
     local r = patch_size/2;
     local cx = r-0.5;
     local sample_res = patch_size/num_bins;
     local weight_x = torch.abs(torch.range(1,patch_size)-cx)/sample_res
     weight_x:apply(function(x) if x <= 1 then return 1-x else return 0 end end)
-    local tw = torch.ger(weight_x, weight_x)
-    print('5',torch.toc(t))
+
+    if fex.verbose then print('5',torch.toc(t)) end
+
     local wx= torch.Tensor(weight_x):resize(weight_x:size(1),1)
 
-    --local I_orientation2 = torch.Tensor(num_angles, hgt-wx:size(1)+1, wid-wx:size(1)+1)
     local I_orientation2 = torch.Tensor():resizeAs(I_orientation)
     for a=1,num_angles do
-        local t = fex.conv2(I_orientation[a],wx,'S')
-        fex.conv2(I_orientation2[a],t,wx:t(),'S')
-        --fex.conv2(I_orientation2[a],I_orientation[a],tw,'S')
+        local t = fex.conv2(I_orientation[a],wx:t(),'S')
+        fex.conv2(I_orientation2[a],t,wx,'S')
+    end
+
+    if fex.verbose then print('6',torch.toc(t)) end
+
+    local sift=sift_sampler(I_orientation2,patch_size,grid_spacing,num_bins)
+
+    collectgarbage()
+    if fex.verbose then print('7',torch.toc(t)) end
+    return sift
+end
+
+-- SIFT feature extractor as implemented in Svetlana Lazebnik's
+-- Pyramid code.
+function fex.dsiftfast(im, params)
+
+    local t = torch.tic()
+
+    params = params or {}
+    local patch_size = params.patch_size or 16
+    local grid_spacing = params.grid_spacing or 8
+
+    local num_angles = 8
+    local num_bins = 4
+    local num_samples = num_bins * num_bins
+    local alpha = 9
+    local sigma_edge = 1
+
+    local angle_step = 2 * math.pi / num_angles
+    local angles = torch.range(0,2*math.pi,angle_step)
+    angles = angles:narrow(1,1,num_angles)
+
+    local I
+    if im:dim() == 3 then
+        I = torch.Tensor(im:size(2),im:size(3))
+        torch.add(I,im[1],im[2])
+        I:add(im[3])
+        I:div(3)
+    elseif im:dim() == 2 then
+        I = im
+    else
+        error('im has to be 2D or 3D')
+    end
+
+    I:div(I:max())
+
+    local hgt,wid = I:size(1),I:size(2)
+    local G_X,G_Y = fex.gendgauss(sigma_edge)
+
+    I:add(-I:mean())
+    print('0',torch.toc(t))
+    local I_X = torch.xcorr2(I, G_X)
+    local I_Y = torch.xcorr2(I, G_Y)
+
+    if fex.verbose then print('1',torch.toc(t)) end
+
+    local I_theta = torch.atan2(I_Y,I_X)
+    I_theta[torch.ne(I_theta,I_theta)]=0
+    I_X:cmul(I_X)
+    I_Y:cmul(I_Y)
+    I_X:add(I_Y)
+    local I_mag = I_X:sqrt()
+    hgt,wid = I_mag:size(1),I_mag:size(2)
+
+    if fex.verbose then print('2',torch.toc(t)) end
+    local I_orientation = torch.Tensor(num_angles, hgt, wid)
+    local cosI = torch.cos(I_theta)
+    local sinI = torch.sin(I_theta)
+
+    if fex.verbose then print('3',torch.toc(t)) end
+
+    local tmp = torch.Tensor(num_angles,hgt,wid)
+    for a=1,num_angles do
+        torch.mul(tmp[a],cosI,math.cos(angles[a]))
+        tmp[a]:add(math.sin(angles[a]),sinI)
+    end
+    tmp:pow(alpha)
+    tmp[torch.le(tmp,0)] = 0
+    local tt = torch.Tensor(I_mag:storage(),I_mag:storageOffset(),num_angles,0,hgt,I_mag:stride(1),wid,I_mag:stride(2))
+    torch.cmul(I_orientation,tmp,tt)
+    if fex.verbose then print('4',torch.toc(t)) end
+
+    local weight_kernel = torch.zeros(patch_size,patch_size)
+    local r = patch_size/2;
+    local cx = r-0.5;
+    local sample_res = patch_size/num_bins;
+    local weight_x = torch.abs(torch.range(1,patch_size)-cx)/sample_res
+    weight_x:apply(function(x) if x <= 1 then return 1-x else return 0 end end)
+    weight_x = torch.Tensor(weight_x):resize(weight_x:size(1),1)
+
+    if fex.verbose then print('5',torch.toc(t)) end
+
+    local I_orientation2 = torch.Tensor(num_angles, hgt-weight_x:size(1)+1, wid-weight_x:size(1)+1)
+    --local I_orientation2 = torch.Tensor():resizeAs(I_orientation)
+    local tim = torch.Tensor(num_angles, I_orientation2:size(2), wid)
+    for a=1,num_angles do
+        torch.conv2(tim[a],I_orientation[a],weight_x)
+        torch.conv2(I_orientation2[a],tim[a],weight_x:t())
     end
-    print('6',torch.toc(t))
-    return I_orientation2
+    if fex.verbose then print('6',torch.toc(t)) end
+    local sift=sift_sampler(I_orientation2,patch_size,grid_spacing,num_bins)
+    collectgarbage()
+    if fex.verbose then print('7',torch.toc(t)) end
+    return sift
 end
 
 

utils.lua

@@ -14,8 +14,8 @@ function fex.conv2(...)
             ro = x.new(x:size())
         end
         local r = torch.conv2(x,k,'F')
-        local shifti = math.max(1,math.floor((r:size(1)-x:size(1))/2))
-        local shiftj = math.max(1,math.floor((r:size(2)-x:size(2))/2))
+        local shifti = 1+math.ceil((r:size(1)-x:size(1))/2)
+        local shiftj = 1+math.ceil((r:size(2)-x:size(2))/2)
         local ii = r:dim()-1
         local jj = r:dim()
         ro:resizeAs(x)
@@ -41,8 +41,8 @@ function fex.xcorr2(...)
             ro = x.new(x:size())
         end
         local r = torch.xcorr2(x,k,'F')
-        local shifti = math.floor((r:size(1)-x:size(1))/2)
-        local shiftj = math.floor((r:size(2)-x:size(2))/2)
+        local shifti = 1+math.ceil((r:size(1)-x:size(1))/2)
+        local shiftj = 1+math.ceil((r:size(2)-x:size(2))/2)
         local ii = r:dim()-1
         local jj = r:dim()
         ro:resizeAs(x)