sift matlab

1. f unction [ SIFTFeatureVector, locationX, locationY ] = DenseSIFT( image, nPatc

hSize, nGridSpacing )

2.

3. i mage = double( image );

4. i mage = mean( image, 3 );

5. i mage = image / max( image( : ) );

6.

7. % parameters

8. n AngleNums = 8;

9. n BinNums = 4;

10. nSampleNums = nBinNums * nBinNums;

11. alpha = 9; %% parameter for attenuation of angles (must be odd)

12.

13. if nargin < 5

14. sigmaGuassian = 1;

15. end

16.

17. angleStep = 2 * pi / nAngleNums;

18. angles = 0 : angleStep : 2 * pi;

19. angles( nAngleNums + 1 ) = [ ]; % bin centers

20.

21. [ nRow nCol ] = size( image );

22.

23. [ gaussianX, gaussianY ] = genDeltaGaussian( sigmaGuassian );

24. imageVerticalEdges = filter2( gaussianX, image, 'same' ); % vertical edges

25. imageHorizontalEdges = filter2( gaussianY, image, 'same' ); % horizontal edg

es

26. imageGradientMagnitude = sqrt( imageVerticalEdges.^2 + imageHorizontalEd

ges.^2 ); % gradient magnitude

27. imageTheta = atan2( imageHorizontalEdges, imageVerticalEdges );

28. imageTheta( isnan( imageTheta ) ) = 0; % replace illegal result with 0

29.

30. % descriptor locations

31. locationX = nPatchSize / 2 : nGridSpacing : nCol - nPatchSize / 2 + 1;

32. locationY = nPatchSize / 2 : nGridSpacing : nRow - nPatchSize / 2 + 1;

33.

34. % make orientation images

35. imageOrientation = zeros( [ nRow, nCol, nAngleNums ], 'single' );

36.

37. % for each histogram angle

38. imageCos = cos( imageTheta );

39. imageSin = sin( imageTheta );

40.

41. for index = 1 : nAngleNums

42. % compute each orientation channel

43. tmp = ( imageCos * cos( angles( index ) ) + imageSin * sin( angles( index ) )

).^ alpha;

44. tmp = tmp .* ( tmp > 0 );

45.

46. % weight by magnitude

47. imageOrientation( :, :, index ) = tmp .* imageGradientMagnitude;

48. end

49.

50. % Convolution formulation:

51. nHalfPatchSize = nPatchSize / 2;

52. nHalfPatchSizeMinusDotFive = nHalfPatchSize - 0.5;

53. sampleResolution = nPatchSize / nBinNums;

54. weightX = abs( ( 1 : nPatchSize ) - nHalfPatchSizeMinusDotFive ) / sampleRe

solution;

55. weightX = ( 1 - weightX ) .* ( weightX <= 1 );

56.

57. for index = 1 : nAngleNums

58. imageOrientation( :, :, index ) = conv2( weightX, weightX', imageOrientation

( :, :, index ), 'same' );

59. end

60.

61. % Sample SIFT bins at valid locations (without boundary artifacts)

62. % find coordinates of sample points (bin centers)

63. [ samplePosX, samplePosY ] = meshgrid( linspace( 1, nPatchSize + 1, nBinN

ums + 1 ) );

64. samplePosX = samplePosX( 1 : nBinNums, 1 : nBinNums ); samplePosX = sa

mplePosX( : ) - nPatchSize / 2;

65. samplePosY = samplePosY( 1 : nBinNums, 1 : nBinNums ); samplePosY = sa

mplePosY( : ) - nPatchSize / 2;

66.

67. SIFTFeatureVector = zeros( [ length( locationY ) length( locationX ) nAngleNu

ms * nSampleNums ] , 'single' );

68.

69. nOffset = 0;

70. for n = 1 : nBinNums * nBinNums

71. SIFTFeatureVector( :, :, nOffset + 1 : nOffset + nAngleNums ) = imageOrien

tation( locationY + samplePosY(n), locationX + samplePosX(n), : );

72. nOffset = nOffset + nAngleNums;

73. end

74.

75. clear imageOrientation

76.

77.

78. % Outputs:

79. [ locationX, locationY ] = meshgrid( locationX, locationY );

80. [ nrows, ncols, cols ] = size( SIFTFeatureVector );

81.

82. % normalize SIFT descriptors

83. SIFTFeatureVector = reshape( SIFTFeatureVector, [nrows * ncols nAngleNum

s * nSampleNums ] );

84. SIFTFeatureVector = SIFTNormalization( SIFTFeatureVector );

85. SIFTFeatureVector = reshape( SIFTFeatureVector, [ nrows ncols nAngleNum

s * nSampleNums] );

86.

87.

88. function [ GX, GY ] = genDeltaGaussian( sigma )

89.

90. % laplacian of size sigma

91.

92. G = genGaussian(sigma);

93. [ GX, GY ] = gradient( G );

94.

95. GX = GX * 2 ./ sum( sum( abs( GX ) ) );

96. GY = GY * 2 ./ sum( sum( abs( GY ) ) );

97.

98. function G = genGaussian( sigma )

99.

100. if all( size( sigma ) == [ 1, 1 ] )

101. % isotropic gaussian

102. filterWindow = 4 * ceil( sigma ) + 1;

103. G = fspecial( 'gaussian', filterWindow, sigma );

104. else

105. % anisotropic gaussian

106. filterWindowX = 2 * ceil( sigma( 1 ) ) + 1;

107. filterWindowY = 2 * ceil( sigma( 2 ) ) + 1;

108. GaussianX = normpdf( -filterWindowX: filterWindowX, 0, sigma( 1 ) ); 109. GaussianY = normpdf( -filterWindowY: filterWindowY, 0, sigma( 2 ) ); 110. G = GaussianY' * GaussianX;

111. end

112.

113. function SIFTFeatureVector = SIFTNormalization( SIFTFeatureVector ) 114. % normalize SIFT descriptors (after Lowe)

115.

116. % find indices of descriptors to be normalized (those whose norm is larger t han 1)

117. tmp = sqrt( sum( SIFTFeatureVector.^2, 2 ) );

118. normalizeIndex = find( tmp > 1 );

119.

120. SiftFeatureVectorNormed = SIFTFeatureVector( normalizeIndex, : ); 121. SiftFeatureVectorNormed = SiftFeatureVectorNormed ./ repmat( tmp( norm alizeIndex, : ), [ 1 size( SIFTFeatureVector, 2 ) ] );

122.

123. % suppress large gradients

124. SiftFeatureVectorNormed( SiftFeatureVectorNormed > 0.2 ) = 0.2;

125.

126. % finally, renormalize to unit length

127. tmp = sqrt( sum( SiftFeatureVectorNormed.^2, 2 ) );

128. SiftFeatureVectorNormed = SiftFeatureVectorNormed ./ repmat( tmp, [ 1 siz e( SIFTFeatureVector, 2 ) ] );

129.

130. SIFTFeatureVector( normalizeIndex, : ) = SiftFeatureVectorNormed;