有沒有大神能告知lic（線積分卷積）具體的代碼解釋以及如何使用，或者推薦一篇詳細解釋lic代碼**具體演算法流程**的論文或網址，看過好多論文都沒有提及具體的演算法流程，和lic代碼對不太上，求大神們指點迷津。下附帶lic代碼
#include <math.h>
#include <stdio.h>
#include<malloc.h>
#include<iostream>
#include<stdlib.h>
#define SQUARE_FLOW_FIELD_SZ 400
#define DISCRETE_FILTER_SIZE 2048 //卷積核的對應項，如果太小的話有些采樣點就沒有被卷積進去，影響影像效果
#define LOWPASS_FILTR_LENGTH 32.00000f
#define LINE_SQUARE_CLIP_MAX 100000.0f //只要大于等于0就行，等于0時會出現黑色線條
#define VECTOR_COMPONENT_MIN 0.050000f
/*----------------------------------------------------------------------------------------------*/
void SyntheszSaddle(int n_xres, int n_yres, float* pVector);
void NormalizVectrs(int n_xres, int n_yres, float* pVector);
void GenBoxFiltrLUT(int LUTsiz, float* p_LUT0, float* p_LUT1);
void MakeWhiteNoise(int n_xres, int n_yres, unsigned char* pNoise);
void FlowImagingLIC(int n_xres, int n_yres, float* pVector, unsigned char* pNoise,
unsigned char* pImage, float* p_LUT0, float* p_LUT1, float krnlen);
void WriteImage2PPM(int n_xres, int n_yres, unsigned char* pImage, char* f_name);
/*----------------------------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------------------------*/
void main()
{
int n_xres = SQUARE_FLOW_FIELD_SZ;
int n_yres = SQUARE_FLOW_FIELD_SZ;
float* pVectr = (float* )malloc( sizeof(float ) * n_xres * n_yres * 2 );
float* p_LUT0 = (float* )malloc( sizeof(float ) * DISCRETE_FILTER_SIZE);
float* p_LUT1 = (float* )malloc( sizeof(float ) * DISCRETE_FILTER_SIZE);
unsigned char* pNoise = (unsigned char*)malloc( sizeof(unsigned char) * n_xres * n_yres );
unsigned char* pImage = (unsigned char*)malloc( sizeof(unsigned char) * n_xres * n_yres );
SyntheszSaddle(n_xres, n_yres, pVectr);
// NormalizVectrs(n_xres, n_yres, pVectr);/*資料不規整一下的話，速度越小的地方噪聲越強
MakeWhiteNoise(n_xres, n_yres, pNoise);
GenBoxFiltrLUT(DISCRETE_FILTER_SIZE, p_LUT0, p_LUT1);
FlowImagingLIC(n_xres, n_yres, pVectr, pNoise, pImage, p_LUT0, p_LUT1, LOWPASS_FILTR_LENGTH);
WriteImage2PPM(n_xres, n_yres, pImage, "LIC.ppm");
free(pVectr); pVectr = NULL;
free(p_LUT0); p_LUT0 = NULL;
free(p_LUT1); p_LUT1 = NULL;
free(pNoise); pNoise = NULL;
free(pImage); pImage = NULL;
}
/*----------------------------------------------------------------------------------------------*/
/*合成資料----------------------------------------------------*/
void SyntheszSaddle(int _row, int _col, float* pVectr)
{
int row = 400, col = 400, index = 0;
float vec_x = 0.0f, vec_y = 0.0f, vcMag = 0.0f, scale = 0.0f;
for(int i = 0; i < row; i++)
{
for(int j = 0; j < col; j++)
{
// 資料生成
index = i * col + j;// 中心矢量場
// index = (row - 1 - i) * col + j;// 馬鞍矢量場
vec_x = -(float)i / row + 0.5f;
vec_y = (float)j / col - 0.5f;
// 資料規整
vcMag = float(sqrt(vec_x * vec_x + vec_y * vec_y));
scale = (vcMag < 0.001f) ? 0.0f : 1.0f / vcMag;
vec_x *= scale;
vec_y *= scale;
pVectr[2*index] = vec_x;
pVectr[2*index + 1] = vec_y;
}
}
}
/*------------------------------------------------------------*/
/*資料規整------------------------------------------------------------------------------------------------------------*/
void NormalizVectrs(int n_xres, int n_yres, float* pVectr)
{
for(int j = 0; j < n_yres; j++)
{
for(int i = 0; i< n_xres; i++)
{
int index = (j * n_xres + i) << 1;
float vcMag = float( sqrt( double(pVectr[index] * pVectr[index] + pVectr[index + 1] * pVectr[index + 1]) ) );
float scale = (vcMag == 0.0f) ? 0.0f : 1.0f / vcMag;
pVectr[index ] = scale;//矢量分量的放大倍數，越大，速度越大，程式運行速度會減慢，但效果會更清晰
pVectr[index + 1] *= scale;
}//for
}//for
}
/*--------------------------------------------------------------------------------------------------------------------*/
/*生成白噪聲----------------------------------------------------*/
void MakeWhiteNoise(int n_xres, int n_yres, unsigned char* pNoise)
{
for(int j = 0; j < n_yres; j++)
{
for(int i = 0; i< n_xres; i++)
{
int r = rand();
r = ( (r & 0xff) + ( (r & 0xff00) >> 8 ) ) & 0xff;//
pNoise[j * n_xres + i] = (unsigned char) r;
}
}
}
/*--------------------------------------------------------------*/
/*卷積核LUTsiz = DISCRETE_FILTER_SIZE ----------------------*/
void GenBoxFiltrLUT(int LUTsiz, float* p_LUT0, float* p_LUT1)
{
for(int i = 0; i < LUTsiz; i++)
{
p_LUT0[i] = p_LUT1[i] = (float)i;
}
}
/*----------------------------------------------------------*/
/*寫檔案----------------------------------------------------------------------*/
void WriteImage2PPM(int n_xres, int n_yres, unsigned char* pImage, char* f_name)
{
FILE* o_file;
if( (o_file = fopen(f_name, "w") ) == NULL )
{
printf("Can't open output file\n");
return;
}
fprintf(o_file, "P6\n%d %d\n255\n", n_xres, n_yres);
for(int j = 0; j < n_yres; j++)
{
for(int i = 0; i< n_xres; i++)
{
unsigned char unchar = pImage[j * n_xres + i];
fprintf(o_file, "%c%c%c", unchar, unchar, unchar);
}
}
fclose(o_file);
o_file = NULL;
}
/*-----------------------------------------------------------------------------*/
/*計算像素點------------------------------------------------------------------------------------------------*/
void FlowImagingLIC(int n_xres, int n_yres, float* pVectr, unsigned char* pNoise, unsigned char* pImage,
float* p_LUT0, float* p_LUT1, float krnlen)
{
int vec_id = 0;
int advDir = 0; //方向
int advcts = 0; //追蹤的步數
int ADVCTS = int(krnlen * 3);//追蹤的最大步數
float vctr_x = 0; //速度x分量
float vctr_y = 0; //速度y分量
float clp0_x = 0; //當前點x坐標
float clp0_y = 0; //當前點y坐標
float clp1_x = 0; //下一點x坐標
float clp1_y = 0; //下一點y坐標
float samp_x = 0; //采樣點x坐標
float samp_y = 0; //采樣點y坐標
float tmpLen = 0; //臨時長度
float segLen = 0; //每段長度
float curLen = 0; //當前的流線長度
float prvLen = 0; //上一條流線長度
float W_ACUM = 0; //計算權重之和
float texVal = 0; //紋理的灰度值
float smpWgt = 0; //當前采樣點的權重值
float t_acum[2]; //輸入紋理對應流線上的灰度值之和
float w_acum[2]; //權重之和
float* wgtLUT = NULL;//權重查找表
float len2ID = (DISCRETE_FILTER_SIZE - 1) / krnlen;//將曲線長度映射到卷積核的某一項
//krnlen等于LOWPASS_FILTR_LENGTH
for(int j = 0; j < n_yres; j++)
{
for(int i = 0; i < n_xres; i++)
{
t_acum[0] = t_acum[1] = w_acum[0] = w_acum[1] = 0.0f;
//追蹤方向，正向流線與反向流線
for(advDir = 0; advDir < 2; advDir++)
{
//**初始化追蹤步數、追蹤長度、種子點位置
advcts = 0;
curLen = 0.0f;
clp0_x = i + 0.5f;
clp0_y = j + 0.5f;
//**獲取相應的卷積核
wgtLUT = (advDir == 0) ? p_LUT0 : p_LUT1;
//**回圈終止條件：流線追蹤的足夠長或者到達了渦流的中心
while(curLen < krnlen && advcts < ADVCTS)
{
//**獲得采樣點的矢量資料
vec_id = ( int(clp0_y) * n_xres + int(clp0_x) ) << 1;
vctr_x = pVectr[vec_id ];
vctr_y = pVectr[vec_id + 1];
//**若為關鍵點即一般情況下為渦流的中心時，跳出本次追蹤
if( vctr_x == 0 && vctr_y == 0)
{
t_acum[advDir] = (advcts == 0) ? 0.0f : t_acum[advDir];
w_acum[advDir] = (advcts == 0) ? 1.0f : w_acum[advDir];
break;
}
//**正向追蹤或反向追蹤
vctr_x = (advDir == 0) ? vctr_x : -vctr_x;
vctr_y = (advDir == 0) ? vctr_y : -vctr_y;
segLen = LINE_SQUARE_CLIP_MAX;
segLen = (vctr_x < -VECTOR_COMPONENT_MIN) ? ( int( clp0_x ) - clp0_x ) / vctr_x : segLen;
segLen = (vctr_x > VECTOR_COMPONENT_MIN) ? ( int( int(clp0_x) + 1.5f ) - clp0_x ) / vctr_x : segLen;
segLen = (vctr_y < -VECTOR_COMPONENT_MIN) ?
( ( (tmpLen = ( int( clp0_y ) - clp0_y ) / vctr_y) < segLen) ? tmpLen : segLen) : segLen;
segLen = (vctr_y > VECTOR_COMPONENT_MIN) ?
( ( (tmpLen = ( int( int(clp0_y) + 1.5f ) - clp0_y ) / vctr_y) < segLen) ? tmpLen : segLen) : segLen;
prvLen = curLen;
curLen += segLen;
segLen += 0.0004f;//如何不增加的話，還是會出現問題
//**判斷長度
segLen = (curLen > krnlen) ? ( (curLen = krnlen) - prvLen ) : segLen;
//**獲取下一個追蹤點位置
clp1_x = clp0_x + vctr_x * segLen;
clp1_y = clp0_y + vctr_y * segLen;
//**計算采樣點位置
samp_x = (clp0_x + clp1_x) * 0.5f;
samp_y = (clp0_y + clp1_y) * 0.5f;
///獲取紋理采樣點的灰度值->這里如果采用插值的話或者效果應該又是另一個樣子
texVal = pNoise[ int(samp_y) * n_xres + int(samp_x) ];
W_ACUM = wgtLUT[ int(curLen * len2ID) ];
smpWgt = W_ACUM - w_acum[advDir];
w_acum[advDir] = W_ACUM;
t_acum[advDir] += texVal * smpWgt;
advcts ++;
clp0_x = clp1_x;
clp0_y = clp1_y;
if( clp0_x < 0.0f || clp0_x >= n_xres || clp0_y < 0.0f || clp0_y >= n_yres)
{
break;
}
}//while
}//for
texVal = (t_acum[0] + t_acum[1]) / (w_acum[0] + w_acum[1]);
texVal = (texVal < 0.0f) ? 0.0f : texVal;
texVal = (texVal > 255.0f) ? 255.0f : texVal;
pImage[j * n_xres + i] = (unsigned char)texVal;
}//for
}//for
}
uj5u.com熱心網友回復：
在百度搜這方面的圖片。uj5u.com熱心網友回復：
老師，能不能推薦一些權威的關于流線可視化的論文和書籍，謝謝轉載請註明出處，本文鏈接：https://www.uj5u.com/gongcheng/133319.html
標籤：圖形處理/算法
上一篇：子執行緒通知主執行緒處理結果