【算法随记二】线卷积积分及其在图像增强和特效方面的应用(一)。

LIC(LineIntegralConvolution)isawell-knowntexturesynthesistechniqueproposedbyCabralandLeedom[33]atLawrenceLivermoreNationalLaboratoryinACMSigGraph93.Itemploysalow-passfiltertoconvolveaninputnoisetextur...

【算法随记二】线卷积积分及其在图像增强和特效方面的应用(一)。

  LIC (Line Integral Convolution) is a well-known texture synthesis technique proposed by Cabral and Leedom[33]atLawrence Livermore National LaboratoryinACM SigGraph 93. It employs a low-pass filter to convolve an input noise texture along pixel-centered symmetrically bi-directional streamlines to exploit spatial correlation in the flow direction. LIC provides a global dense representation of the flow, emulating what happens when a rectangular area of massless fine sand is blown by strong wind (Figure 1).

以上内容和图片摘自http://www.zhanpingliu.org/Research/FlowVis/LIC/LIC.htm,这是一个中国人的网站,并且还是一个很老的网站。

  关于LIC的实现代码,网络上流传的最为广泛的也是这里的代码,详见:http://www.zhanpingliu.org/Research/FlowVis/LIC/LIC_Source.htm

按照个人的理解,LIC算法他就是把一幅矢量场数据用图像的方式可视化出来,那么对于某一个固定位置的点,他其实是只有当前点的矢量值的,一个带有方向信息的Vector(一般都是归一化后的数据,即矢量的长度为1)。要把该点转换为一个像素值,那么我们需要首先有个基点,对于全局来说就是一幅图,通常情况下,我们会产生一幅和矢量数据大小一样的白噪声图来作为基点图,之后,我们采用卷积的方法,沿着当前点的矢量方向前进一定的距离D,得到新的坐标位置,记录下这个位置在基点图中的像素值,并累加,之后,这个新位置也有他对应的矢量方向,沿着这个矢量方向继续前进,并执行相同的累加操作,直到前进了指定的距离后,再计算累加后的平均值最为可视化后的像素值。

  整个流程其实看起来就是沿着某一条线,对线上相关离散点进行卷积,所以严格的说可以叫做折线卷积。同时,为了保证对称性,我们会在卷积起点时也会沿着矢量相反的方向进行卷积,很明显,这个反响卷积的路线和正向卷积的一半来说不会时对称的。

  那么谈到代码,我把上面zhanpingliu的方案整理成了一个可运行的C 方案,并未做任何的优化,对一幅512*512的灰度图做测试,由矢量数据生成图像的耗时约为145ms,还是有点慢的,那么我们看下他的代码怎么写的:

  1 ///  flow imaging (visualization) through Line Integral Convolution  ///  2 void FlowImagingLIC(int  Width, int  Height, float*  pVectr, unsigned char*  pNoise, unsigned char*  pImage,  3  float*  p_LUT0, float*  p_LUT1, floatkrnlen)  4 {  5  int  vec_id;///ID in the VECtor buffer (for the input flow field)  6  int  advDir;///ADVection DIRection (0: positive;  1: negative)  7  int  advcts;///number of ADVeCTion stepS per direction (a step counter)  8  int  ADVCTS = int(krnlen * 3); ///MAXIMUM number of advection steps per direction to break dead loops   9  10  float vctr_x;///x-component  of the VeCToR at the forefront point 11  float vctr_y;///y-component  of the VeCToR at the forefront point 12  float clp0_x;///x-coordinate of CLiP point 0 (current) 13  float clp0_y;///y-coordinate of CLiP point 0 (current) 14  float clp1_x;///x-coordinate of CLiP point 1 (next) 15  float clp1_y;///y-coordinate of CLiP point 1 (next) 16  float samp_x;///x-coordinate of the SAMPle in the current pixel 17  float samp_y;///y-coordinate of the SAMPle in the current pixel 18  float tmpLen;///TeMPorary LENgth of a trial clipped-segment 19  float segLen;///SEGmentLENgth 20  float curLen;///CURrentLENgth of the streamline 21  float prvLen;///PReVious  LENgth of the streamline   22  float W_ACUM;///ACcuMulated Weight from the seed to the current streamline forefront 23  float texVal;///TEXture VALue 24  float smpWgt;///WeiGhT of the current SaMPle 25  float t_acum[2];  ///two ACcUMulated composite Textures for the two directions, perspectively 26  float w_acum[2];  ///two ACcUMulated Weighting valuesfor the two directions, perspectively 27  float* wgtLUT = NULL; ///WeiGhT Look Up Table pointing to the target filter LUT 28  float len2ID = (DISCRETE_FILTER_SIZE - 1) / krnlen; ///map a curve LENgth TO an ID in the LUT 29  ///for each pixel in the 2D output LIC image/// 30  for (int j = 0; j < Height; j  ) 31  for (int i = 0; i < Width; i  ) 32  { 33///init the composite texture accumulators and the weight accumulators/// 34t_acum[0] = t_acum[1] = w_acum[0] = w_acum[1] = 0.0f; 35///for either advection direction/// 36for (advDir = 0; advDir < 2; advDir  ) 37   { 38 ///init the step counter, curve-length measurer, and streamline seed/// 39 advcts = 0; 40 curLen = 0.0f; 41 clp0_x = i0.5f; 42 clp0_y = j0.5f; 43  44 ///access the target filter LUT/// 45 wgtLUT = (advDir == 0) ? p_LUT0 : p_LUT1; 46  47 ///until the streamline is advected long enough or a tightly  spiralling center / focus is encountered/// 48 while (curLen < krnlen && advcts < ADVCTS) 49 { 50  ///access the vector at the sample/// 51  vec_id = (int(clp0_y) * Widthint(clp0_x)) << 1; 52  vctr_x = pVectr[vec_id]; 53
源文地址:https://www.guoxiongfei.cn/cntech/16242.html