// #include <math.h> /***马鞍型矢量场卷积白噪声**稀疏白噪声的产生*/
// #include <stdio.h>
// #include <stdlib.h>
// #include <malloc.h>
// #include <ml.h>
// #include <cxcore.h>
// #include <highgui.h>
// #include <iostream>
// using namespace std;
//
// #define SQUARE_FLOW_FIELD_SZ 400
// #define DISCRETE_FILTER_SIZE 2048
// #define LOWPASS_FILTR_LENGTH 10.00000f
// #define LINE_SQUARE_CLIP_MAX 100000.0f
// #define VECTOR_COMPONENT_MIN 0.050000f
//
//
// void SyntheszSaddle(int n_xres, int n_yres, float* pVectr);
// void NormalizVectrs(int n_xres, int n_yres, float* pVectr);
// void GenBoxFiltrLUT(int LUTsiz, float* p_LUT0, float* p_LUT1);
// void MakeWhiteNoise(int n_xres, int n_yres, float* pNoise);
// void FlowImagingLIC(int n_xres, int n_yres, float* pVectr, float* pNoise,
// float* pImage, float* p_LUT0, float* p_LUT1, float krnlen);
// void WriteImage2PPM(int n_xres, int n_yres, float* pImage, char* f_name);
//
//
// void main()
// {
// //int n_xres = SQUARE_FLOW_FIELD_SZ;
// //int n_yres = SQUARE_FLOW_FIELD_SZ;
//
// //SyntheszSaddle(n_xres, n_yres, pVectr);
// int info[2];
// FILE* fp = fopen("vector256^256.dat", "rb");
// fread(info,sizeof(int),2,fp);
// int col=info[0];
// int row = info[1];
// float* pVectr = (float* ) malloc( sizeof(float ) * col * row * 2 );
// float* p_LUT0 = (float* ) malloc( sizeof(float ) * DISCRETE_FILTER_SIZE);
// float* p_LUT1 = (float* ) malloc( sizeof(float ) * DISCRETE_FILTER_SIZE);
// float* pNoise = (float* ) malloc( sizeof(float) * col * row );
// float* pImage = (float* ) malloc( sizeof(float) * col * row );
//
// cout<<"row="<<row<<";"<<"col="<<col<<endl;
// //float * pVectr = new float[row*col*2];
// fread(pVectr,sizeof(float),row*col*2,fp);
//
//
// //NormalizVectrs(n_xres, n_yres, pVectr);
// MakeWhiteNoise( row,col, pNoise);
// GenBoxFiltrLUT(DISCRETE_FILTER_SIZE, p_LUT0, p_LUT1);
// FlowImagingLIC(row,col, pVectr, pNoise, pImage, p_LUT0, p_LUT1, LOWPASS_FILTR_LENGTH);
// WriteImage2PPM(row, col, 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;
// fclose(fp);
// fp=NULL;
// //system("pause");
// }
//
//
// /// synthesize a saddle-shaped vector field ///
// void SyntheszSaddle(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 = ( (n_yres - 1 - j) * n_xres + i ) << 1;
// pVectr[index ] = - ( j / (n_yres - 1.0f) - 0.5f );
// pVectr[index + 1] = i / (n_xres - 1.0f) - 0.5f;
// }
//
// }
//
//
// /// normalize the vector field ///
// 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;
// }
// }
//
//
// /// make white noise as the LIC input texture ///
// void MakeWhiteNoise(int n_xres, int n_yres, float* pNoise)
// {
// IplImage * NoiseImg=cvCreateImage(cvSize(n_xres,n_yres),IPL_DEPTH_8U,1);
// CvScalar s;
//
// for(int j = 0; j < n_xres; j =j++)
// for(int i = 0; i < n_yres; i=i++)//产生白噪声，改变j = j+10;i = i+10
// {
// int r = rand();
// r = ( (r & 0xff) + ( (r & 0xff00) >> 8 ) ) & 0xff;//参数8的改变
// pNoise[j * n_xres + i] = (float) r;
// s = cvGet2D(NoiseImg,i,j);
// s.val[0]=r;
// s.val[1]=r;
// s.val[2]=r;
// cvSet2D(NoiseImg,i,j,s);
// }
// // cvShowImage("NoiseImg",NoiseImg);
// // cvWaitKey(0);
// // cvDestroyWindow("NoiseImg");
// // cvReleaseImage(&NoiseImg);
// }
//
//
// /// generate box filter LUTs ///
// void GenBoxFiltrLUT(int LUTsiz, float* p_LUT0, float* p_LUT1)
// {
// for(int i = 0; i < LUTsiz; i ++) p_LUT0[i] = p_LUT1[i] = i;
// }
//
//
// /// write the LIC image to a PPM file ///
// void WriteImage2PPM(int n_xres, int n_yres, float* 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_xres; j ++)
// for(int i = 0; i < n_yres; i ++)
// {
// unsigned char unchar = pImage[j * n_xres + i];
// fprintf(o_file, "%c%c%c", unchar, unchar, unchar);
// }
//
// fclose (o_file); o_file = NULL;
// }
//
//
// /// flow imaging (visualization) through Line Integral Convolution ///
// void FlowImagingLIC(int n_xres, int n_yres, float* pVectr, float* pNoise, float* pImage,
// float* p_LUT0, float* p_LUT1, float krnlen)
// {
// int vec_id; ///ID in the VECtor buffer (for the input flow field)
// int advDir; ///ADVection DIRection (0: positive; 1: negative)
// int advcts; ///number of ADVeCTion stepS per direction (a step counter)
// int ADVCTS = int(krnlen * 3); ///MAXIMUM number of advection steps per direction to break dead loops
//
// float vctr_x; ///x-component of the VeCToR at the forefront point
// float vctr_y; ///y-component of the VeCToR at the forefront point
// float clp0_x; ///x-coordinate of CLiP point 0 (current)
// float clp0_y; ///y-coordinate of CLiP point 0 (current)
// float clp1_x; ///x-coordinate of CLiP point 1 (next )
// float clp1_y; ///y-coordinate of CLiP point 1 (next )
// float samp_x; ///x-coordinate of the SAMPle in the current pixel
// float samp_y; ///y-coordinate of the SAMPle in the current pixel
// float tmpLen; ///TeMPorary LENgth of a trial clipped-segment
// float segLen; ///SEGment LENgth
// float curLen; ///CURrent LENgth of the streamline
// float prvLen; ///PReVious LENgth of the streamline
// float W_ACUM; ///ACcuMulated Weight from the seed to the current streamline forefront
// float texVal; ///TEXture VALue
// float smpWgt; ///WeiGhT of the current SaMPle
// float t_acum[2]; ///two ACcUMulated composite Textures for the two directions, perspectively
// float w_acum[2]; ///two ACcUMulated Weighting values for the two directions, perspectively
// float* wgtLUT = NULL; ///WeiGhT Look Up Table pointing to the target filter LUT
// float len2ID = (DISCRETE_FILTER_SIZE - 1) / krnlen; ///map a curve LENgth TO an ID in the LUT
//
// ///for each pixel in the 2D output LIC image///
// for(int j = 0; j < n_xres; j ++)
// for(int i = 0; i < n_yres; i ++)
// {
// ///init the composite texture accumulators and the weight accumulators///
// t_acum[0] = t_acum[1] = w_acum[0] = w_acum[1] = 0.0f;
//
// ///for either advection direction///
// for(advDir = 0; advDir < 2; advDir ++)
// {
// ///init the step counter, curve-length measurer, and streamline seed///
// advcts = 0;
// curLen = 0.0f;
// clp0_x = j + 0.5f;
// clp0_y = i + 0.5f;
//
// ///access the target filter LUT///
// wgtLUT = (advDir == 0) ? p_LUT0 : p_LUT1;
//
// ///until the streamline is advected long enough or a tightly spiralling center / focus is encountered///
// while( curLen < krnlen && advcts < ADVCTS )
// {
// ///access the vector at the sample///
// vec_id = ( int(clp0_y) * n_xres + int(clp0_x) ) << 1;
// vctr_x = pVectr[vec_id ];
// vctr_y = pVectr[vec_id + 1];
//
// ///in case of a critical point///
// if( vctr_x == 0.0f && vctr_y == 0.0f )
// {
// t_acum[advDir] = (advcts == 0) ? 0.0f : t_acum[advDir]; ///this line is indeed unnecessary
// w_acum[advDir] = (advcts == 0) ? 1.0f : w_acum[advDir];
// break;
// }
//
// ///negate the vector for the backward-advection case///
// vctr_x = (advDir == 0) ? vctr_x : -vctr_x;
// vctr_y = (advDir == 0) ? vctr_y : -vctr_y;
//
// ///clip the segment against the pixel boundaries --- find the shorter from the two clipped segments///
// ///replace all if-statements whenever possible as they might affect the computational speed///
// 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;
//
// ///update the curve-length measurers///
// prvLen = curLen;
// curLen+= segLen;
// segLen+= 0.0004f;
//
// ///check if the filter has reached either end///
// segLen = (curLen > krnlen) ? ( (curLen = krnlen) - prvLen ) : segLen;
//
// ///obtain the next clip point///
// clp1_x = clp0_x + vctr_x * segLen;
// clp1_y = clp0_y + vctr_y * segLen;
//
// ///obtain the middle point of the segment as the texture-contributing sample///
// samp_x = (clp0_x + clp1_x) * 0.5f;
// samp_y = (clp0_y + clp1_y) * 0.5f;
//
// ///obtain the texture value of the sample///
// texVal = pNoise[ int(samp_y) * n_xres + int(samp_x) ];
//
// ///update the accumulated weight and the accumulated composite texture (texture x weight)///
// W_ACUM = wgtLUT[ int(curLen * len2ID) ];
// smpWgt = W_ACUM - w_acum[advDir];
// w_acum[advDir] = W_ACUM;
// t_acum[advDir] += texVal * smpWgt;
//
// ///update the step counter and the "current" clip point///
// advcts ++;
// clp0_x = clp1_x;
// clp0_y = clp1_y;
//
// ///check if the streamline has gone beyond the flow field///
// if( clp0_x < 0.0f || clp0_x >= n_xres || clp0_y < 0.0f || clp0_y >= n_yres) break;
// }
// }
//
// ///normalize the accumulated composite texture///
// texVal = (t_acum[0] + t_acum[1]) / (w_acum[0] + w_acum[1]);
//
// ///clamp the texture value against the displayable intensity range [0, 255]
// texVal = (texVal < 0.0f) ? 0.0f : texVal;
// texVal = (texVal > 255.0f) ? 255.0f : texVal;
// pImage[i * n_xres + j] = (float) texVal;
// }
// }
/************中心环形矢量场*马鞍矢量场*****卷积白噪声纹理***********/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <cv.h>
#include <ml.h>
#include <highgui.h>
#include <iostream>
using namespace std;
//#define SQUARE_FLOW_FIELD_SZ 400
#define DISCRETE_FILTER_SIZE 2048 //离散的滤波尺寸,
#define LOWPASS_FILTR_LENGTH 10.00000f //低通滤波长度，滤波核kernel
#define LINE_SQUARE_CLIP_MAX 100000.0f //线性平方夹
#define VECTOR_COMPONENT_MIN 0.050000f //矢量分量最小值
void CenterFiled(int n_xres, int n_yres, float* pVectr);
void NormalizVectrs(int n_xres, int n_yres, float* pVectr,float* vecSize,float*normMag);
void GenBoxFiltrLUT(int LUTsiz, float* p_LUT0, float* p_LUT1);
void MakeWhiteNoise(int n_xres, int n_yres, float* pNoise);
void nNoise(int n_xres, int n_yres, float* pNoise,float*pVectr,float*newNoise,float* normMag);
void FlowImagingLIC(int n_xres, int n_yres, float* pVectr, float* newNoise,
float* pImage, float* p_LUT0, float* p_LUT1, float krnlen);
void WriteImage2PPM(int n_xres, int n_yres, float* pImage, char* f_name);
void gray(int n_xres,int n_yres,float * pImage);
void color(int n_xres, int n_yres,float *pImage,float* vecSize);
double maxvecmag;
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);
// float* pNoise = (float* ) malloc( sizeof(float) * n_xres * n_yres );
// float* pImage = (float* ) malloc( sizeof(float) * n_xres * n_yres );
//CenterFiled(n_xres, n_yres, pVectr);
//NormalizVectrs(n_xres, n_yres, pVectr);
int info[2];
FILE* fp = fopen("fieldfile/Flow_16.vec", "rb"); //打开矢量场数据文件
fread(info,sizeof(int),2,fp);
int n_yres=info[0];
int n_xres = info[1];
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);
float* pNoise = (float* ) malloc( sizeof(float) *n_xres*n_yres );
float* pImage = (float* ) malloc( sizeof(float) * n_xres*n_yres );
float* vecSize = (float* ) malloc( sizeof(float) * n_xres*n_yres );
float* newNoise = (float* ) malloc( sizeof(float) *n_xres*n_yres );
float* normMag = (float* ) malloc( sizeof(float) *n_xres*n_yres );
// cout<<"row="<<n_xres<<";"<<"col="<<n_yres<<endl;
// //float * pVectr = new float[row*col*2];
fread(pVectr,sizeof(float),n_xres*n_yres*2,fp);//读入矢量场数据，是未归一化的矢量场
NormalizVectrs(n_xres, n_yres, pVectr,vecSize,normMag);//利用刘占平的数据文件必须有归一化
MakeWhiteNoise(n_xres, n_yres, pNoise);
nNoise(n_xres,n_yres,pNoise,pVectr,newNoise,normMag);
GenBoxFiltrLUT(DISCRETE_FILTER_SIZE, p_LUT0, p_LUT1);
FlowImagingLIC(n_xres, n_yres, pVectr, newNoise, pImage, p_LUT0, p_LUT1, LOWPASS_FILTR_LENGTH);
gray(n_xres,n_yres,pImage);
color(n_xres, n_yres,pImage,vecSize);
//WriteImage2PPM(n_xres, n_yres, pImage, "LIC.ppm");
//system("pause");
fclose(fp);
fp=NULL;
free(pVectr); pVectr = NULL;
free(p_LUT0); p_LUT0 = NULL;
free(p_LUT1); p_LUT1 = NULL;
free(pNoise); pNoise = NULL;
free(pImage); pImage = NULL;
}
/// 中心环形矢量场图形 synthesize a saddle-shaped vector field ///
void CenterFiled(int n_xres, int n_yres, float* pVectr)
{
float vec_x = 0.0f;
float vec_y = 0.0f;
float vcMag = 0.0f;
float scale = 0.0f;
for(int i = 0; i < n_xres; i ++)
{
for(int j = 0; j < n_yres; j ++)
{
int index = i*n_yres+j;
vec_x = -(float)i/n_xres+0.5f;
vec_y = (float)j/n_yres-0.5f;
// vec_x = -i/n_xres+0.5f;//如果不进行float形式转化，生成的矢量场不正确
// vec_y = j/n_yres-0.5f;
vcMag = 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;
// int index = ( (n_yres - 1 - j) * n_xres + i ) << 1;//向左移动一位
// pVectr[index ] = - ( j / (n_yres - 1.0f) - 0.5f );
// cout<<"pVectr[index ]="<<pVectr[index ];
// pVectr[index + 1] = i / (n_xres - 1.0f) - 0.5f;
// cout<<"pVectr[index ]="<<pVectr[index +1 ];
}
}
}
/// normalize the vector field ///
void NormalizVectrs(int n_xres, int n_yres, float* pVectr,float* vecSize,float*normMag)
{
// float* vecSize = (float* ) malloc( sizeof(float) * n_xres*n_yres );
double magind;
double scale;
double mag;
double x=10;;
矢量归一化**与**矢量大小归一化不一样
//矢量归一化如下：
for(int j = 0; j < n_yres; j ++)
for(int i = 0; i < n_xres; i ++)//图像时竖着绘制的，一列一列的画的
{
int index = (j * n_xres + i) << 1;//0,2,4,6,8...向左移一位，index值加2，因为有vec_x,vec_y代表一个像素点的矢量值
float vcMag = float( sqrt( double(pVectr[index] * pVectr[index] + pVectr[index + 1] * pVectr[index + 1]) ) );
vecSize[j * n_xres + i]=vcMag;//存储矢量原始大小
//cout<<"vcmag="<<vcMag<<endl;
if (vcMag>maxvecmag)
{
maxvecmag=vcMag;
}
//cout<<vcMag<<endl;
float scale = (vcMag == 0.0f) ? 0.0f : 1.0f / vcMag;//矢量大小归一化后的矢量值
//cout<<"scale"<<scale<<endl;
//normMag[j*n_yres+i] = scale;
pVectr[index ]=pVectr[index ]*scale;
pVectr[index + 1] *= scale;
}
for(int j = 0; j < n_yres; j ++)
for(int i = 0; i < n_xres; i ++)//图像时竖着绘制的，一列一列的画的
{
//归一化矢量大小
normMag[j * n_xres + i] = (float)vecSize[j * n_xres + i]/maxvecmag;
//cout<< normMag[j * n_xres + i]<<endl;
}
}
/// make white noise as the LIC input texture ///
void MakeWhiteNoise(int n_xres, int n_yres, float* pNoise)
{
IplImage * NoiseImg=cvCreateImage(cvSize(n_yres,n_xres),IPL_DEPTH_8U,1);
CvScalar s;
for(int j = 0; j < n_yres; j ++)
{
for(int i = 0; i < n_xres; i ++)
// for(int j = 0; j < n_yres; j=j +10)//产生稀疏白噪声
// {
// for(int i = 0; i < n_xres; i=i+10)
{
int r = rand();
r = ( (r & 0xff) + ( (r & 0xff00) >> 8 ) ) & 0xff;
pNoise[j * n_xres + i] = (float) r;
s = cvGet2D(NoiseImg,i,j);
s.val[0]=r;
s.val[1]=r;
s.val[2]=r;
cvSet2D(NoiseImg,i,j,s);
}
}
}
void nNoise(int n_xres, int n_yres, float*pNoise,float* pVectr,float*newNoise,float* normMag)
{
for(int j = 0; j < n_yres; j ++)
{
for(int i = 0; i < n_xres; i ++)
{
newNoise[j*n_yres+i]=255*normMag[j*n_yres+i]+pNoise[j*n_yres+i]-128;
}
}
}
/// generate box filter LUTs ///
void GenBoxFiltrLUT(int LUTsiz, float* p_LUT0, float* p_LUT1)
{
for(int i = 0; i < LUTsiz; i ++) p_LUT0[i] = p_LUT1[i] = i;
}
/// flow imaging (visualization) through Line Integral Convolution ///
void FlowImagingLIC(int n_xres, int n_yres, float* pVectr, float* newNoise, float* pImage,
float* p_LUT0, float* p_LUT1, float krnlen)
{
int vec_id; ///ID in the VECtor buffer (for the input flow field)
int advDir; ///ADVection DIRection (0: positive; 1: negative)
int advcts; ///number of ADVeCTion stepS per direction (a step counter)
int ADVCTS = int(krnlen * 3); ///MAXIMUM number of advection steps per direction to break dead loops //跳出死循环的条件
float vctr_x; ///x-component of the Vector at the forefront point
float vctr_y; ///y-component of the Vector at the forefront point
float clp0_x; ///x-coordinate of CLiP point 0 (current)
float clp0_y; ///y-coordinate of CLiP point 0 (current)
float clp1_x; ///x-coordinate of CLiP point 1 (next )
float clp1_y; ///y-coordinate of CLiP point 1 (next )
float samp_x; ///x-coordinate of the Sample in the current pixel
float samp_y; ///y-coordinate of the Sample in the current pixel
float tmpLen; ///Temporary Length of a trial clipped-segment实验剪辑片段的临时长度
float segLen; ///Segment Length
float curLen; ///Current Length of the streamline
float prvLen; ///Previous Length of the streamline
float W_ACUM; ///Accumulated Weight from the seed to the current streamline forefront
float texVal; ///Texture Value
float smpWgt; ///Weight of the current Sample
float t_acum[2]; ///two Accumulated composite Textures for the two directions, perspectively 两个方向的纹理卷积累加和
float w_acum[2]; ///two Accumulated Weighting values for the two directions, perspectively 两个方向的权重和
float* wgtLUT = NULL; ///WeiGhT Look Up Table pointing to the target filter LUT权重查找表
float len2ID = (DISCRETE_FILTER_SIZE - 1) / krnlen; ///map a curve Length To an ID in the LUT
double scale;//颜色映射表比例
double maxmag;
double magind;
double mag;
double x = 0.1;//x为非线性映射因子，且x！=1
IplImage * licImage = cvCreateImage(cvSize(n_yres,n_xres),IPL_DEPTH_8U,3);
CvScalar s;
///for each pixel in the 2D output LIC image///对输出图像的每一个像素
for(int j = 0; j < n_yres; j ++)
{
for(int i = 0; i < n_xres; i ++)
{
///init the composite texture accumulators and the weight accumulators///每一个像素点为起始点，初始化一次权重和卷积和
t_acum[0] = t_acum[1] = w_acum[0] = w_acum[1] = 0.0f;//初始化正反方向卷积和及权重和
///for either advection direction///分别计算正反方向的卷积和及权重和
for(advDir = 0; advDir < 2; advDir ++)
{
///init the step counter, curve-length measurer, and streamline seed///
//初始化当前方向上前进的步数和当前流线的总长
advcts = 0;//前进的步数
curLen = 0.0f;
clp0_x = i + 0.5f;//当前点的坐标
clp0_y = j + 0.5f;
///access the target filter LUT///LUT显示查找表
wgtLUT = (advDir == 0) ? p_LUT0 : p_LUT1;//当前噪声点所对应的权重系数
///until the streamline is advected long enough or a tightly spiralling center / focus is encountered///
while( curLen < krnlen && advcts < ADVCTS ) //??????
{
///access the vector at the sample///
vec_id = ( int(clp0_y) * n_xres + int(clp0_x) ) << 1;//vec_id相当于index
vctr_x = pVectr[vec_id ];//clp0_y相当于当前像素列坐标，clp0_x相当于当前像素的横坐标
vctr_y = pVectr[vec_id + 1];
///in case of a critical point///遇到零向量，结束循环
if( vctr_x == 0.0f && vctr_y == 0.0f )
{
t_acum[advDir] = (advcts == 0) ? 0.0f : t_acum[advDir]; ///this line is indeed unnecessary
w_acum[advDir] = (advcts == 0) ? 1.0f : w_acum[advDir];
break;
}
///negate the vector for the backward-advection case///相反的方向取相反的方向
vctr_x = (advDir == 0) ? vctr_x : -vctr_x;//因为矢量是用x,y方向的值合成的，所以反向的值就是负的
vctr_y = (advDir == 0) ? vctr_y : -vctr_y;
//这儿可以就计算矢量大小、归一化运算？？？？？？？？？？？？前面已经归一化了
//。。。。。。//
// mag= sqrt(vctr_x*vctr_x+vctr_y*vctr_y);
// if (mag>maxmag)
// {
// maxmag=mag;
// }
///clip the segment against the pixel boundaries --- find the shorter from the two clipped segments///
///replace all if-statements whenever possible as they might affect the computational speed///影响算法计算速度
segLen = LINE_SQUARE_CLIP_MAX;
segLen = (vctr_x < -VECTOR_COMPONENT_MIN) ? ( int( clp0_x ) - clp0_x ) / vctr_x : segLen;//int（0.5）=0
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;
///update the curve-length measurers///
prvLen = curLen;
curLen+= segLen;//segLen是单步步长
segLen+= 0.0004f;//步长？？？
// segLen+= 0.0001f;//步长？？？
///check if the filter has reached either end///
segLen = (curLen > krnlen) ? ( (curLen = krnlen) - prvLen ) : segLen;
///obtain the next clip point///
clp1_x = clp0_x + vctr_x * segLen;
clp1_y = clp0_y + vctr_y * segLen;
///obtain the middle point of the segment as the texture-contributing sample///
samp_x = (clp0_x + clp1_x) * 0.5f;
samp_y = (clp0_y + clp1_y) * 0.5f;
///obtain the texture value of the sample///
texVal = newNoise[ int(samp_y) * n_xres + int(samp_x) ];
///update the accumulated weight and the accumulated composite texture (texture x weight)///
W_ACUM = wgtLUT[ int(curLen * len2ID) ];
smpWgt = W_ACUM - w_acum[advDir];
w_acum[advDir] = W_ACUM;
t_acum[advDir] += texVal * smpWgt;//当前噪声点的权重系数
///update the step counter and the "current" clip point///
advcts ++;
clp0_x = clp1_x;
clp0_y = clp1_y;
///check if the streamline has gone beyond the flow field///
if( clp0_x < 0.0f || clp0_x >= n_xres || clp0_y < 0.0f || clp0_y >= n_yres) break;
}
}
///normalize the accumulated composite texture///
texVal = (t_acum[0] + t_acum[1]) / (w_acum[0] + w_acum[1]);
///clamp the texture value against the displayable intensity range [0, 255]
texVal = (texVal < 0.0f) ? 0.0f : texVal/255;
//cout<<"texval="<<texVal<<endl;
texVal = (texVal > 255.0f) ? 255.0f : texVal;
//cout<<texVal<<endl;
pImage[j * n_xres + i] = (float) texVal;
}
}
}
void gray(int n_xres,int n_yres,float *pImage)
{
CvScalar s;
IplImage * GrayImage = cvCreateImage(cvSize(n_yres,n_xres),IPL_DEPTH_8U,1);
for (int i = 0;i<GrayImage->width;i++)
{
for (int j = 0;j<GrayImage->height;j++)
{
s.val[0] = pImage[i*n_xres+j]*255;
s.val[1] = pImage[i*n_xres+j]*255;
s.val[2] = pImage[i*n_xres+j]*255;
//cout<<s.val[0]<<endl;
cvSet2D(GrayImage,i,j,s);
}
}
cvSaveImage("gray.jpg",GrayImage,0);
cvWaitKey(0);
IplImage *gray = cvLoadImage("gray.jpg",1);
//直方图尺寸
int r_bins =256, b_bins = 256;
CvHistogram* hist;
{
int hist_size[] = { r_bins, b_bins };
float r_ranges[] = { 0, 255 }; // hue is [0,180]
float b_ranges[] = { 0, 255 };
float* ranges[] = { r_ranges,b_ranges };
hist = cvCreateHist( 2, hist_size, CV_HIST_ARRAY, ranges,1);
}
cvReleaseImage(&GrayImage);
}
// void color(int n_xres,int n_yres,float *pImage,float* vecSize)
// {
// IplImage * licImage = cvCreateImage(cvSize(n_yres,n_xres),IPL_DEPTH_8U,3);
// IplImage* img = cvLoadImage("11.jpg",1);
// //IplImage* destImg = cvCreateImage(cvSize(n_yres,n_xres),IPL_DEPTH_8U,3);
// //int length=img->width;
// int k = 0;
//
// double magind;
// double mag;
// double newMag;
// double x = 0.01;//x为非线性映射因子，且x！=1
//
// CvScalar colorTable[500];
// CvScalar s,s1;
// //cout<<"h="<<img->height<<";"<<"w"<<img->width<<endl;
// //system("pause");
// for (int i = 0;i < img->width;i++)
// {
// s = cvGet2D(img,1,i);
// colorTable[i] =s;
// //cout<<"colorTable[i]="<<colorTable[i].val[0]<<endl;
// }
// for (int i = 0;i<n_xres;i++)
// {
// for (int j= 0; j <n_yres;j++)
// {
//
// if (k>=img->width)
// {
// k=0;
// }
// double scale= pImage[j * n_xres + i];
//
//
// mag = vecSize[j * n_xres + i];
// //cout<<"mag="<<mag<<endl;
// //********矢量大小归一化******
// magind = mag/maxvecmag;
// //cout<<"maxvecmag="<<maxvecmag<<endl;
// //cout<<"magind="<<magind<<endl;
// //非线性颜色增强LIC
// newMag = (pow(x,magind)-1)/(x-1);
// //cout<<"newMag="<<newMag<<endl;
// // cout<<"scale="<<scale<<endl;
//
// //最终的输出颜色值计算如下
// /*color = table[newMag];*/
//
// s = cvGet2D(licImage,i,j);
//
// //渐变颜色映射表
// int k = int(newMag*446);
// s1.val[0]=colorTable[k].val[0]*(k+1-newMag*446)+colorTable[k+1].val[0]*(newMag*446-k);
// s1.val[1]=colorTable[k].val[1]*(k+1-newMag*446)+colorTable[k+1].val[1]*(newMag*446-k);
// s1.val[2]=colorTable[k].val[2]*(k+1-newMag*446)+colorTable[k+1].val[2]*(newMag*446-k);
// s1.val[0]*=scale;
// s1.val[1]*=scale;
// s1.val[2]*=scale;
//
//
// //离散型颜色映射表
// // int destcolorIndext=int(446*mag*100);
// // s1.val[0]=colorTable[destcolorIndext].val[0]*scale;
// // s1.val[1]=colorTable[destcolorIndext].val[1]*scale;
// // s1.val[2]=colorTable[destcolorIndext].val[2]*scale;
//
// cvSet2D(licImage,i,j,s1);
// //cout<<"s1.val[0]="<<s1.val[0]<<endl;
// }
// }
// cvSaveImage("color.jpg",licImage,0);
//
// cvNamedWindow("lic_three channles",0 );
// cvShowImage("lic_three channles",licImage);
// cvWaitKey(0);
// system("pause");
// cvDestroyWindow("lic_three channles");
//
// cvReleaseImage(&licImage);
// }
/// write the LIC image to a PPM file ///
void color(int n_xres,int n_yres,float *pImage,float* vecSize)
{
IplImage * licImage = cvCreateImage(cvSize(n_yres,n_xres),IPL_DEPTH_8U,3);
IplImage* img = cvLoadImage("11.jpg",1);
//IplImage* destImg = cvCreateImage(cvSize(n_yres,n_xres),IPL_DEPTH_8U,3);
//int length=img->width;
int k = 0;
double magind;
double mag;
double newMag;
double x = 0.01;//x为非线性映射因子，且x！=1
CvScalar colorTable[500];
CvScalar s,s1;
//cout<<"h="<<img->height<<";"<<"w"<<img->width<<endl;
//system("pause");
for (int i = 0;i < img->width;i++)
{
s = cvGet2D(img,1,i);
colorTable[i] =s;
//cout<<"colorTable[i]="<<colorTable[i].val[0]<<endl;
}
for (int i = 0;i<n_xres;i++)
{
for (int j= 0; j <n_yres;j++)
{
if (k>=img->width)
{
k=0;
}
double scale= pImage[j * n_xres + i];
mag = vecSize[j * n_xres + i];
//cout<<"mag="<<mag<<endl;
//********矢量大小归一化******
magind = mag/maxvecmag;
//cout<<"maxvecmag="<<maxvecmag<<endl;
//cout<<"magind="<<magind<<endl;
//非线性颜色增强LIC
newMag = (pow(x,magind)-1)/(x-1);
//cout<<"newMag="<<newMag<<endl;
// cout<<"scale="<<scale<<endl;
//最终的输出颜色值计算如下
/*color = table[newMag];*/
s = cvGet2D(licImage,i,j);
//渐变颜色映射表
int k = int(newMag*446);
s1.val[0]=colorTable[k].val[0]*(k+1-newMag*446)+colorTable[k+1].val[0]*(newMag*446-k);
s1.val[1]=colorTable[k].val[1]*(k+1-newMag*446)+colorTable[k+1].val[1]*(newMag*446-k);
s1.val[2]=colorTable[k].val[2]*(k+1-newMag*446)+colorTable[k+1].val[2]*(newMag*446-k);
s1.val[0]*=scale;
s1.val[1]*=scale;
s1.val[2]*=scale;
//离散型颜色映射表
// int destcolorIndext=int(446*mag*100);
// s1.val[0]=colorTable[destcolorIndext].val[0]*scale;
// s1.val[1]=colorTable[destcolorIndext].val[1]*scale;
// s1.val[2]=colorTable[destcolorIndext].val[2]*scale;
cvSet2D(licImage,i,j,s1);
//cout<<"s1.val[0]="<<s1.val[0]<<endl;
}
}
cvSaveImage("color.jpg",licImage,0);
cvNamedWindow("lic_three channles",0 );
cvShowImage("lic_three channles",licImage);
cvWaitKey(0);
system("pause");
cvDestroyWindow("lic_three channles");
cvReleaseImage(&licImage);
}
void WriteImage2PPM(int n_xres, int n_yres, float* 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];//某点像素的灰度纹理值
CvScalar s;
/*
s.val[0]
s.val[1]
s.val[2]
fprintf(o_file, "%c%c%c", s.val[0], s.val[1], s.val[2]);//
*/
fprintf(o_file, "%c%c%c", unchar, unchar, unchar);//
}
fclose (o_file); o_file = NULL;
}
11-06
790
790
“相关推荐”对你有帮助么？
-
非常没帮助 -
没帮助 -
一般 -
有帮助 -
非常有帮助
提交
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
