精华内容
下载资源
问答
  • 教你一步一步C语言实现sift算法、上
    原文:http://blog.csdn.net/v_july_v/article/details/6245939 
    

    引言:
        在我写的关于sift算法的前倆篇文章里头,已经对sift算法有了初步的介绍:九、图像特征提取与匹配之SIFT算法,而后在:九(续)、sift算法的编译与实现里,我也简单记录下了如何利用OpenCV,gsl等库编译运行sift程序。
        但据一朋友表示,是否能用c语言实现sift算法,同时,尽量不用到opencv,gsl等第三方库之类的东西。而且,Rob Hess维护的sift 库,也不好懂,有的人根本搞不懂是怎么一回事。
        那么本文,就教你如何利用c语言一步一步实现sift算法,同时,你也就能真正明白sift算法到底是怎么一回事了。

        ok,先看一下,本程序最终运行的效果图,sift 算法分为五个步骤(下文详述),对应以下第二--第六幅图

     

    sift算法的步骤
        要实现一个算法,首先要完全理解这个算法的原理或思想。咱们先来简单了解下,什么叫sift算法:
        sift,尺度不变特征转换,是一种电脑视觉的算法用来侦测与描述影像中的局部性特征,它在空间尺度中寻找极值点,并提取出其位置、尺度、旋转不变量,此算法由 David Lowe 在1999年所发表,2004年完善总结。

        所谓,Sift算法就是用不同尺度(标准差)的高斯函数对图像进行平滑,然后比较平滑后图像的差别,
    差别大的像素就是特征明显的点。

        以下是sift算法的五个步骤:
        一、建立图像尺度空间(或高斯金字塔),并检测极值点

        首先建立尺度空间,要使得图像具有尺度空间不变形,就要建立尺度空间,sift算法采用了高斯函数来建立尺度空间,高斯函数公式为:
           
        上述公式G(x,y,e),即为尺度可变高斯函数。

        而,一个图像的尺度空间L(x,y,e) ,定义为原始图像I(x,y)与上述的一个可变尺度的2维高斯函数G(x,y,e) 卷积运算。
        即,原始影像I(x,y)在不同的尺度e下,与高斯函数G(x,y,e)进行卷积,得到L(x,y,e),如下:
           

        以上的(x,y)是空间坐标, e,是尺度坐标,或尺度空间因子,e的大小决定平滑程度,大尺度对应图像的概貌特征,小尺度对应图像的细节特征。大的e值对应粗糙尺度(低分辨率),反之,对应精细尺度(高分辨率)。

        尺度,受e这个参数控制的表示。而不同的L(x,y,e)就构成了尺度空间,具体计算的时候,即使连续的高斯函数,都被离散为(一般为奇数大小)(2*k+1) *(2*k+1)矩阵,来和数字图像进行卷积运算。
        随着e的变化,建立起不同的尺度空间,或称之为建立起图像的高斯金字塔。

        像上述L(x,y,e) = G(x,y,e)*I(x,y)的操作,在进行高斯卷积时,整个图像就要遍历所有的像素进行卷积(边界点除外),于此,就造成了时间和空间上的很大浪费。
        为了更有效的在尺度空间检测到稳定的关键点,也为了缩小时间和空间复杂度,对上述的操作作了一个改建:即,提出了高斯差分尺度空间(DOG scale-space)。利用不同尺度的高斯差分与原始图像I(x,y)相乘 ,卷积生成。
           
        DOG算子计算简单,是尺度归一化的LOG算子的近似。

        ok,耐心点,咱们再来总结一下上述内容:
        1、高斯卷积
        在组建一组尺度空间后,再组建下一组尺度空间,对上一组尺度空间的最后一幅图像进行二分之一采样,得到下一组尺度空间的第一幅图像,然后进行像建立第一组尺度空间那样的操作,得到第二组尺度空间,公式定义为
             L(x,y,e) = G(x,y,e)*I(x,y)

        图像金字塔的构建:图像金字塔共O组,每组有S层,下一组的图像由上一组图像降采样得到,效果图,图A如下(左为上一组,右为下一组):

       
        2、高斯差分
        在尺度空间建立完毕后,为了能够找到稳定的关键点,采用高斯差分的方法来检测那些在局部位置的极值点,即采用俩个相邻的尺度中的图像相减,即公式定义为:
            D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y)
                     = L(x,y,ke) - L(x,y,e)

        效果图,图B


        SIFT的精妙之处在于采用图像金字塔的方法解决这一问题,我们可以把两幅图像想象成是连续的,分别以它们作为底面作四棱锥,就像金字塔,那么每一个 截面与原图像相似,那么两个金字塔中必然会有包含大小一致的物体的无穷个截面,但应用只能是离散的,所以我们只能构造有限层,层数越多当然越好,但处理时 间会相应增加,层数太少不行,因为向下采样的截面中可能找不到尺寸大小一致的两个物体的图像。

        咱们再来具体阐述下构造D(x,y,e)的详细步骤:
        1、首先采用不同尺度因子的高斯核对图像进行卷积以得到图像的不同尺度空间,将这一组图像作为金子塔图像的第一层。
        2、接着对第一层图像中的2倍尺度图像(相对于该层第一幅图像的2倍尺度)以2倍像素距离进行下采样来得到金子塔图像的第二层中的第一幅图像,对该图像采用不同尺度因子的高斯核进行卷积,以获得金字塔图像中第二层的一组图像。
        3、再以金字塔图像中第二层中的2倍尺度图像(相对于该层第一幅图像的2倍尺度)以2倍像素距离进行下采样来得到金字塔图像的第三层中的第一幅图像,对该图像采用不同尺度因子的高斯核进行卷积,以获得金字塔图像中第三层的一组图像。这样依次类推,从而获得了金字塔图像的每一层中的一组图像,如下图所示:

        4、对上图得到的每一层相邻的高斯图像相减,就得到了高斯差分图像,如下述第一幅图所示。下述第二幅图中的右列显示了将每组中相邻图像相减所生成的高斯差分图像的结果,限于篇幅,图中只给出了第一层和第二层高斯差分图像的计算(下述俩幅图统称为图2):

    图像金字塔的建立:对于一幅图像I,建立其在不同尺度(scale)的图像,也成为子八度(octave),这是为了scale-invariant,也就是在任何尺度都能够有对应的特征点,第一个子八度的scale为原图大小,后面每个octave为上一个octave降采样的结果,即原图的1/4(长宽分别减半),构成下一个子八度(高一层金字塔)

        5、因为高斯差分函数是归一化的高斯拉普拉斯函数的近似,所以可以从高斯差分金字塔分层结构提取出图像中的极值点作为候选的特征点。对DOG 尺度空间每个点与相邻尺度和相邻位置的点逐个进行比较,得到的局部极值位置即为特征点所处的位置和对应的尺度。

        二、检测关键点
        为了寻找尺度空间的极值点,每一个采样点要和它所有的相邻点比较,看其是否比它的图像域和尺度域的相邻点大或者小。如下图,图3所示,中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较,以确保在尺度空间和二维图像空间都检测到极值点。

        因为需要同相邻尺度进行比较,所以在一组高斯差分图像中只能检测到两个尺度的极值点(如上述第二幅图中右图的五角星标识),而其它尺度的极值点检测则需要在图像金字塔的上一层高斯差分图像中进行。依次类推,最终在图像金字塔中不同层的高斯差分图像中完成不同尺度极值的检测。
        当然这样产生的极值点并不都是稳定的特征点,因为某些极值点响应较弱,而且DOG算子会产生较强的边缘响应。

        三、关键点方向的分配
        为了使描述符具有旋转不变性,需要利用图像的局部特征为给每一个关键点分配一个方向。利用关键点邻域像素的梯度及方向分布的特性,可以得到梯度模值和方向如下:

       
        其中,尺度为每个关键点各自所在的尺度。
        在以关键点为中心的邻域窗口内采样,并用直方图统计邻域像素的梯度方向。梯度直方图的范围是0~360度,其中每10度一个方向,总共36个方向。
    直方图的峰值则代表了该关键点处邻域梯度的主方向,即作为该关键点的方向。

        在计算方向直方图时,需要用一个参数等于关键点所在尺度1.5倍的高斯权重窗对方向直方图进行加权,上图中用蓝色的圆形表示,中心处的蓝色较重,表示权值最大,边缘处颜色潜,表示权值小。如下图所示,该示例中为了简化给出了8方向的方向直方图计算结果,实际sift创始人David Lowe的原论文中采用36方向的直方图。

        方向直方图的峰值则代表了该特征点处邻域梯度的方向,以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性,只保留峰值大于主方向峰值80%的方向作为该关键点的辅方向。因此,对于同一梯度值的多个峰值的关键点位置,在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15%的关键点被赋予多个方向,但可以明显的提高关键点匹配的稳定性。
           
        至此,图像的关键点已检测完毕,每个关键点有三个信息:位置、所处尺度、方向。由此可以确定一个SIFT特征区域。

        四、特征点描述符
        通过以上步骤,对于每一个关键点,拥有三个信息:位置、尺度以及方向。接下来就是为每个关键点建立一个描述符,使其不随各种变化而改变,比如光照变化、视角变化等等。并且描述符应该有较高的独特性,以便于提高特征点正确匹配的概率。 
    首先将坐标轴旋转为关键点的方向,以确保旋转不变性。

        接下来以关键点为中心取8×8的窗口。
        上图,图5中左部分的中央黑点为当前关键点的位置,每个小格代表关键点邻域所在尺度空间的一个像素,箭头方向代表该像素的梯度方向,箭头长度代表梯度模值,图中蓝色的圈代表高斯加权的范围(越靠近关键点的像素梯度方向信息贡献越大)。
        然后在每4×4的小块上计算8个方向的梯度方向直方图,绘制每个梯度方向的累加值,即可形成一个种子点,如图5右部分所示。此图中一个关键点由2×2共4个种子点组成,每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力,同时对于含有定位误差的特征匹配也提供了较好的容错性。

        实际计算过程中,为了增强匹配的稳健性,Lowe建议对每个关键点使用4×4共16个种子点来描述,这样对于一个关键点就可以产生128个数据,即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响,再继续将特征向量的长度归一化,则可以进一步去除光照变化的影响。


        五、最后一步:当两幅图像的SIFT特征向量生成后,下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取上图中,图像A中的某个关键点,并找出其与图像B中欧式距离最近的前两个关键点,在这两个关键点中,如果最近的距离除以次近的距离少于某个比例阈值,则接受这一对匹配点。降低这个比例阈值,SIFT匹配点数目会减少,但更加稳定。关于sift 算法的更多理论介绍请参看此文:http://blog.csdn.net/abcjennifer/article/details/7639681

     

    sift算法的逐步c实现
        ok,上文搅了那么多的理论,如果你没有看懂它,咋办列?没关系,下面,咱们来一步一步实现此sift算法,即使你没有看到上述的理论,慢慢的,你也会明白sift算法到底是怎么一回事,sift算法到底是怎么实现的...。
       yeah,请看:

    前期工作:
    在具体编写核心函数之前,得先做几个前期的准备工作:

    0、头文件:

    1. #ifdef _CH_  
    2. #pragma package <opencv>  
    3. #endif  
    4.   
    5. #ifndef _EiC  
    6. #include <stdio.h>  
    7.   
    8. #include "stdlib.h"  
    9. #include "string.h"   
    10. #include "malloc.h"   
    11. #include "math.h"   
    12. #include <assert.h>  
    13. #include <ctype.h>  
    14. #include <time.h>  
    15. #include <cv.h>  
    16. #include <cxcore.h>  
    17. #include <highgui.h>  
    18. #include <vector>  
    19. #endif  
    20.   
    21. #ifdef _EiC  
    22. #define WIN32  
    23. #endif  

    1、定义几个宏,及变量,以免下文函数中,突然冒出一个变量,而您却不知道怎么一回事:
    1. #define NUMSIZE 2  
    2. #define GAUSSKERN 3.5  
    3. #define PI 3.14159265358979323846  
    4.   
    5. //Sigma of base image -- See D.L.'s paper.  
    6. #define INITSIGMA 0.5  
    7. //Sigma of each octave -- See D.L.'s paper.  
    8. #define SIGMA sqrt(3)//1.6//  
    9.   
    10. //Number of scales per octave.  See D.L.'s paper.  
    11. #define SCALESPEROCTAVE 2  
    12. #define MAXOCTAVES 4  
    13. int     numoctaves;  
    14.   
    15. #define CONTRAST_THRESHOLD   0.02  
    16. #define CURVATURE_THRESHOLD  10.0  
    17. #define DOUBLE_BASE_IMAGE_SIZE 1  
    18. #define peakRelThresh 0.8  
    19. #define LEN 128  
    20.   
    21. // temporary storage  
    22. CvMemStorage* storage = 0;   

    2、然后,咱们还得,声明几个变量,以及建几个 数据结构(数据结构是一切程序事物的基础麻,:D。):
    1. //Data structure for a float image.  
    2. typedef struct ImageSt {        /*金字塔每一层*/  
    3.    
    4.  float levelsigma;  
    5.  int levelsigmalength;  
    6.  float absolute_sigma;  
    7.  CvMat *Level;       //CvMat是OPENCV的矩阵类,其元素可以是图像的象素值         
    8. } ImageLevels;  
    9.   
    10. typedef struct ImageSt1 {      /*金字塔每一阶梯*/  
    11.  int row, col;          //Dimensions of image.   
    12.  float subsample;  
    13.  ImageLevels *Octave;                
    14. } ImageOctaves;  
    15.   
    16. ImageOctaves *DOGoctaves;        
    17. //DOG pyr,DOG算子计算简单,是尺度归一化的LoG算子的近似。  
    18.     
    19. ImageOctaves *mag_thresh ;  
    20. ImageOctaves *mag_pyr ;  
    21. ImageOctaves *grad_pyr ;  
    22.   
    23. //keypoint数据结构,Lists of keypoints are linked by the "next" field.  
    24. typedef struct KeypointSt   
    25. {  
    26.  float row, col; /* 反馈回原图像大小,特征点的位置 */  
    27.  float sx,sy;    /* 金字塔中特征点的位置*/  
    28.  int octave,level;/*金字塔中,特征点所在的阶梯、层次*/  
    29.    
    30.  float scale, ori,mag; /*所在层的尺度sigma,主方向orientation (range [-PI,PI]),以及幅值*/  
    31.  float *descrip;       /*特征描述字指针:128维或32维等*/  
    32.  struct KeypointSt *next;/* Pointer to next keypoint in list. */  
    33. } *Keypoint;  
    34.   
    35. //定义特征点具体变量  
    36. Keypoint keypoints=NULL;      //用于临时存储特征点的位置等  
    37. Keypoint keyDescriptors=NULL; //用于最后的确定特征点以及特征描述字  

    3、声明几个图像的基本处理函数:

    1. CvMat * halfSizeImage(CvMat * im);     //缩小图像:下采样  
    2. CvMat * doubleSizeImage(CvMat * im);   //扩大图像:最近临方法  
    3. CvMat * doubleSizeImage2(CvMat * im);  //扩大图像:线性插值  
    4. float getPixelBI(CvMat * im, float col, float row);//双线性插值函数  
    5. void normalizeVec(float* vec, int dim);//向量归一化    
    6. CvMat* GaussianKernel2D(float sigma);  //得到2维高斯核  
    7. void normalizeMat(CvMat* mat) ;        //矩阵归一化  
    8. float* GaussianKernel1D(float sigma, int dim) ; //得到1维高斯核  
    9.   
    10. //在具体像素处宽度方向进行高斯卷积  
    11. float ConvolveLocWidth(float* kernel, int dim, CvMat * src, int x, int y) ;    
    12. //在整个图像宽度方向进行1D高斯卷积  
    13. void Convolve1DWidth(float* kern, int dim, CvMat * src, CvMat * dst) ;         
    14. //在具体像素处高度方向进行高斯卷积  
    15. float ConvolveLocHeight(float* kernel, int dim, CvMat * src, int x, int y) ;   
    16. //在整个图像高度方向进行1D高斯卷积  
    17. void Convolve1DHeight(float* kern, int dim, CvMat * src, CvMat * dst);       
    18. //用高斯函数模糊图像    
    19. int BlurImage(CvMat * src, CvMat * dst, float sigma) ;            
                
    算法核心
       本程序中,sift算法被分为以下五个步骤及其相对应的函数(可能表述与上,或与前俩篇文章有所偏差,但都一个意思):
    1. //SIFT算法第一步:图像预处理  
    2. CvMat *ScaleInitImage(CvMat * im) ;                  //金字塔初始化  
    3.   
    4. //SIFT算法第二步:建立高斯金字塔函数  
    5. ImageOctaves* BuildGaussianOctaves(CvMat * image) ;  //建立高斯金字塔  
    6.   
    7. //SIFT算法第三步:特征点位置检测,最后确定特征点的位置  
    8. int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);  
    9. void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr);  
    10.   
    11. //SIFT算法第四步:计算高斯图像的梯度方向和幅值,计算各个特征点的主方向  
    12. void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);  
    13.   
    14. int FindClosestRotationBin (int binCount, float angle);  //进行方向直方图统计  
    15. void AverageWeakBins (double* bins, int binCount);       //对方向直方图滤波  
    16. //确定真正的主方向  
    17. bool InterpolateOrientation (double left, double middle,double right, double *degreeCorrection, double *peakValue);  
    18. //确定各个特征点处的主方向函数  
    19. void AssignTheMainOrientation(int numoctaves, ImageOctaves *GaussianPyr,ImageOctaves *mag_pyr,ImageOctaves *grad_pyr);  
    20. //显示主方向  
    21. void DisplayOrientation (IplImage* image, ImageOctaves *GaussianPyr);  
    22.   
    23. //SIFT算法第五步:抽取各个特征点处的特征描述字  
    24. void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);  
    25.   
    26. //为了显示图象金字塔,而作的图像水平、垂直拼接  
    27. CvMat* MosaicHorizen( CvMat* im1, CvMat* im2 );  
    28. CvMat* MosaicVertical( CvMat* im1, CvMat* im2 );  
    29.   
    30. //特征描述点,网格    
    31. #define GridSpacing 4  


    主体实现
        ok,以上所有的工作都就绪以后,那么接下来,咱们就先来编写main函数,因为你一看主函数之后,你就立马能发现sift算法的工作流程及其原理了。
        (主函数中涉及到的函数,下一篇文章:一、教你一步一步用c语言实现sift算法、下,咱们自会一个一个编写):

    1. int main( void )  
    2. {  
    3.  //声明当前帧IplImage指针  
    4.  IplImage* src = NULL;   
    5.  IplImage* image1 = NULL;   
    6.  IplImage* grey_im1 = NULL;   
    7.  IplImage* DoubleSizeImage = NULL;  
    8.    
    9.  IplImage* mosaic1 = NULL;   
    10.  IplImage* mosaic2 = NULL;   
    11.    
    12.  CvMat* mosaicHorizen1 = NULL;  
    13.  CvMat* mosaicHorizen2 = NULL;  
    14.  CvMat* mosaicVertical1 = NULL;  
    15.    
    16.  CvMat* image1Mat = NULL;  
    17.  CvMat* tempMat=NULL;  
    18.    
    19.  ImageOctaves *Gaussianpyr;  
    20.  int rows,cols;  
    21.   
    22. #define Im1Mat(ROW,COL) ((float *)(image1Mat->data.fl + image1Mat->step/sizeof(float) *(ROW)))[(COL)]  
    23.    
    24.  //灰度图象像素的数据结构  
    25. #define Im1B(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3]  
    26. #define Im1G(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+1]  
    27. #define Im1R(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+2]  
    28.    
    29.  storage = cvCreateMemStorage(0);   
    30.    
    31.  //读取图片  
    32.  if( (src = cvLoadImage( "street1.jpg", 1)) == 0 )  // test1.jpg einstein.pgm back1.bmp  
    33.   return -1;  
    34.   
    35.  //为图像分配内存   
    36.  image1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,3);  
    37.  grey_im1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,1);  
    38.  DoubleSizeImage = cvCreateImage(cvSize(2*(src->width), 2*(src->height)),  IPL_DEPTH_8U,3);  
    39.   
    40.  //为图像阵列分配内存,假设两幅图像的大小相同,tempMat跟随image1的大小  
    41.  image1Mat = cvCreateMat(src->height, src->width, CV_32FC1);  
    42.  //转化成单通道图像再处理  
    43.  cvCvtColor(src, grey_im1, CV_BGR2GRAY);  
    44.  //转换进入Mat数据结构,图像操作使用的是浮点型操作  
    45.  cvConvert(grey_im1, image1Mat);  
    46.    
    47.  double t = (double)cvGetTickCount();  
    48.  //图像归一化  
    49.  cvConvertScale( image1Mat, image1Mat, 1.0/255, 0 );  
    50.    
    51.  int dim = min(image1Mat->rows, image1Mat->cols);  
    52.  numoctaves = (int) (log((double) dim) / log(2.0)) - 2;    //金字塔阶数  
    53.  numoctaves = min(numoctaves, MAXOCTAVES);  
    54.   
    55.  //SIFT算法第一步,预滤波除噪声,建立金字塔底层  
    56.  tempMat = ScaleInitImage(image1Mat) ;  
    57.  //SIFT算法第二步,建立Guassian金字塔和DOG金字塔  
    58.  Gaussianpyr = BuildGaussianOctaves(tempMat) ;  
    59.    
    60.  t = (double)cvGetTickCount() - t;  
    61.  printf( "the time of build Gaussian pyramid and DOG pyramid is %.1f/n", t/(cvGetTickFrequency()*1000.) );  
    62.    
    63. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
    64.  //显示高斯金字塔  
    65.  for (int i=0; i<numoctaves;i++)  
    66.  {  
    67.   if (i==0)  
    68.   {  
    69.    mosaicHorizen1=MosaicHorizen( (Gaussianpyr[0].Octave)[0].Level, (Gaussianpyr[0].Octave)[1].Level );  
    70.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
    71.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[0].Octave)[j].Level );  
    72.    for ( j=0;j<NUMSIZE;j++)  
    73.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    74.   }  
    75.   else if (i==1)  
    76.   {  
    77.    mosaicHorizen2=MosaicHorizen( (Gaussianpyr[1].Octave)[0].Level, (Gaussianpyr[1].Octave)[1].Level );  
    78.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
    79.     mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (Gaussianpyr[1].Octave)[j].Level );  
    80.    for ( j=0;j<NUMSIZE;j++)  
    81.     mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
    82.    mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
    83.   }  
    84.   else  
    85.   {  
    86.    mosaicHorizen1=MosaicHorizen( (Gaussianpyr[i].Octave)[0].Level, (Gaussianpyr[i].Octave)[1].Level );  
    87.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
    88.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[i].Octave)[j].Level );  
    89.    for ( j=0;j<NUMSIZE;j++)  
    90.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    91.    mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
    92.   }  
    93.  }  
    94.  mosaic1 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
    95.  cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0, 0 );  
    96.  cvConvertScaleAbs( mosaicVertical1, mosaic1, 1, 0 );  
    97.    
    98.  //  cvSaveImage("GaussianPyramid of me.jpg",mosaic1);  
    99.  cvNamedWindow("mosaic1",1);  
    100.  cvShowImage("mosaic1", mosaic1);  
    101.  cvWaitKey(0);  
    102.  cvDestroyWindow("mosaic1");  
    103.  //显示DOG金字塔  
    104.  for ( i=0; i<numoctaves;i++)  
    105.  {  
    106.   if (i==0)  
    107.   {  
    108.    mosaicHorizen1=MosaicHorizen( (DOGoctaves[0].Octave)[0].Level, (DOGoctaves[0].Octave)[1].Level );  
    109.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
    110.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[0].Octave)[j].Level );  
    111.    for ( j=0;j<NUMSIZE;j++)  
    112.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    113.   }  
    114.   else if (i==1)  
    115.   {  
    116.    mosaicHorizen2=MosaicHorizen( (DOGoctaves[1].Octave)[0].Level, (DOGoctaves[1].Octave)[1].Level );  
    117.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
    118.     mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (DOGoctaves[1].Octave)[j].Level );  
    119.    for ( j=0;j<NUMSIZE;j++)  
    120.     mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
    121.    mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
    122.   }  
    123.   else  
    124.   {  
    125.    mosaicHorizen1=MosaicHorizen( (DOGoctaves[i].Octave)[0].Level, (DOGoctaves[i].Octave)[1].Level );  
    126.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
    127.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[i].Octave)[j].Level );  
    128.    for ( j=0;j<NUMSIZE;j++)  
    129.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
    130.    mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
    131.   }  
    132.  }  
    133.  //考虑到DOG金字塔各层图像都会有正负,所以,必须寻找最负的,以将所有图像抬高一个台阶去显示  
    134.  double min_val=0;  
    135.  double max_val=0;  
    136.  cvMinMaxLoc( mosaicVertical1, &min_val, &max_val,NULL, NULL, NULL );  
    137.  if ( min_val<0.0 )  
    138.   cvAddS( mosaicVertical1, cvScalarAll( (-1.0)*min_val ), mosaicVertical1, NULL );  
    139.  mosaic2 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
    140.  cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0/(max_val-min_val), 0 );  
    141.  cvConvertScaleAbs( mosaicVertical1, mosaic2, 1, 0 );  
    142.    
    143.  //  cvSaveImage("DOGPyramid of me.jpg",mosaic2);  
    144.  cvNamedWindow("mosaic1",1);  
    145.  cvShowImage("mosaic1", mosaic2);  
    146.  cvWaitKey(0);  
    147.    
    148.  //SIFT算法第三步:特征点位置检测,最后确定特征点的位置  
    149.  int keycount=DetectKeypoint(numoctaves, Gaussianpyr);  
    150.  printf("the keypoints number are %d ;/n", keycount);  
    151.  cvCopy(src,image1,NULL);  
    152.  DisplayKeypointLocation( image1 ,Gaussianpyr);  
    153.    
    154.  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
    155.  cvNamedWindow("image1",1);  
    156.  cvShowImage("image1", DoubleSizeImage);  
    157.  cvWaitKey(0);    
    158.  cvDestroyWindow("image1");  
    159.    
    160.  //SIFT算法第四步:计算高斯图像的梯度方向和幅值,计算各个特征点的主方向  
    161.  ComputeGrad_DirecandMag(numoctaves, Gaussianpyr);  
    162.  AssignTheMainOrientation( numoctaves, Gaussianpyr,mag_pyr,grad_pyr);  
    163.  cvCopy(src,image1,NULL);  
    164.  DisplayOrientation ( image1, Gaussianpyr);  
    165.    
    166.  //  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
    167.  cvNamedWindow("image1",1);  
    168.  //  cvResizeWindow("image1", 2*(image1->width), 2*(image1->height) );  
    169.  cvShowImage("image1", image1);  
    170.  cvWaitKey(0);  
    171.    
    172.  //SIFT算法第五步:抽取各个特征点处的特征描述字  
    173.  ExtractFeatureDescriptors( numoctaves, Gaussianpyr);  
    174.  cvWaitKey(0);  
    175.    
    176.  //销毁窗口  
    177.  cvDestroyWindow("image1");  
    178.  cvDestroyWindow("mosaic1");  
    179.  //释放图像  
    180.  cvReleaseImage(&image1);  
    181.  cvReleaseImage(&grey_im1);  
    182.  cvReleaseImage(&mosaic1);  
    183.  cvReleaseImage(&mosaic2);  
    184.  return 0;  
    185. }  

    更多见下文:一、教你一步一步用c语言实现sift算法、下。本文完。

    展开全文
  • 九之再续:教你一步一步c语言实现sift算法、下

    万次阅读 热门讨论 2011-03-13 13:10:00
    教你一步一步c语言实现sift算法、下作者:July、二零一一年三月十二日出处:http://blog.csdn.net/v_JULY_v。...------------------------本文接上,教你一步一步c语言实现sift算法、上,而来:函数编写

                          教你一步一步用c语言实现sift算法、下


    作者:July、二零一一年三月十二日
    出处:http://blog.csdn.net/v_JULY_v
    参考:Rob Hess维护的sift 库
    环境:windows xp+vc6.0
    条件:c语言实现。
    说明:本BLOG内会陆续一一实现所有经典算法。
    ------------------------


    本文接上,教你一步一步用c语言实现sift算法、上,而来:
    函数编写
        ok,接上文,咱们一个一个的来编写main函数中所涉及到所有函数,这也是本文的关键部分:

    五个步骤

        ok,接下来,进入重点部分,咱们依据上文介绍的sift算法的几个步骤,来一一实现这些函数。
        为了版述清晰,再贴一下,主函数,顺便再加强下对sift 算法的五个步骤的认识:

    1、SIFT算法第一步:图像预处理
    CvMat *ScaleInitImage(CvMat * im) ;                  //金字塔初始化
    2、SIFT算法第二步:建立高斯金字塔函数
    ImageOctaves* BuildGaussianOctaves(CvMat * image) ;  //建立高斯金字塔
    3、SIFT算法第三步:特征点位置检测,最后确定特征点的位置
    int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);
    4、SIFT算法第四步:计算高斯图像的梯度方向和幅值,计算各个特征点的主方向
    void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);
    5、SIFT算法第五步:抽取各个特征点处的特征描述字
    void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);

    ok,接下来一一具体实现这几个函数:
    SIFT算法第一步
        SIFT算法第一步:扩大图像,预滤波剔除噪声,得到金字塔的最底层-第一阶的第一层:
      

    SIFT算法第二步
        SIFT第二步,建立Gaussian金字塔,给定金字塔第一阶第一层图像后,计算高斯金字塔其他尺度图像,
    每一阶的数目由变量SCALESPEROCTAVE决定,给定一个基本图像,计算它的高斯金字塔图像,返回外部向量是阶梯指针,内部向量是每一个阶梯内部的不同尺度图像。
      

    SIFT算法第三步
        SIFT算法第三步,特征点位置检测,最后确定特征点的位置检测DOG金字塔中的局部最大值,找到之后,还要经过两个检验才能确认为特征点:一是它必须有明显的差异,二是他不应该是边缘点,(也就是说,在极值点处的主曲率比应该小于某一个阈值)。

    SIFT算法第四步

    SIFT算法第五步
        SIFT算法第五步:抽取各个特征点处的特征描述字,确定特征点的描述字。描述字是Patch网格内梯度方向的描述,旋转网格到主方向,插值得到网格处梯度值。
    一个特征点可以用2*2*8=32维的向量,也可以用4*4*8=128维的向量更精确的进行描述。

     

    ok,为了版述清晰,再贴一下上文所述的主函数(注,上文已贴出,此是为了版述清晰,重复造轮):

     

    最后,再看一下,运行效果(图中美女为老乡+朋友,何姐08年照):

    完。 

    updated

        有很多朋友都在本文评论下要求要本程序的完整源码包(注:本文代码未贴全,复制粘贴编译肯定诸多错误),但由于时隔太久,这份代码我自己也找不到了,不过,我可以提供一份sift + KD + BBF,且可以编译正确的代码供大家参考学习,有pudn帐号的朋友可以前去下载:http://www.pudn.com/downloads340/sourcecode/graph/texture_mapping/detail1486667.html没有pudn账号的同学请加群:169056165,验证信息:sift,至群共享下载,然后用两幅不同的图片做了下匹配(当然,运行结果显示是不匹配的),效果还不错:http://weibo.com/1580904460/yDmzAEwcV#1348475194313! July、二零一二年十月十一日。

    展开全文
  • [纯C语言 + Win32 API]一步一步写个围棋程序之三:个菜单
  • 九之再续:教你一步一步c语言实现sift算法、上

    万次阅读 多人点赞 2011-03-13 09:32:00
    教你一步一步c语言实现sift算法、上作者:July、二零一一年三月十二日出处:http://blog.csdn.net/v_JULY_v参考:Rob Hess维护的sift 库环境:windows xp+vc6.0条件:c语言实现。说明:本BLOG内会陆续一一实现所有...

                         教你一步一步用c语言实现sift算法、上

    作者:July、二零一一年三月十二日
    出处:http://blog.csdn.net/v_JULY_v
    参考:Rob Hess维护的sift 库
    环境:windows xp+vc6.0
    条件:c语言实现。
    说明:本BLOG内会陆续一一实现所有经典算法。
    ------------------------

    引言:
        在我写的关于sift算法的前倆篇文章里头,已经对sift算法有了初步的介绍:九、图像特征提取与匹配之SIFT算法,而后在:九(续)、sift算法的编译与实现里,我也简单记录下了如何利用opencv,gsl等库编译运行sift程序。
        但据一朋友表示,是否能用c语言实现sift算法,同时,尽量不用到opencv,gsl等第三方库之类的东西。而且,Rob Hess维护的sift 库,也不好懂,有的人根本搞不懂是怎么一回事。
        那么本文,就教你如何利用c语言一步一步实现sift算法,同时,你也就能真正明白sift算法到底是怎么一回事了。

        ok,先看一下,本程序最终运行的效果图,sift 算法分为五个步骤(下文详述),对应以下第二--第六幅图

     

    sift算法的步骤
        要实现一个算法,首先要完全理解这个算法的原理或思想。咱们先来简单了解下,什么叫sift算法:
        sift,尺度不变特征转换,是一种电脑视觉的算法用来侦测与描述影像中的局部性特征,它在空间尺度中寻找极值点,并提取出其位置、尺度、旋转不变量,此算法由 David Lowe 在1999年所发表,2004年完善总结。

        所谓,Sift算法就是用不同尺度(标准差)的高斯函数对图像进行平滑,然后比较平滑后图像的差别,
    差别大的像素就是特征明显的点。

        以下是sift算法的五个步骤:
        一、建立图像尺度空间(或高斯金字塔),并检测极值点

        首先建立尺度空间,要使得图像具有尺度空间不变形,就要建立尺度空间,sift算法采用了高斯函数来建立尺度空间,高斯函数公式为:
           
        上述公式G(x,y,e),即为尺度可变高斯函数。

        而,一个图像的尺度空间L(x,y,e) ,定义为原始图像I(x,y)与上述的一个可变尺度的2维高斯函数G(x,y,e) 卷积运算。
        即,原始影像I(x,y)在不同的尺度e下,与高斯函数G(x,y,e)进行卷积,得到L(x,y,e),如下:
           

        以上的(x,y)是空间坐标, e,是尺度坐标,或尺度空间因子,e的大小决定平滑程度,大尺度对应图像的概貌特征,小尺度对应图像的细节特征。大的e值对应粗糙尺度(低分辨率),反之,对应精细尺度(高分辨率)。

        尺度,受e这个参数控制的表示。而不同的L(x,y,e)就构成了尺度空间,具体计算的时候,即使连续的高斯函数,都被离散为(一般为奇数大小)(2*k+1) *(2*k+1)矩阵,来和数字图像进行卷积运算。
        随着e的变化,建立起不同的尺度空间,或称之为建立起图像的高斯金字塔。

        像上述L(x,y,e) = G(x,y,e)*I(x,y)的操作,在进行高斯卷积时,整个图像就要遍历所有的像素进行卷积(边界点除外),于此,就造成了时间和空间上的很大浪费。
        为了更有效的在尺度空间检测到稳定的关键点,也为了缩小时间和空间复杂度,对上述的操作作了一个改建:即,提出了高斯差分尺度空间(DOG scale-space)。利用不同尺度的高斯差分与原始图像I(x,y)相乘 ,卷积生成。
           
        DOG算子计算简单,是尺度归一化的LOG算子的近似。

        ok,耐心点,咱们再来总结一下上述内容:
        1、高斯卷积
        在组建一组尺度空间后,再组建下一组尺度空间,对上一组尺度空间的最后一幅图像进行二分之一采样,得到下一组尺度空间的第一幅图像,然后进行像建立第一组尺度空间那样的操作,得到第二组尺度空间,公式定义为
             L(x,y,e) = G(x,y,e)*I(x,y)

        图像金字塔的构建:图像金字塔共O组,每组有S层,下一组的图像由上一组图像降采样得到,效果图,图A如下(左为上一组,右为下一组):

       
        2、高斯差分
        在尺度空间建立完毕后,为了能够找到稳定的关键点,采用高斯差分的方法来检测那些在局部位置的极值点,即采用俩个相邻的尺度中的图像相减,即公式定义为:
            D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y)
                     = L(x,y,ke) - L(x,y,e)

        效果图,图B


        SIFT的精妙之处在于采用图像金字塔的方法解决这一问题,我们可以把两幅图像想象成是连续的,分别以它们作为底面作四棱锥,就像金字塔,那么每一个 截面与原图像相似,那么两个金字塔中必然会有包含大小一致的物体的无穷个截面,但应用只能是离散的,所以我们只能构造有限层,层数越多当然越好,但处理时 间会相应增加,层数太少不行,因为向下采样的截面中可能找不到尺寸大小一致的两个物体的图像。

        咱们再来具体阐述下构造D(x,y,e)的详细步骤:
        1、首先采用不同尺度因子的高斯核对图像进行卷积以得到图像的不同尺度空间,将这一组图像作为金子塔图像的第一层。
        2、接着对第一层图像中的2倍尺度图像(相对于该层第一幅图像的2倍尺度)以2倍像素距离进行下采样来得到金子塔图像的第二层中的第一幅图像,对该图像采用不同尺度因子的高斯核进行卷积,以获得金字塔图像中第二层的一组图像。
        3、再以金字塔图像中第二层中的2倍尺度图像(相对于该层第一幅图像的2倍尺度)以2倍像素距离进行下采样来得到金字塔图像的第三层中的第一幅图像,对该图像采用不同尺度因子的高斯核进行卷积,以获得金字塔图像中第三层的一组图像。这样依次类推,从而获得了金字塔图像的每一层中的一组图像,如下图所示:

        4、对上图得到的每一层相邻的高斯图像相减,就得到了高斯差分图像,如下述第一幅图所示。下述第二幅图中的右列显示了将每组中相邻图像相减所生成的高斯差分图像的结果,限于篇幅,图中只给出了第一层和第二层高斯差分图像的计算(下述俩幅图统称为图2):

    图像金字塔的建立:对于一幅图像I,建立其在不同尺度(scale)的图像,也成为子八度(octave),这是为了scale-invariant,也就是在任何尺度都能够有对应的特征点,第一个子八度的scale为原图大小,后面每个octave为上一个octave降采样的结果,即原图的1/4(长宽分别减半),构成下一个子八度(高一层金字塔)

        5、因为高斯差分函数是归一化的高斯拉普拉斯函数的近似,所以可以从高斯差分金字塔分层结构提取出图像中的极值点作为候选的特征点。对DOG 尺度空间每个点与相邻尺度和相邻位置的点逐个进行比较,得到的局部极值位置即为特征点所处的位置和对应的尺度。

        二、检测关键点
        为了寻找尺度空间的极值点,每一个采样点要和它所有的相邻点比较,看其是否比它的图像域和尺度域的相邻点大或者小。如下图,图3所示,中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较,以确保在尺度空间和二维图像空间都检测到极值点。

        因为需要同相邻尺度进行比较,所以在一组高斯差分图像中只能检测到两个尺度的极值点(如上述第二幅图中右图的五角星标识),而其它尺度的极值点检测则需要在图像金字塔的上一层高斯差分图像中进行。依次类推,最终在图像金字塔中不同层的高斯差分图像中完成不同尺度极值的检测。
        当然这样产生的极值点并不都是稳定的特征点,因为某些极值点响应较弱,而且DOG算子会产生较强的边缘响应。

        三、关键点方向的分配
        为了使描述符具有旋转不变性,需要利用图像的局部特征为给每一个关键点分配一个方向。利用关键点邻域像素的梯度及方向分布的特性,可以得到梯度模值和方向如下:

       
        其中,尺度为每个关键点各自所在的尺度。
        在以关键点为中心的邻域窗口内采样,并用直方图统计邻域像素的梯度方向。梯度直方图的范围是0~360度,其中每10度一个方向,总共36个方向。
    直方图的峰值则代表了该关键点处邻域梯度的主方向,即作为该关键点的方向。

        在计算方向直方图时,需要用一个参数等于关键点所在尺度1.5倍的高斯权重窗对方向直方图进行加权,上图中用蓝色的圆形表示,中心处的蓝色较重,表示权值最大,边缘处颜色潜,表示权值小。如下图所示,该示例中为了简化给出了8方向的方向直方图计算结果,实际sift创始人David Lowe的原论文中采用36方向的直方图。

        方向直方图的峰值则代表了该特征点处邻域梯度的方向,以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性,只保留峰值大于主方向峰值80%的方向作为该关键点的辅方向。因此,对于同一梯度值的多个峰值的关键点位置,在相同位置和尺度将会有多个关键点被创建但方向不同。仅有15%的关键点被赋予多个方向,但可以明显的提高关键点匹配的稳定性。
           
        至此,图像的关键点已检测完毕,每个关键点有三个信息:位置、所处尺度、方向。由此可以确定一个SIFT特征区域。

        四、特征点描述符
        通过以上步骤,对于每一个关键点,拥有三个信息:位置、尺度以及方向。接下来就是为每个关键点建立一个描述符,使其不随各种变化而改变,比如光照变化、视角变化等等。并且描述符应该有较高的独特性,以便于提高特征点正确匹配的概率。 
    首先将坐标轴旋转为关键点的方向,以确保旋转不变性。

        接下来以关键点为中心取8×8的窗口。
        上图,图5中左部分的中央黑点为当前关键点的位置,每个小格代表关键点邻域所在尺度空间的一个像素,箭头方向代表该像素的梯度方向,箭头长度代表梯度模值,图中蓝色的圈代表高斯加权的范围(越靠近关键点的像素梯度方向信息贡献越大)。
        然后在每4×4的小块上计算8个方向的梯度方向直方图,绘制每个梯度方向的累加值,即可形成一个种子点,如图5右部分所示。此图中一个关键点由2×2共4个种子点组成,每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力,同时对于含有定位误差的特征匹配也提供了较好的容错性。

        实际计算过程中,为了增强匹配的稳健性,Lowe建议对每个关键点使用4×4共16个种子点来描述,这样对于一个关键点就可以产生128个数据,即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响,再继续将特征向量的长度归一化,则可以进一步去除光照变化的影响。


        五、最后一步:当两幅图像的SIFT特征向量生成后,下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取上图中,图像A中的某个关键点,并找出其与图像B中欧式距离最近的前两个关键点,在这两个关键点中,如果最近的距离除以次近的距离少于某个比例阈值,则接受这一对匹配点。降低这个比例阈值,SIFT匹配点数目会减少,但更加稳定。关于sift 算法的更多理论介绍请参看此文:http://blog.csdn.net/abcjennifer/article/details/7639681

     

    sift算法的逐步c实现
        ok,上文搅了那么多的理论,如果你没有看懂它,咋办列?没关系,下面,咱们来一步一步实现此sift算法,即使你没有看到上述的理论,慢慢的,你也会明白sift算法到底是怎么一回事,sift算法到底是怎么实现的...。
       yeah,请看:

    前期工作:
    在具体编写核心函数之前,得先做几个前期的准备工作:

    0、头文件:
    1、定义几个宏,及变量,以免下文函数中,突然冒出一个变量,而您却不知道怎么一回事:


    2、然后,咱们还得,声明几个变量,以及建几个数据结构(数据结构是一切程序事物的基础麻,:D。):

    3、声明几个图像的基本处理函数:
                
    算法核心
       本程序中,sift算法被分为以下五个步骤及其相对应的函数(可能表述与上,或与前俩篇文章有所偏差,但都一个意思):


    主体实现
        ok,以上所有的工作都就绪以后,那么接下来,咱们就先来编写main函数,因为你一看主函数之后,你就立马能发现sift算法的工作流程及其原理了。
        (主函数中涉及到的函数,下一篇文章:一、教你一步一步用c语言实现sift算法、下,咱们自会一个一个编写):

    更多见下文:一、教你一步一步用c语言实现sift算法、下。本文完。

    展开全文
  • struct sha256 { u32 block[16]; //加密measage 512bit u32 hash[8]; //hash的结果 u64 hash_length;//总共hash的byte数 u8 offset; //一个update未对齐Word(4字节)的字节数 u8 index; //当前已经写到block的...
  • DES算法加密C语言实现

    千次阅读 多人点赞 2018-11-24 19:34:57
    这几天把老师布置作业写了写,主要是DES算法加密,随便写了写,这里是...代码每一步都有解释,可能实现的有些麻烦,但是易懂,代码如下: ///DES算法加密 #include<cstdio> #include<cstring> const ...

    这几天把老师布置作业写了写,主要是DES算法加密,随便写了写,这里是实现对明文为8字节的加密以及对加密产生的密文进行解密,而不能直接输入密文。。。

    什么是DES算法,请自行百度学习这里只给出实现代码。

    代码每一步都有解释,可能实现的有些麻烦,但是易懂,代码如下:

    ///DES算法加密
    #include<cstdio>
    #include<cstring>
    const int maxn=2e4;
    char mingwen[maxn],miwen[maxn],miyao[maxn];///输入的明文密文密钥
    int smbit[maxn],miyaobit[maxn],mibit[maxn],zhongbit[maxn];///明文,密文,密钥的二进制,与转换中需要作为媒介的
    int mip[maxn],rr[maxn],SYa[maxn],er[10];///mip经过ip初始置换得到rrE盒扩展后得到,SYas盒压缩得到,er在s盒的转换过程中使用
    int Left[maxn],Right[maxn];///左右两部分
    int miyaopc1[maxn],miyaopc2[maxn],c[maxn],d[maxn];///c,d密钥置换得到
    int zimi[16][maxn];///产生的子密钥
    int LL[17][64],RR[17][64],cnt;///用来记录每一次的左右部分
    int times,xz[]= {1,1,2,2,2,2,2,2,1,2,2,2,2,2,2,1}; ///循环左移位,times用来计算是第几次移位
    void init()///初始化
    {
        times=cnt=0;
        memset(c,0,sizeof(c));
        memset(d,0,sizeof(d));
        memset(LL,0,sizeof(LL));
        memset(RR,0,sizeof(RR));
        memset(mip,0,sizeof(mip));
        memset(SYa,0,sizeof(SYa));
        memset(zimi,0,sizeof(zimi));
        memset(Left,0,sizeof(Left));
        memset(Right,0,sizeof(Right));
        memset(miyao,0,sizeof(miyao));
        memset(mingwen,0,sizeof(mingwen));
        memset(miyaopc1,0,sizeof(miyaopc1));
        memset(miyaopc2,0,sizeof(miyaopc2));
        memset(miyaobit,0,sizeof(miyaobit));
    }
    const int IP[]=///ip置换表,用1维存的原因是方便置换。下面的都一样
    {
        58,50,42,34,26,18,10,2,
        60,52,44,36,28,20,12,4,
        62,54,46,38,30,22,14,6,
        64,56,48,40,32,24,16,8,
        57,49,41,33,25,17, 9,1,
        59,51,43,35,27,19,11,3,
        61,53,45,37,29,21,13,5,
        63,55,47,39,31,23,15,7,
    };
    const int E[]=///扩展E盒
    {
        32,1,2,3,4,5,
        4,5,6,7,8,9,
        8,9,10,11,12,13,
        12,13,14,15,16,17,
        16,17,18,19,20,21,
        20,21,22,23,24,25,
        24,25,26,27,28,29,
        28,29,30,31,32,1,
    };
    const int P[]=///p盒置换表
    {
        16,7,20,21,29,12,28,17,1,15,23,26,5,18,31,10,
        2,8,24,14,32,27,3,9,19,13,30,6,22,11,4,25,
    };
    const int PC_1[]=///密钥置换选择PC_1盒
    {
        57,49,41,33,25,17,9,1,58,50,42,34,26,18,
        10,2,59,51,43,35,27,19,11,3,60,52,44,36,
        63,55,47,39,31,23,15,7,62,54,46,38,30,22,
        14,6,61,53,45,37,29,21,13,5,28,20,12,4,
    };
    const int PC_2[]=///密钥置换选择PC_2盒
    {
        14,17,11,24,1,5,3,28,15,6,21,10,
        23,19,12,4,26,8,16,7,27,20,13,2,
        41,52,31,37,47,55,30,40,51,45,33,48,
        44,49,39,56,34,53,46,42,50,36,29,32
    };
    const int S[8][4][16]=///8个s盒
    {
        ///s1
        14,4,13,1,2,15,11,8,3,10,6,12,5,9,0,7,
        0,15,7,4,14,2,13,1,10,6,12,11,9,5,3,8,
        4,1,14,8,13,6,2,11,15,12,9,7,3,10,5,0,
        15,12,8,2,4,9,1,7,5,11,3,14,10,0,6,13,
        ///s2
        15,1,8,14,6,11,3,4,9,7,2,13,12,0,5,10,
        3,13,4,7,15,2,8,14,12,0,1,10,6,9,11,5,
        0,14,7,11,10,4,13,1,5,8,12,6,9,3,2,15,
        13,8,10,1,3,15,4,2,11,6,7,12,0,5,14,9,
        ///s3
        10,0,9,14,6,3,15,5,1,13,12,7,11,4,2,8,
        13,7,0,9,3,4,6,10,2,8,5,14,12,11,15,1,
        13,6,4,9,8,15,3,0,11,1,2,12,5,10,14,7,
        1,10,13,0,6,9,8,7,4,15,14,3,11,5,2,12,
        ///s4
        7,13,14,3,0,6,9,10,1,2,8,5,11,12,4,15,
        13,8,11,5,6,15,0,3,4,7,2,12,1,10,14,9,
        10,6,9,0,12,11,7,13,15,1,3,14,5,2,8,4,
        3,15,0,6,10,1,13,8,9,4,5,11,12,7,2,14,
        ///s5
        2,12,4,1,7,10,11,6,8,5,3,15,13,0,14,9,
        14,11,2,12,4,7,13,1,5,0,15,10,3,9,8,6,
        4,2,1,11,10,13,7,8,15,9,12,5,6,3,0,14,
        11,8,12,7,1,14,2,13,6,15,0,9,10,4,5,3,
        ///s6
        12,1,10,15,9,2,6,8,0,13,3,4,14,7,5,11,
        10,15,4,2,7,12,9,5,6,1,13,14,0,11,3,8,
        9,14,15,5,2,8,12,3,7,0,4,10,1,13,11,6,
        4,3,2,12,9,5,15,10,11,14,1,7,6,0,8,13,
        ///s7
        4,11,2,14,15,0,8,13,3,12,9,7,5,10,6,1,
        13,0,11,7,4,9,1,10,14,3,5,12,2,15,8,6,
        1,4,11,13,12,3,7,14,10,15,6,8,0,5,9,2,
        6,11,13,8,1,4,10,7,9,5,0,15,14,2,3,12,
        ///s8
        13,2,8,4,6,15,11,1,10,9,3,14,5,0,12,7,
        1,15,13,8,10,3,7,4,12,5,6,11,0,14,9,2,
        7,11,4,1,9,12,14,2,0,6,10,13,15,3,5,8,
        2,1,14,7,4,10,8,13,15,12,9,0,3,5,6,11,
    };
    const int IP_1[]=///逆初始置换表
    {
        40,8,48,16,56,24,64,32,
        39,7,47,15,55,23,63,31,
        38,6,46,14,54,22,62,30,
        37,5,45,13,53,21,61,29,
        36,4,44,12,52,20,60,28,
        35,3,43,11,51,19,59,27,
        34,2,42,10,50,18,58,26,
        33,1,41,9,49,17,57,25,
    };
    //int smbit[]={0,0,0,1,0,1,0,0, 0,1,0,0,1,1,1,0, 1,1,0,1,0,1,0,0, 1,1,1,0,1,1,1,1,0,0,0,1,0,0,0,1, 0,0,1,0,1,1,0,0, 1,1,1,0,0,1,1,0, 0,0,0,0,1,1,1,0};
    void change(char s[],int flag)///转换为二进制,(验证正确)
    {
        int l=0;
        for(int i=0; i<8; i++)
        {
            int n=s[i],time=0;
            memset(zhongbit,0,sizeof(zhongbit));
            while(n!=0)
            {
                zhongbit[time++]=n%2;
                n/=2;
            }
            if(flag==1)
            {
                for(int j=7; j>=0; j--)
                    smbit[l++]=zhongbit[j];
            }
            else if(flag==2)
            {
                for(int j=7; j>=0; j--)
                    miyaobit[l++]=zhongbit[j];
            }
        }
    }
    void fen()///将64bit的minip分为左右两部分(验证正确)
    {
        for(int i=0; i<32; i++)
        {
            Left[i]=mip[i];
            Right[i]=mip[32+i];
        }
        for(int i=0; i<32; i++)
        {
            LL[cnt][i]=Left[i];
            RR[cnt][i]=Right[i];
        }
        cnt++;
    }
    void ipzhihuan()///第一次循环使用的ip置换,ip置换(验证正确)
    {
        for(int i=0; i<64; i++)
            mip[i]=smbit[IP[i]-1];///下标从0开始
        fen();
        memset(smbit,0,sizeof(smbit));
    }
    void EK()///E盒扩展(验证正确)
    {
        memset(rr,0,sizeof(rr));
        for(int i=0; i<48; i++)
            rr[i]=Right[E[i]-1];///下标从0开始
    }
    void SY(int k)///异或及S盒压缩(验证正确)
    {
        for(int i=0; i<48; i++)
            rr[i]^=zimi[k][i];
        int c=0;
        for(int i=0; i<48; i++)
        {
            if((i+1)%6==0)///进入下个s盒
            {
                int w=(i+1)/6,h,l;///哪个盒,哪一行那一列
                h=rr[i-5]*2+rr[i]*1;
                l=rr[i-4]*8+rr[i-3]*4+rr[i-2]*2+rr[i-1]*1;
                int ans=S[w-1][h][l];
                int tt=0;
                memset(er,0,sizeof(er));
                while(ans)
                {
                    er[tt++]=ans%2;
                    ans/=2;
                }
                for(int j=3; j>=0; j--)
                    SYa[c++]=er[j];
            }
        }
    }
    void PHe()///P盒置换
    {
        for(int i=0; i<32; i++)
            Right[i]=Left[i]^SYa[P[i]-1];
        for(int i=0; i<32; i++)
            Left[i]=mip[32+i];
        for(int i=0; i<32; i++)
        {
            mip[i]=Left[i];
            mip[32+i]=Right[i];
        }
        fen();///分为左右两部分
    }
    void nizhihuan()///逆初始置换表(验证正确)
    {
        for(int i=0; i<64; i++)
            smbit[i]=mip[IP_1[i]-1];
    }
    void PC_1zhi()///置换选择pc1盒及循环移位(验证正确)
    {
        if(times==0)///只有当是第一次的时候才会进行PC置换盒1
        {
            for(int i=0; i<56; i++)
                miyaopc1[i]=miyaobit[PC_1[i]-1];///下标从0开始
        }
        for(int i=0; i<28; i++) ///分为两部分,循环移位
        {
            int ans=i+xz[times];
            if(ans>27)
                ans-=28;
            c[i]=miyaopc1[ans];///得到前半部分
            d[i]=miyaopc1[28+ans];///得到后半部分
        }
        for(int i=0; i<28; i++)
        {
            miyaopc1[i]=c[i];
            miyaopc1[28+i]=d[i];
        }
        times++;
    }
    void PC_2zhi()///16轮置换,pc2盒得到子密钥(验证正确)
    {
        for(int k=0; k<16; k++)
        {
            PC_1zhi();
            for(int i=0; i<48; i++)
                zimi[k][i]=miyaopc1[PC_2[i]-1];///将每一次产生的子密钥存下来
        }
    }
    void zhuan()///将二进制转换成16进制密文
    {
        int ans=0,l=0;
        for(int i=0; i<64; i++)
        {
            if((i+1)%4==0)
            {
                ans=smbit[i-3]*8+smbit[i-2]*4+smbit[i-1]*2+smbit[i];
                if(ans>9)
                    miwen[l++]='A'+ans-10;
                else
                    miwen[l++]=ans+'0';
                ans=0;
            }
        }
    }
    void jiami()///加密算法
    {
        PC_2zhi();///子密钥产生(正确)
        ipzhihuan();///ip置换(正确)
        for(int i=0; i<16; i++)///一样的方式一共循环15次
        {
            EK();///E盒扩展(正确)
            SY(i);///S盒压缩(正确)
            PHe();///p盒置换
        }
        for(int i=0; i<32; i++)
        {
            mip[32+i]=LL[16][i];
            mip[i]=RR[16][i];
        }
        nizhihuan();
        zhuan();
    }
    void jiemi()
    {
        cnt=0;
        memset(LL,0,sizeof(LL));
        memset(RR,0,sizeof(RR));
        ipzhihuan();///ip置换(正确)
        for(int i=0; i<16; i++)///一样的方式一共循环15次
        {
            EK();///E盒扩展(正确)
            SY(15-i);///S盒压缩(正确)
            PHe();///p盒置换
        }
        for(int i=0; i<32; i++)
        {
            mip[32+i]=LL[16][i];
            mip[i]=RR[16][i];
        }
        nizhihuan();
    }
    int main()
    {
        init();
        printf("请输入8位明文:");
        scanf("%s",mingwen);
        change(mingwen,1);
        printf("\n请输入8位密钥:");
        scanf("%s",miyao);
        change(miyao,2);
        jiami();
        printf("加密完成,得到十六进制密文:\n");
        printf("%s\n",miwen);
        jiemi();
        printf("\n解密完成,得到明文:\n\n");
        for(int i=0;i<64;i++)
        {
            if((i+1)%8==0)
            {
                int ans=smbit[i-7]*128+smbit[i-6]*64+smbit[i-5]*32+smbit[i-4]*16+smbit[i-3]*8+smbit[i-2]*4+smbit[i-1]*2+smbit[i];
                printf("%c",ans);
            }
        }
        printf("\n");
        return 0;
    }
    

    运行结果:

    某网站在线加密运行结果:

    结果除了大小写以及尾串不一样,其它相同,加密成功

    展开全文
  • 计算机中快捷键ctrl什么是返回上一步发布时间:2021-06-10 11:33:52来源:亿速云阅读:67作者:小新这篇文章主要为大家展示了“计算机中快捷键ctrl什么是返回上一步”,内容简而易懂,条理清晰,希望能够帮助...
  • C/C++】C语言特性总结

    万次阅读 多人点赞 2019-08-10 16:21:28
    已经有大约半年的时间没有碰C语言了,当时学习的时候记录了很多的笔记,但是都是特别混乱,后悔那个时候,不懂得写博客,这里凭借记忆和零零散散的笔记记录,尝试系统性地复习一下C语言。 之前都是在Windows环境下...
  • 如何从RNN起步,一步一步通俗理解LSTM

    万次阅读 多人点赞 2019-05-06 23:47:54
    如何从RNN起步,一步一步通俗理解LSTM 前言 提到LSTM,之前学过的同学可能最先想到的是ChristopherOlah的博文《理解LSTM网络》,这篇文章确实厉害,网上流传也相当之广,而且当你看过了网上很多关于LSTM的文章...
  • DES一轮加密(C语言

    千次阅读 2015-01-23 20:58:59
    输出每一步 的结果。 已知R0=00000000 11111111 00000110 10000011, 子密钥K1=00111101 10001111 11001101 00110111 00111111 00000110 实现E盒扩展、与轮密钥异或、S盒代换(8个一样的S盒,都为S1盒)和P盒置换...
  • } //这个函数能够 把 字符串数组中的 pos1到pos2位置的 数字 转化成int 类型 方便减乘除的计算 int f(int pos1,int pos2) { int sum=0; for(int i=pos1 ; i; i++) { if(s[i][i]>='0') sum=sum*10+(s[i]-'0');...
  • 程序员是如何一步一步被诈骗的?

    万次阅读 多人点赞 2020-05-01 09:50:59
    她说会有一个她们公司的技术部门的人专门来指导我销户,只需要几分钟的时间,让我一下客服的 QQ,说是全程指导销户(我居然就了扣扣,也是傻得可以)。加上所谓的工作人员的 QQ 以后,我进开始一步一步陷入他们...
  • C语言程序从编写到运行历经的几个阶段一 ...但实际上C语言程序从编写到运行,这期间的经历并不是这么简单,那么现在小编就带领大家探索,这期间具体有几个步骤?一 过程简介从上图可知从C源码到可执行程序,我们会...
  • AES/密解密算法 (C)

    2010-12-25 15:53:30
    能够对一个明文分组进行解密,并能够显示每一步的步骤
  • C语言基础第五课写在前面:一、字符输出函数putchar二、字符输入函数getchar三、字符型数据与整型数据可以进行混合运算四、强制类型转换五、顺序结构程序设计(按顺序一步一步来) 写在前面: 简单谈谈写程序吧,写...
  • 最近简单的学习了一些常用的数据处理算法,在这里我简单的做一个系列笔记和大家一起学习国内外一些常见的数据处理算法,之前也做了java的实现,但是由于考虑到性能问题最终还是采用的C语言实现(文章末尾附带VC6.0 ...
  • c简单实现Feistel解密

    千次阅读 2019-03-22 17:09:06
    而解密时最后一步也是如此。 #include #include #include #include int main() { int i=0,j=0; char key[50]; //key最大位数50位,但是还是只取16位 char plaintext[100]; //明文 char pGroup[10][10]; ...
  • c语言

    2019-09-29 22:43:14
    原创 C语言在开发中的应用博文汇总贴 置顶 ...
  • 为小学的广大学子写一个,减乘除法做题系统,思路简单清晰,欢迎品尝 文章目录为小学的广大学子写一个,减乘除法做题系统,思路简单清晰,欢迎品尝代码应解决的问题:一、代码如下二、对代码进行测试总结: ...
  • Linux网络编程基础和一步一步

    千次阅读 多人点赞 2014-01-14 12:24:25
    ·Linux网络编程 基础(一) ·Linux网络编程 基础(二) ·Linux网络编程 基础(三) ...• Linux网络编程一步一步学-简单客户端编写  • Linux网络编程一步一步学-绑定IP和端口 • Linux网络编程一步
  • 初学者无需纠结应该学种语言。当你一种入门了,其他语言无非语法不同,要实现的业务功能是一样的。所以,一步一步认真学,一定要踏实,多练多思考。 Java特点 • Java是跨平台的 • Java是简单的 • Java是安全的...
  • 一步一步写ARM汇编(一)

    万次阅读 多人点赞 2018-06-05 21:18:48
    本博文开始学习一步一步写ARM汇编程序。 一、重要概念理解1. 立即数1)把数据转换成二进制形式,从低到高写成 4位1组的形式,最高位一组不够4位的前面补02)数1的个数,如果大于8个【可能也是立即数,取反】不是...
  • 作者:C 更新时间:2012-8-20  平衡二叉树(Balanced Binary Tree)是二叉查找树的一个进化体,也是第一个引入平衡概念的二叉树。1962年,G.M. Adelson-Velsky 和 E.M. Landis发明了这棵树,所以它又叫AVL树...
  • gdb 一步一步调试程序

    千次阅读 2015-01-23 11:41:51
    -g 之后可执行文件中会记录编译路径这些信息,list 时 gdb 去读编译路径下的源码,如果编译路径下的源码删除或者没有源码,那么 list 是显示不了代码的。 (gdb) info source Current source file is main.c...
  • 高精度加法利用竖式计算的原理,通过处理进位,来计算高精度数据的和。 核心代码(处理进位):...void add(int a[],int b[]) //a,b,c都为数组,分别存储被数、数、结果 { int i=1,x=0; //x是进位 while ((i数组
  • ThinkPad E430c加装内存和SSD固态硬盘硬件准备1、光驱位硬盘托架2、SSD固态硬盘3、内存条操作步骤:1、拆下电池2、松动后盖的三个固定螺丝注意:螺丝是拿不下来的! 这种设计很好,防止螺丝丢失或用错位置!3、打开...
  • 一个类型里会出现很多运算符,他们也像普通的表达式一样,有优先级,其优先级和运算优先级一样,所以我总结了一下其原则:从变量名处起,根据运算符优先级结合,一步一步分析.下面让我们先从简单的类型开始慢慢分析吧: int...
  • 下面是基于一个简单例子,一步一步注释翻译 A look at a very simple neural network inTensorFlow 让我们来看一个基于tensorFlow创建一个非常简单的神经网络,实现线性回归的功能。   This is an ...
  • C语言结构体(struct)常见使用方法

    万次阅读 多人点赞 2014-04-14 01:51:57
    注意:盗版是不会得到修正和更新的! 今天复习一下struct,顺便挖掘一下以前没注意的小...(因为C++和C有共通之处,但是在结构体上的某些机制又有所不同,所以后边提了一下,不喜欢可以略过) 结构体定义: ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 304,194
精华内容 121,677
关键字:

哪一步加c