devbox
IT小白
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

关于GVF的一些理解_distance transforms of sampled functions

作者：IT小白 | 2024-04-11 19:52:30

赞

踩

distance transforms of sampled functions

Distance Transforms of Sampled Functions 

基于采样距离变换算法

https://blog.csdn.net/frozenspring/article/details/78357858

距离变换及其数学推导

https://blog.csdn.net/abcjennifer/article/details/7617883

距离变换

https://blog.csdn.net/lingyunxianhe/article/details/73555520

二值图距离变换

https://blog.csdn.net/qq_34784753/article/details/68951918

【OpenCV图像处理】十三、图像的距离变换

Matlab部分知识了解

1--浅析image,imagesc,imshow的用法

https://blog.csdn.net/zhupananhui/article/details/16340345

magesc: 在figure2中我们用imagesc来显示图像与figure1相比能较好的显示出来，同样我们也得搞明白调用imagesc时矩阵的数和colormap中索引的对应关系，与image不同的是imagesc采用的不是直接映射而是线性映射，至于什么是线性映射，我粗略的说一下，比如把区间A = [0,a]映射到区间B = [0,b]我们对A中的元素做A/a*b就可以了，矩阵的数到colormap索引的线性映射大概就是这样，Matlab会自动获取矩阵中数的最小值和最大值，并把区间[Cmin,Cmax]映射到colormap的[最小索引，最大索引]比如[1,64],然后再根据这个对应关系把图像显示出来，具体的算法细节是Matlab确定的，当然也可以自己指定显示范围，比如一副索引图像I范围为[27,218]，而我只想显示[1 64 ],使用命令imagesc(I,[1 64])就可以了，如果你把上面程序中的imagesc(I)换成imagesc(I,[1,64])，那么figure2中的效果就和figure1中一样了，因为只是把[1,64]这个范围映射到色图，超过的都被认为是64。关于映射，我截图Matlab中imagesc的help页给大家看看，这里要自己慢慢体会哦，使用imagesc(I)这种线性映射就可以用到整个色图从而将图像较好的显示出来，这就是figure2中的显示效果比figure1中好的原因。

 

  1. A = [4 7 2 9 8;3 9 1 4 3;1 5 9 6 4;8 3 7 1 0]
  2. [i,j] = ind2sub(A,5)
  3. %3×3矩阵的第 5个元素的全下标；单下标转变为多下标
  4. ind=sub2ind(size(A),[1,2;3,4],[1,1;2,2])
  5. %矩阵第三行、第三列元素的序号
  6.  
  7. A =
  8.  
  9.      4     7     2     9     8
  10.      3     9     1     4     3
  11.      1     5     9     6     4
  12.      8     3     7     1     0
  13.  
  14.  
  15. i =
  16.  
  17.      1
  18.  
  19.  
  20. j =
  21.  
  22.      2
  23.  
  24.  
  25. ind =
  26.  
  27.      1     2
  28.      7     8
  1. ind=sub2ind(siz,I,J)：siz表示要转换的矩阵的行列数，I是要转换矩阵的行标，J是要转换矩阵的列标。I,J的行列数必须相同。
  2. 可以看出，在矩阵A中，下标（1,1）的索引值为1，下标（2,1）的索引值为2，下标（3,2）的索引值为7，下标（4,2）的索引值为8

2--matlab的double和single类型 

 代码部分

  1. % load image
  2. im = imread(fullfile(vl_root, 'data', 'roofs1.jpg')) ;
  3. im = im(100:200,100:200,:) ;
  4. imSize = [size(im,1) size(im,2)] ;
  5.  
  6. % creates an edge map
  7. edges = zeros(imSize) + inf;
  8. edges(edge(rgb2gray(im), 'canny')) = 0 ;
  9.  
  10. % compute distance transform
  11. [distanceTransform, neighbors] = vl_imdisttf(single(edges)) ;
  12.  
  13. % plot
  14. [u,v] = meshgrid(1:imSize(2),1:imSize(1)) ;
  15. [v_,u_] = ind2sub(imSize, neighbors) ;
  16.  
  17. % avoid cluttering the plot too much
  18. u = u(1:3:end,1:3:end) ;
  19. v = v(1:3:end,1:3:end) ;
  20. u_ = u_(1:3:end,1:3:end) ;
  21. v_ = v_(1:3:end,1:3:end) ;
  22.  
  23. figure(1) ; clf ; imagesc(im) ; axis off image ;
  24. figure(2) ; clf ; imagesc(edges) ; axis off image ;
  25. figure(3) ; clf ; imshow(edges) ; axis off image ;
  26. hold on ; h = quiver(u,v,u_-u,v_-v,0) ; colormap gray ;
  27. figure(4) ; clf ; imagesc(sqrt(distanceTransform)) ; axis off image ;

vl_imdisttf函数：

输入为：二值图canny结果

输出为：[distanceTransform, neighbors]

其中neighbors的值表示每个像素点与图片中某个梯度点相邻，值是该点的位置

输入图像大小为101*101  （10201）

[v_,u_] = ind2sub(imSize, neighbors) ;

v_ 和 u_ 构成了neighbors对应的的坐标矩阵

比如104对用（2，3）----AC是怎么构成的看后面的解析

其中distanceTransform的值表示每个像素点与图片中相邻梯度点的距离值

定义：二值图像，任意点到最近背景点的距离，一般为非零点到最近零点的距离。

(1,1)最近的最近背景点是索引为：104

转换为下表为:(2,3)

(1,1)到(2,3)的距离为：5=(2-1)^2+(3-1)^2

MATLAB中数据是按列的方式存储的，索引是一维向量。，而非多行1列的二维向量

  1. 在命令窗口输入：
  2. >> A=[4 7 2 9 8;3 9 1 4 3;1 5 9 6 4;8 3 7 1 0]
  3. A =
  4.      4     7     2     9     8
  5.      3     9     1     4     3
  6.      1     5     9     6     4
  7.      8     3     7     1     0
  8. 则A中每个元素对应的索引如下：
  9.      1     5     9    13    17
  10.      2     6    10    14    18
  11.      3     7    11    15    19
  12.      4     8    12    16    20

我们可以通过：AC=reshape(1:10201,101,101)

生成顺序矩阵，按列的方式存储索引(即matlab的方式)

或者按行生成顺序矩阵

AC=reshape(1:10201,101,101)

AB = AC’

[distanceTransform, neighbors] = vl_imdisttf(single(edges)) ;

 

  1. /** @internal
  2.  ** @file    vl_imdisttf.c
  3.  ** @brief   vl_imdisttf - MEX definition
  4.  ** @author  Andrea Vedaldi
  5.  **/
  6.  
  7. /*
  8. Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
  9. All rights reserved.
  11. This file is part of the VLFeat library and is made available under
  12. the terms of the BSD license (see the COPYING file).
  13. */
  14.  
  15. #include <mexutils.h>
  16.  
  17. #include <vl/generic.h>
  18. #include <vl/mathop.h>
  19. #include <vl/imopv.h>
  20.  
  21. void
  22. mexFunction(int nout, mxArray *out[],
  23.             int nin, const mxArray *in[])
  24. {
  25.   vl_size M, N ;
  26.   enum {IN_I = 0, IN_PARAM, IN_END} ;
  27.   enum {OUT_DT = 0, OUT_INDEXES} ;
  28.   vl_uindex * indexes = NULL ;
  29.   mxClassID classId ;
  30.   double const defaultParam [] = {1.0, 0.0, 1.0, 0.0} ;
  31.   double const * param = defaultParam ;
  32.  
  33.   /* -----------------------------------------------------------------
  34.    *                                               Check the arguments
  35.    * -------------------------------------------------------------- */
  36.  
  37.   if (nin < 1) {
  38.     vlmxError(vlmxErrNotEnoughInputArguments, NULL) ;
  39.   }
  40.   if (nin > 2) {
  41.     vlmxError(vlmxErrTooManyInputArguments, NULL) ;
  42.   }
  43.   if (nout > 2) {
  44.     vlmxError(vlmxErrTooManyOutputArguments, NULL) ;
  45.   }
  46.   classId = mxGetClassID(IN(I)) ;
  47.   if (! vlmxIsMatrix(IN(I), -1, -1) ||
  48.       (classId != mxSINGLE_CLASS && classId != mxDOUBLE_CLASS)) {
  49.     vlmxError(vlmxErrInvalidArgument,
  50.              "I is not a SINGLE or DOUBLE matrix.") ;
  51.   }
  52.   if (nin == 2) {
  53.     if (! vlmxIsPlainVector(IN(PARAM), 4)) {
  54.       vlmxError(vlmxErrInvalidArgument,
  55.                "PARAM is not a 4-dimensional vector.") ;
  56.     }
  57.     param = mxGetPr (IN(PARAM)) ;
  58.     if (param[0] < 0.0 ||
  59.         param[2] < 0.0) {
  60.       vlmxError(vlmxErrInvalidArgument,
  61.                 "Either PARAM[0] or PARAM[2] is negative.") ;
  62.     }
  63.   }
  64.  
  65.   M = mxGetM (IN(I)) ;
  66.   N = mxGetN (IN(I)) ;
  67.  
  68.   OUT(DT) = mxCreateNumericMatrix (M, N, classId, mxREAL) ;
  69.   if (nout > 1) {
  70.     vl_uindex i ;
  71.     OUT(INDEXES) = mxCreateDoubleMatrix (M, N, mxREAL) ;
  72.     indexes = mxMalloc(sizeof(vl_uindex) * M * N) ;
  73.     for (i = 0 ; i < M * N ; ++i) indexes[i] = i + 1 ;
  74.   }
  75.  
  76.   /* -----------------------------------------------------------------
  77.    *                                                        Do the job
  78.    * -------------------------------------------------------------- */
  79.  
  80.   switch (classId) {
  81.     case mxSINGLE_CLASS:
  82.       vl_image_distance_transform_f((float const*)mxGetData(IN(I)),
  83.                                     M, N,
  84.                                     1, M,
  85.                                     (float*)mxGetPr(OUT(DT)),
  86.                                     indexes,
  87.                                     param[2],
  88.                                     param[3]) ;
  89.  
  90.       vl_image_distance_transform_f((float*)mxGetPr(OUT(DT)),
  91.                                     N, M,
  92.                                     M, 1,
  93.                                     (float*)mxGetPr(OUT(DT)),
  94.                                     indexes,
  95.                                     param[0],
  96.                                     param[1]) ;
  97.       break ;
  98.  
  99.     case mxDOUBLE_CLASS:
  100.       vl_image_distance_transform_d((double const*)mxGetData(IN(I)),
  101.                                     M, N,
  102.                                     1, M,
  103.                                     (double*)mxGetPr(OUT(DT)),
  104.                                     indexes,
  105.                                     param[2],
  106.                                     param[3]) ;
  107.  
  108.       vl_image_distance_transform_d((double*)mxGetPr(OUT(DT)),
  109.                                     N, M,
  110.                                     M, 1,
  111.                                     (double*)mxGetPr(OUT(DT)),
  112.                                     indexes,
  113.                                     param[0],
  114.                                     param[1]) ;
  115.       break;
  116.  
  117.     default:
  118.       abort() ;
  119.   }
  120.  
  121.   if (indexes) {
  122.     vl_uindex i ;
  123.     double * pt = mxGetPr(OUT(INDEXES)) ;
  124.     for (i = 0 ; i < M * N ; ++i) pt[i] = indexes[i] ;
  125.     mxFree(indexes) ;
  126.   }
  127. }

 

  1. VL_EXPORT void
  2. vl_image_distance_transform_f (float const * image,
  3.                                vl_size numColumns,
  4.                                vl_size numRows,
  5.                                vl_size columnStride,
  6.                                vl_size rowStride,
  7.                                float * distanceTransform,
  8.                                vl_uindex * indexes,
  9.                                float coeff,
  10.                                float offset) ;
  1. VL_EXPORT void
  2. VL_XCAT(vl_image_distance_transform_,SFX)
  3. (T const * image,
  4.  vl_size numColumns,
  5.  vl_size numRows,
  6.  vl_size columnStride,
  7.  vl_size rowStride,
  8.  T * distanceTransform,
  9.  vl_uindex * indexes,
  10.  T coeff,
  11.  T offset)
  12. {
  13.   /* Each image pixel corresponds to a parabola. The algorithm scans
  14.    such parabolas from left to right, keeping track of which
  15.    parabolas belong to the lower envelope and in which interval. There are
  16.    NUM active parabolas, FROM stores the beginning of the interval
  17.    for which a certain parabola is part of the envoelope, and WHICH store
  18.    the index of the parabola (that is, the pixel x from which the parabola
  19.    originated).
  20. 每个图像像素对应一个抛物线。算法扫描这样的抛物线从左到右，
  21. 跟踪哪个抛物线属于下包层和哪个区间。 有NUM个有效抛物线，
  22. FROM存储某个抛物线是envoelope的一部分的间隔的开始，
  23. WHICH存储抛物线的指数（即抛物线起源的像素x）。
  24.    */
  25.   vl_uindex x, y ;
  26.   T * from = vl_malloc (sizeof(T) * (numColumns + 1)) ;
  27.   T * base = vl_malloc (sizeof(T) * numColumns) ;
  28.   vl_uindex * baseIndexes = vl_malloc (sizeof(vl_uindex) * numColumns) ;
  29.   vl_uindex * which = vl_malloc (sizeof(vl_uindex) * numColumns) ;
  30.   vl_uindex num = 0 ;
  31.  
  32.   for (y = 0 ; y < numRows ; ++y) {
  33.     num = 0 ;
  34.     for (x = 0 ; x < numColumns ; ++x) {
  35.       T r = image[x  * columnStride + y * rowStride] ;
  36.       T x2 = x * x ;
  37. #if (FLT == VL_TYPE_FLOAT)
  38.       T from_ = - VL_INFINITY_F ;
  39. #else
  40.       T from_ = - VL_INFINITY_D ;
  41. #endif
  42.  
  43.       /*
  44.        Add next parabola (there are NUM so far). The algorithm finds
  45.        intersection INTERS with the previously added parabola. If
  46.        the intersection is on the right of the "starting point" of
  47.        this parabola, then the previous parabola is kept, and the
  48.        new one is added to its right. Otherwise the new parabola
  49.        "eats" the old one, which gets deleted and the check is
  50.        repeated with the parabola added before the deleted one.
  51. 添加下一个抛物线（到目前为止有NUM个）。 算法找到交叉口INTERS与先前添加的抛物线。 
  52. 如果交叉点位于“起点”的右侧这个抛物线，然后保持先前的抛物线，并且右边添加了新的。 
  53. 否则新的抛物线“吃掉”旧的，将被删除，并在删除之前添加抛物线重复检查。
  54.        */
  55.  
  56.       while (num >= 1) {
  57.         vl_uindex x_ = which[num - 1] ;
  58.         T x2_ = x_ * x_ ;
  59.         T r_ = image[x_ * columnStride + y * rowStride] ;
  60.         T inters ;
  61.         if (r == r_) {
  62.           /* handles the case r = r_ = \pm inf */
  63.           inters = (x + x_) / 2.0 + offset ;
  64.         }
  65. #if (FLT == VL_TYPE_FLOAT)
  66.         else if (coeff > VL_EPSILON_F)
  67. #else
  68.         else if (coeff > VL_EPSILON_D)
  69. #endif
  70.         {
  71.           inters = ((r - r_) + coeff * (x2 - x2_)) / (x - x_) / (2*coeff) + offset ;
  72.         } else {
  73.           /* If coeff is very small, the parabolas are flat (= lines).
  74.            In this case the previous parabola should be deleted if the current
  75.            pixel has lower score 
  76.           如果coeff非常小，抛物线是平的（=线）。 在这种情况下，如果当前像素的得分较低，则应删除 
  77.           先前的抛物线
  78.           */
  79. #if (FLT == VL_TYPE_FLOAT)
  80.           inters = (r < r_) ? - VL_INFINITY_F : VL_INFINITY_F ;
  81. #else
  82.           inters = (r < r_) ? - VL_INFINITY_D : VL_INFINITY_D ;
  83. #endif
  84.         }
  85.         if (inters <= from [num - 1]) {
  86.           /* delete a previous parabola */
  87.           -- num ;
  88.         } else {
  89.           /* accept intersection */
  90.           from_ = inters ;
  91.           break ;
  92.         }
  93.       }
  94.  
  95.       /* add a new parabola */
  96.       which[num] = x ;
  97.       from[num] = from_ ;
  98.       base[num] = r ;
  99.       if (indexes) baseIndexes[num] = indexes[x  * columnStride + y * rowStride] ;
  100.       num ++ ;
  101.     } /* next column */
  102.  
  103. #if (FLT == VL_TYPE_FLOAT)
  104.     from[num] = VL_INFINITY_F ;
  105. #else
  106.     from[num] = VL_INFINITY_D ;
  107. #endif
  108.  
  109.     /* fill in */
  110.     num = 0 ;
  111.     for (x = 0 ; x < numColumns ; ++x) {
  112.       double delta ;
  113.       while (x >= from[num + 1]) ++ num ;
  114.       delta = (double) x - (double) which[num] - offset ;
  115.       distanceTransform[x  * columnStride + y * rowStride]
  116.       = base[num] + coeff * delta * delta ;
  117.       if (indexes) {
  118.         indexes[x  * columnStride + y * rowStride]
  119.         = baseIndexes[num] ;
  120.       }
  121.     }
  122.   } /* next row */
  123.  
  124.   vl_free (from) ;
  125.   vl_free (which) ;
  126.   vl_free (base) ;
  127.   vl_free (baseIndexes) ;
  128. }
  1. /** @internal @brief IEEE single precision infinity constant */
  2. static union { vl_uint32 raw ; float value ; }
  3.   const vl_infinity_f =
  4.     { 0x7F800000UL } ;
  5.  
  6. /** @internal @brief IEEE double precision infinity constant */
  7. static union { vl_uint64 raw ; double value ; }
  8.   const vl_infinity_d =
  9. #ifdef VL_COMPILER_MSC
  10.     { 0x7FF0000000000000ui64 } ;
  11. #else
  12.     { 0x7FF0000000000000ULL } ;
  13. #endif

 

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/IT小白/article/detail/407060
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号