noss_haru 2015-10-18 12:55 采纳率: 0%
浏览 3688
已结题

LIC(线积分卷积)源代码，怎么结合具体图像怎么使用？

运行没有错误，刚接触图形学，不知道如何使用这段代码。
/*----------------------------------------------------------*/
/* 代码参考如下链接: /
/ http://www.zhanpingliu.org/Research/FlowVis/LIC/LIC.htm /
/ /
/----------------------------------------------------------*/

#include
#include
#include

#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

}
/*----------------------------------------------------------------------------------------------------------*/

  • 写回答

4条回答

  • threenewbee 2015-10-18 15:42
    关注

    有注释，有函数，直接调用呗。

    评论

报告相同问题？

悬赏问题

  • ¥15 深度学习根据CNN网络模型，搭建BP模型并训练MNIST数据集
  • ¥15 lammps拉伸应力应变曲线分析
  • ¥15 C++ 头文件/宏冲突问题解决
  • ¥15 用comsol模拟大气湍流通过底部加热（温度不同）的腔体
  • ¥50 安卓adb backup备份子用户应用数据失败
  • ¥20 有人能用聚类分析帮我分析一下文本内容嘛
  • ¥15 请问Lammps做复合材料拉伸模拟，应力应变曲线问题
  • ¥30 python代码，帮调试，帮帮忙吧
  • ¥15 #MATLAB仿真#车辆换道路径规划
  • ¥15 java 操作 elasticsearch 8.1 实现 索引的重建