벡터필드 : 풍속 데이터로 바람 그리기

풍속 데이터로 바람을 표현해보자. 이 글은 두 개의 레퍼런스와 본 글로 구성된다.

이 글에서는 OpenGL을 이용해서 그리는 방식을 설명한다. 필요한 데이터와 기본적인 원리들을 설명하고 있으므로, 다른 방식으로 그리려면 글 내용을 이해한 후 스스로 시도해봐야 한다. 그리기 위한 전체 코드를 설명하지는 않고 핵심적인 부분만 다루므로 OpenGL의 기본적인 부분들은 알고 있어야 한다.

 

예시 1 : Cameron Beccario의 <earth> 

 

아래 사이트에 들어가면 둥근 지구 표면 위에서 바람이 흐르는 모습을 볼 수 있다. (바람 이외에도 여러가지 기상 정보를 볼 수 있다)

earth :: a global map of wind, weather, and ocean conditions

 

우리나라 부분을 캡쳐해보면 아래 화면처럼 보인다.

Cameron Beccario 는 대기, 해수, 화학물질 등 여러가지 데이터를 바탕으로 지구 위에 다양한 정보들을 그려내고 있다. 그 중 바람 시각화는 NOAA(National Oceanic and Atmospheric Administration)에서 업데이트하는 3시간 간격의 데이터를 바탕으로 만들었다고 한다. (https://earth.nullschool.net/ko/about.html)

 

이 데이터를 보면, 경위도 1도 간격으로 바람의 속도가 2차원 벡터로 표현되어 있다. 이런 종류의 데이터 시각화는 랜덤하게 한번 그린 점들의 다음 위치를, 지구 전체에 깔린 벡터필드(vector field)를 바탕으로 연산하여 만들어내게 된다. 

 

시각화 구현의 핵심적인 내용들은 github의 <implemention note> 에 적혀 있다.(https://github.com/cambecc/earth) d3.js를 사용했으며 cpu 를 이용한 연산을 하고 있다.

 

 

예시 2 : Vladimir Agafonkin 의 WebGL 방식

 

아래 링크를 보면 Mapbox 엔지니어인 Vladimir Agafonkin이 비슷한 방식의 시각화를,  WebGL로 구현한 내용을 설명하고 있다.

How I built a wind map with WebGL

GPU 연산을 통해 비슷한 결과물을 만들어내고 있으며, 아래 사이트에서 최대 58만여개의 파티클들로 만들어진 시각화를 볼 수 있다.

 

WebGL wind simulation

여기서 설명하는 파티클 시각화 방식은 대부분의 파티클 시뮬레이션에서 사용하는 일반적인 방식이다. 두 개의 버퍼(A, B라고 하자)를 이용하는데, 0번 프레임에서는 A 버퍼를 이용해서 그 다음 위치를 계산한 후 B 버퍼에 쓰고, 1번 프레임에서는 B 버퍼를 읽어서 다음 위치를 계산한 후 다시 A 버퍼에 쓴다.

많은 파티클들을 동시에 계산하기 때문에 파티클 간의 상호작용을 계산할 경우 하나의 버퍼에서 읽고 바로 쓰게 되면, 어떤 입자는 다른 입자가 이미 이동한 위치와의 관계속에서 계산 될 수도 있다. 그러므로 프레임을 철저히 분리해서 한 프레임은 이번 계산 이전 시점의 위치, 다른 프레임은 이번 계산이 반영된 시점의 위치로 활용하게 된다.

 

그런데 이 바람 시각화에서는 입자 상호간의 관계를 연산하는 것이 아닌데, 왜 두 개의 버퍼를 사용했지???? 잘 모르겠다. WebGL이 지닌 한계 때문에 그러한 것인가?? 코드를 살펴보진 않았지만, 어쨌든 OpenGL에서는 버퍼 하나만 사용해서 바로 읽은 후 그 버퍼의 자리에 업데이트 된 위치를 써도 아무런 문제가 없다! 버퍼 두 개를 swap 해가면서 쓰는 것은 입자 상호간의 관계를 계산할 시간 전후 순서가 엉키지 않도록 하기 위해서 필요하기 때문이다.

 

 

 

 

이제 직접 만들어보자

 

여기서는 비슷한 결과물을 OpenGL 방식으로 만들어 본 과정을 서술해보겠다. 시점을 움직일 때 화면의 흐름이 끊기고 다시 점을 생성하는 위의 구현 방식들과는 달리, 시점을 옮기더라도 바람의 흐름이 그대로 이어질 수 있도록 해봤다. 물론 매 프레임마다 루프를 돌면서 선형 shape 을 단계적으로 계산해줘야 하므로 GPU 연산 비용은 많이 들어간다. (이 방식은 개선에 대한 여지가 있는데, 마지막에 설명해보겠다)

 

 

 

 

데이터 받기

 

 

데이터는 아래 링크에서 받을 수 있다. 전 세계의 기상 관측과 예측 데이터를 제공한다.

 

 

NCEP Data Products GFS and GDAS

 

 

그런데 그냥 이용하기에는 다소 번거로우므로 다시 아래의 링크를 들어가보자. 위의 데이터를 이용하기 편리한 방식으로 변환해준다. 

 

ERDDAP - NOAA/NCEP Global Forecast System (GFS) Atmospheric Model - Make A Graph

 

Graph Type은 vectors로 설정.

vector X 는 ugrd10m, vector Y는 vgrd10m 로 바꿔준다. 고도 10m, 즉 지표면의 바람이며 동쪽, 북쪽 방향의 벡터로 표현했다. 흔히 보는 xy 평면 방향의 벡터를 동서남북 방향과 일치시켰다고 생각하면 된다.

 

영역은 전 지구를 받아야 하니, 위의 그림에서 설정된대로(기본값이다) 둔다. 경도 360도와 0도는 같으므로 겹치지 않도록 359.5도까지만 제공한다.

마지막에 형식은 tsv0 으로 두었다. 헤더 없는 탭 구분자의 데이터다.

 

3시간 단위의 관측/예측 데이터를 한 벌씩 받을 수 있는데, 좀 더 긴 기간이 필요할 경우 오른쪽 상단의 Data Access Form 을 클릭해보자. 아래와 같은 화면이 보인다.

여기서 시간을 설정하면 된다. 한 시점의 데이터가 텍스트 기준으로 13MB정도 되며 수백MB 데이터까지도 받을 수 있다. 더 큰 용량은 시도해보지 않았다.

 

 

 

 

데이터

 

 

데이터는 위의 그림처럼 생겼다.

UTC +0 시간대 기준이며, 위도, 경도(0~395.5), 그리고 서풍 방향(ugrd10m), 남풍 방향(vgrd10m)의 속력 데이터가 있다. 흔히 사용하는 xy축의 데카르트 좌표계를 생각하면 된다. (서풍은 서쪽에서 동쪽으로 부는 바람이다. 헷갈리지 말자!) 이 두 축의 속력을 각각 벡터의 인자로 간주하여 합성하면 방향과 양을 지닌 속도 벡터가 나오게 된다.

 

 

한 시점의 데이터를 사용할 생각이라면, 4번, 5번 열만 이용하면 된다. 미리 배열에 규칙을 정해놓고 데이터를 잘 배열한다면 2번, 3번 열은 없애버릴 수 있다.

 

struct WIND {
	float ugrd;
	float vgrd;
};

std::vector<vector<WIND>> windvec;

구조체는 위와 같이 잡고, std::vector 에 데이터를 넣었다. vector<WIND>에 한 시점의 경위도에 따른 벡터를 넣고, 다시 시간대별로 구분되는 vector들을 담는 vector를 정의했다. 그래서 이중의 vector가 되었다.

 

 

데이터는 https://github.com/ben-strasser/fast-cpp-csv-parser 라이브러리를 이용해서 아래와 같이 읽었다.

io::CSVReader<5, io::trim_chars<' '>, io::no_quote_escape<'\t'>> table(fileName);
table.set_header("time", "latitude", "longitude", "ugrd10m", "vgrd10m");

char* timechr;
double lat, lon;
timeZero = getTimeT("2012-08-19T00:00:00") + 32400;

while (table.read_row(timechr, lat, lon, wind.ugrd, wind.vgrd))
{
    int time = getTimeT(timechr) - timeZero;
    int hID= time / (3600 * 3);
    int x = (int)round((lon * 10.0) / 5);
    int y = (int)round((lat + 90.0) * 10.0 / 5);
    int index = x + y * 720;

    windvec[hID][index] = wind;		
}

3시간 단위의 데이터이므로, 첫 시점의 hID가 0이 되도록 변환시켰다. 그 다음 3시간 후는 hID가 1이 된다.

getTimeT 함수는 "2019-09-06T00:00:00Z" 와 같은 형식의 문자열을 유닉스타임으로 바꿔주는 함수다. timeZero에는 이 시각화에서 사용하는 첫 시간대 유닉스타임 값이 들어가 있다. 구현한 함수의 특징에 따라 9시간 시차인 32400초를 더해주어 시간대를 제대로 맞추었다.

 

 

텍스쳐에 넣기 위해 위경도를 인덱스로 변환한다. 부동소수점을 신경써야하므로 먼저 round 함수로 반올림해주고, 그 후에 정수 부분을 끊어내서 인덱스값을 만든다. 

특별한 방식은 아니고, 본래 1차원으로 된 데이터를 2차원처럼 간주해서 아래처럼 인덱스를 부여하는 흔한 방식이다. 

세계 지도를 펼쳐놓은 형상에 그대로 대응된다고 보면 된다.

 

 

 

 

 

Texture 정의

 

데이터는 텍스쳐에 넣는다. 위의 포맷으로 집어넣은 텍스쳐를 그대로 밀어넣으면 된다. 좌측 하단이 0,0이다. 

glCreateTextures(GL_TEXTURE_2D, 1, &windBO.windTexture0);  //GLenum target, 생성할 개수, id
glTextureStorage2D(windBO.windTexture0, 	//텍스쳐 id
    1,			// 밉맵 레벨 개수
    GL_RG32F,		// internal format
    360 * 2, 180 * 2 + 1); 	// 가로세로 크기
glTextureParameteri(windBO.windTexture0, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
glTextureParameteri(windBO.windTexture0, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

셰이더 스토리지 버퍼에 넣을 수도 있지만 이 경우 네이티브 선형 보간을 이용할 수 없으므로, 이 방식을 지원하는 텍스쳐를 사용하기로 한다. 예를 들어, 현재 데이터가 0.5도 간격으로 있는데, 그 사이 어딘가의 풍속 데이터를 구하려면, 이중선형 보간 연산을 직접 해줘야 한다. 공식은 복잡하지 않지만, 어쨌든 연산 비용이 들어간다. 그런데 텍스쳐의 경우 매핑 이미지를 확대하더라도 이미지가 부드럽게 그려지도록 하드웨어에서 이 연산을 지원한다. 계산 방식은 근본적으로 같을테지만, 훨씬 더 빨라지게 된다.

 

텍스쳐의 internal format은 RG32F로 둔다. red와 green 의 두 채널을 32비트 부동소수점 정확도로 다루게 된다. 각각의 채널에는 바로 앞에서 정리한 가로세로 벡터값을 넣는다. 크기 설정에서 180 * 2 + 1을 한 이유는 명확하다. -90도에서 90도까지 0.5도 간격으로 늘어놓으면 총 361개가 되기 때문이다.

 

텍스쳐 파라미터에서 여러가지 설정을 해 줘야 한다. 일단 텍스쳐 확대 축소시 GL_LINEAR를 셋팅한다. 그래야 확대될 때 그 사이사이 값을 선형으로 보간하게 된다. 이 부분이 바로 셰이더 스토리지 버퍼 대신 텍스쳐 버퍼를 선택한 가장 중요한 이유다.

 

 

 

드로우 콜

int arr01 = 0;

glTextureSubImage2D(windBO.windTexture0,	//텍스쳐 id
		0,			//밉맵레벨
		0, 0,		//x&y 오프셋
		360 * 2, 180 * 2 + 1,	// width & height
		GL_RG,		// format
		GL_FLOAT,	// type
		&windvec[arr01][0]);


int xCnt = 720;
int yCnt = 360 + 1;

glLineWidth(0.5);
glUniformMatrix4fv(soWindBasic.pmv, 1, GL_FALSE, glm::value_ptr(pmvM));
glBindTextureUnit(3, windBO.windTexture0);
glDrawArrays(GL_POINTS, 0, xCnt * yCnt);

드로우 콜은 위 처럼 보낸다. 일단 앞에서 읽어들인 여러가지 날짜 중 가장 앞의 날짜 한 시점의 데이터를 텍스쳐로 밀어넣었다. 

vertex Array 를 사용하지 않고, 그냥 그려야할 개수만 정의한다. vertex shader에서 gl_VertexID를 x y 인덱스로 변환시킬 예정이다.

 

 

Shader 구성

 

vertex shader 는 아래처럼 구성한다.

out Vertex
{	
	vec2 pos;
	int index;
} vertex;

void main() {
	
	int index = gl_VertexID;

	float lon = (float(index % 720)) / 2.0f;
	float lat = (float(index / 720) / 2.0f - 90.0f);

	vertex.pos = vec2(lon, lat);
	vertex.index = index;
}

gl_VertexID는 0번부터 259,919(720*(360+1)= 259,920) 번까지 발생하게 되는데, 각각의 vertexID를 텍스쳐 데이터 각각에 대응할 수 있도록 위경도 체계로 바꿔준다. 위와 같은 변환을 거치면 lon은 0~359.5까지 0.5도 간격으로, lat은 -90.0부터 90.0까지 0.5도 간격으로 변환된다.

 

geometry shader의 첫 부분은 아래처럼 구성한다.

#version 430 core
#define M_PI 3.14159265358979323846
#define WGS84_EQUATORIAL_RADIUS 6378137.0
uniform mat4 pmv;

layout (points) in;
layout (triangle_strip) out;
layout (max_vertices = 3) out;

in Vertex
{	
	vec2 pos;
	int index;
} vertex[];

out GS
{
	vec4 color;		
	vec2 txcoord;
} gs;

layout (binding = 3) uniform sampler2D wdVel0;

pmv에는 projection x modelview 매트릭스를 곱한 결과를 받아온다. 흔히 사용하는 방식이다.

 

일단 처음에는 바람을 그리기 전에 데이터가 텍스쳐로 제대로 정의되어 있는지 확인하려고 한다. 삼각형 핀 모양을 그려볼 것이므로 max_vertices = 3정도로만 해 둔다.

밖에서 3번으로 바인딩시킨 텍스쳐를 wdVel0 변수로 받아온다.

 

 

그 다음에는, main 함수에서 사용할 서브 함수들을 미리 정의한다. 

vec3 geodeticToCartesian(vec2 lonlat, float metersElevation) {

	float cosLat = cos(lonlat.y * (M_PI / 180)); 
	float sinLat = sin(lonlat.y * (M_PI / 180)); 
	float cosLon = cos(lonlat.x * (M_PI / 180));
	float sinLon = sin(lonlat.x * (M_PI / 180));

	//모두 장반경으로 계산한 약식임
	float x = (WGS84_EQUATORIAL_RADIUS + metersElevation) * cosLat * cosLon;
	float y = (WGS84_EQUATORIAL_RADIUS + metersElevation) * cosLat * sinLon;
	float z = (WGS84_EQUATORIAL_RADIUS + metersElevation) * sinLat;	

	return vec3(x,y,z);
}

mat4 rotationMatrix(vec3 axis, float angle)
{
    axis = normalize(axis);
    float s = sin(angle);
    float c = cos(angle);
    float oc = 1.0 - c;
    
    return mat4(oc * axis.x * axis.x + c,           oc * axis.x * axis.y - axis.z * s,  oc * axis.z * axis.x + axis.y * s,  0.0,
                oc * axis.x * axis.y + axis.z * s,  oc * axis.y * axis.y + c,           oc * axis.y * axis.z - axis.x * s,  0.0,
                oc * axis.z * axis.x - axis.y * s,  oc * axis.y * axis.z + axis.x * s,  oc * axis.z * axis.z + c,           0.0,
                0.0,                                0.0,                                0.0,                                1.0);
}


vec2 getVelocity(vec2 lonlat)
{
	vec2 texCoord = vec2(lonlat.x / 360.0f, (lonlat.y + 90.0f)/180.0f);
	vec2 velocity = texture(wdVel0, texCoord).rg;
	return velocity;
}

geodeticToCartesian( )은 경위도좌표와 고도를 넣으면 구면체의 미터법 좌표로 바꿔주는 함수다. 지구 중심이 (0,0,0)이 된다.

 

rotationMatrix( ) 는 특정 축을 중심으로 설정한 각도만큼 회전시키는 매트릭스를 만드는 일반적인 함수다.

 

getVelocity( ) 는 위경도를 넣으면 텍스쳐에서 해당 위치의 풍속 벡터를 가져오는 함수다. texture 함수가 작동하면서 텍스쳐에 정의된 값들 사이를 찌르면, 하드웨어 레벨에서 보간된 사이 값을 가져오게 된다. 앞에서 GL_LINEAR 로 정의내린 부분과 관계된다.

 

메인 함수는 아래처럼 구성한다.

void main() 
{	
	float elev = 3000.0f;
	float strength = 1000; //이번 프레임에서 진행시킬 부분
    
    //vertex shader에서 받아온 부분을 다시 정리
	vec2 posWgs84 = vertex[0].pos;	
	float lon = posWgs84.x;
	float lat = posWgs84.y;

	//경위도값을 지구 구면체 미터법 좌표로 바꿔준다.
	vec3 posCartesian = geodeticToCartesian(posWgs84,elev);
    
    //속도를 가져와서 이번 프레임에서 진행될 정도를 계산한다.
	vec2 velocity = getVelocity(posWgs84);	
	vec2 posnow = velocity * strength; //기준점에서의 속도

	//기준점에서의 바람 위치를 회전시켜서 해당 경위도 자리로 가져다놓는다.
	mat4 rotLon = rotationMatrix(vec3(0,0,-1), posWgs84.x * M_PI / 180.0f);
	mat4 rotLat = rotationMatrix(vec3(0,1,0), posWgs84.y* M_PI / 180.0f);		
	vec4 posRot = rotLon * rotLat * vec4(WGS84_EQUATORIAL_RADIUS + elev,
				posnow,
				1.0);
                
	//바람 방향의 수직 방향 벡터를 구한다.
	vec3 windNormal = normalize(posRot.xyz - posCartesian );
	mat4 rotPL = rotationMatrix(posCartesian, M_PI/2);
	vec4 pl = rotPL * vec4(windNormal,1.0);

	//바람 벡터 크기는 일단 일괄적으로 동일하게 둔다.
	float speed = 18000;	

	gs.color = vec4(1,1,1,1);
	gs.txcoord = vec2(-0.5,0.5);
	gl_Position = pmv * vec4(posCartesian - (windNormal * speed) + (pl.xyz * speed * 0.1),1.0);
	EmitVertex();
	gs.txcoord = vec2(0.5,0.5);
	gl_Position = pmv * vec4(posCartesian - (windNormal * speed) - (pl.xyz * speed * 0.1),1.0);
	EmitVertex();
	gs.txcoord = vec2(0,0);
	gl_Position = pmv * vec4(posCartesian + (windNormal * speed),1.0);
	EmitVertex();	
	EndPrimitive();
	
} //main

쓰면서 약간 복잡한 것 같다는 느낌이 들었는데, 좀 더 단순한 방법이 있는지는 다시 생각해봐야겠다. 각각의 단계는 주석으로 달아놓았다.

동서남북 0.5도 간격의 위경도 좌표 중심으로 바람 방향을 가리키는 화살표 머리 같은 삼각형을 그리기 위한 코드다.

 

 

 

 

풍속 벡터 검증

 

 

이렇게 그려보면 다음과 같이 지구 구면체에 바람 방향을 그려볼 수 있다.

위경도 0.5도 간격으로 풍속 벡터가 제대로 표현되고 있다.

 

이제 방향 이외에 속력도 반영해보자. geometry shader에서 speed 부분을 아래와 같이 수정해보자.

//float speed = 18000;
float speed = length(velocity) * 3000;

 

이제 다시 그려보면 풍속이 반영되었다.

 

 

지구 전체를 둘러보자. 텍스쳐가 이어진 부분들이 잘 되었는지, 하나로 모이는 극 지방이 제대로 되었는지 살펴보자.

역시나 이상하다. 위의 캡쳐는 북극이다. 고위도 부분이 촘촘한 것은 당연한 이치인데, 위도 90도 지점, 그러니까 북극의 벡터가 정확히 하나로 보이지 않는다.

 

데이터에는 이상이 없다. 아래 그림처럼 서로 다른 경도의 위도90도 벡터들은, 경도의 변화에 따라 조금씩 방향을 틀고 있다. 하나하나 그려보면 720개의 벡터가 정확히 한 방향으로 일치해야 한다.

 

이는, 텍스쳐 매핑의 가장자리 방식의 정의 때문이다. 

glTextureParameteri(windBO.windTexture0, GL_TEXTURE_WRAP_S, GL_REPEAT); //initial value
glTextureParameteri(windBO.windTexture0, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER); //이렇게 끊어주어야 극지방이 에러 안난다.

텍스쳐 방향은 S, T 로 정의하는데, 각각, x,y(가로,세로) 방향을 말한다. 디폴트 값은 GL_REPEAT 인데, 가로 방향의 경우 문제가 없다. 경도 359.5도 다음에 다시 0도가 와야 그 사이값들이 보간되기 때문이다. 그런데, 위도 90도와 -90도가 만나면 문제가 생긴다. 실제로 이어지지 않았고, 90도의 값이 -90도 값에 영향을 받으면 안되기 때문이다.

 

그래서 이 부분은 GL_CLAMP_TO_BORDER로 딱 잘라 끊어주어야 한다. wrapping 방식은 아래 그림으로 차이를 쉽게 이해할 수 있다.

이미지 출처 : https://open.gl/textures&amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;nbsp; (여기도 원 출처는 아닌 듯)

 

이제 다시 그려보면 제대로 되었음을 확인할 수 있다.

북극에서 화살표 하나로 보이는 벡터는 사실 720개가 겹쳐진 결과다.

 

 

 

벡터 필드의 완성

 

이번에는 사이사이의 보간이 제대로 되는지 확인해보자.

위경도 0.125도 간격으로 깔아보자. 4배의 해상도가 된다. 

 

먼저 드로우 콜에서 point 개수를 다음과 같이 바꿔준다.

int xCnt = 720 * 4;
int yCnt = (360 * 4) + 1;

 

 

vertex shader의 인덱싱 부분도 이에 상응하도록 바꿔준다.

	int index = gl_VertexID;

	float lon = (float(index % (720*4))) / (2.0f*4);
	float lat = (float(index / (720*4)) / (2.0f*4) - 90.0f);

 

 

이제 그려보면 아래 그림처럼 보인다.

지구 구면 전체에 415만개의 벡터들이 그려지는 셈인데, 제대로 보간되는 것 같다.

여기까지로 해서, 말 그대로 "vector field"가 제대로 깔렸음을 확인한 셈이다.

 

 

 

 

랜덤하게 발생되는 개체의 정의

 

이제는 벡터 필드 위에 바람처럼 보이게 해 줄 개체들을 발생시켜보자.

 

기본 원리는 이렇다. 랜덤한 위치에서 개체를 그리고, 그 다음 프레임은 그 개체가 놓인 위치의 벡터 필드에서 꺼낸 풍속 벡터을 더해서 결정한다. 기울기 근사법같은 것을 사용하지는 않고, 그냥  [위치 + 속도 = 다음 위치] 의 간단한 식으로 계산한다. 이렇게 되면 같은 위치에 발생한 점들은 항상 같은 경로를 따라서 움직이지만, 랜덤하게 발생한 개체가 이동 후 소멸하는 장면들이 계속되다보면 마치 바람이 부는 것처럼 보이게 된다.

 

머리와 꼬리가 뾰족하고 몸통이 통통한 개체를 정의하면 아래 그림처럼 경로를 따라 흐르게 된다.

 

vertex shader는 다음과 같이 다시 써본다.

#version 430 core
#define M_PI 3.14159265358979323846
#define STATIC_DURATION 360
#define LIFETIME 30

uniform int timeSimul;

out Vertex
{	
	vec2 pos;
	int index;
} vertex;

float rand(float x, float y) { 
    float t = dot(vec2(12.9898, 78.233), vec2(x,y)); 
    return fract(sin(t) * (4375.85453 + t)); 
}

void main() {
	
	int index =gl_VertexID;

	//동일 개체 기준으로 볼 때 (STATIC_DURATION * LIFETIME) 초마다 값이 변화함
	int tt = timeSimul/STATIC_DURATION + index;
	float seedAdd = float(tt/LIFETIME)*0.00001;

	float seed0 = (float(index % 2880)) / 8.0f;
	float seed1 = (float(index / 2880) / 8.0f - 90.0f);
	
	float randReturn = rand(seed0+seedAdd, seed1+seedAdd);
	float lon = randReturn * 360.0f;
	float lat = acos(2.0*rand(seed1+seedAdd, seed0+seedAdd)-1.0) /M_PI * 180.0-90.0f;// 180.0f -90.0f;

	vertex.pos = vec2(lon, lat);
	vertex.index = index;

}

랜덤함수는 인터넷에서 구했다. 변수 두개를 넣어 하나의 난수를 받은 함수다. 함수 안에 적힌 숫자들은 아무런 숫자들처럼 보이지만, mod 연산과 소수를 고려해서 정교하게 선별된 수들이다. 임의로 바꿔보면 리턴 값이 랜덤하게 발생되지 않고 일부 구간에 치우쳐서 발생되는 것을 볼 수 있다.

 

외부에서 초단위의 유닉스 타임을 timeSimul 변수로 받아온다. 이 값을 임의로 정의한 STATIC_DURATION로 나누어준 값을 인덱스에 더해주게 되는데, STATIC_DURATION는 변화하는 유닉스 타임 값을 일정 시간동안 동일하게 유지해 줌으로써 결과적으로 개체의 속도를 정의한다. 숫자를 키우면 천천히 지나가고, 숫자를 줄이면 순식간에 나타났다가 사라지게 된다. 물론 이 부분은 geometry shader에서 STATIC_DURATION을 사용한 방식과 관련 있으므로 함께 보아야 이해할 수 있다.

 

다시 그 값(tt)를 LIFETIME으로 나누어 준 후 0.00001을 곱한 값을 결과적으로 seed로 사용한다.

(STATIC_DURATION * LIFETIME)에 해당하는 시간동안 하나의 개체가 나타났다가 사라지게 된다. 물론 이 경우도 뒤의 geometry shader와 연관되어 있다.

0.00001처럼 충분히 작은 수를 곱해주지 않고 곧바로 seed로 사용하면, 랜덤 함수를 거쳐도 같은 자리에 발생되는 점들이 많아져서 겹쳐지게 된다. 

 

아무런 수도 곱해주지 않으면 아래와 같다. 진한 하얀색은 수많은 결과값이 겹쳐진 결과다.

 

 

0.01을 곱해주면 아래와 같다. 여전히 겹친 곳들이 많다.

 

 

0.00001을 곱해준 값을 seed로 넣으면 이제 고르게 랜덤 함수가 작동한다.

 

 

앞서 벡터 필드를 그려볼 때는 gl_VertedID값을 곧바로 위경도 0.5도 간격의 위치로 결정하고 끝냈는데, 이번에는 그렇게 고르게 배열시킨 위치들을, 랜덤함수에 넣을 seed로 활용한다. 0.0~1.0 사이의 랜덤 함수의 결과를 받고 다시 그 수를 위 경도값으로 변환하는 부분이 float lon, float lat을 정의하는 부분이다.

 

위도의 경우 아래처럼 acos 함수를 다시 한번 거쳤다.

float lat = acos(2.0*rand(seed1+seedAdd, seed0+seedAdd)-1.0) /M_PI * 180.0-90.0f;

그냥 난수로 발생시켜 줄 경우, -90도~90도까지 고르게 분포하게 되는데, 이렇게 되면 지구 구면체의 생김새 상 북극과 남극 부분, 즉 경도 1도 사이가 좁아지는 구간에 발생값들이 상대적으로 많이 몰리게 된다. acos 함수를 거치면, 극지방 부분의 값들을 조금 더 표면적이 넓은 적도 근처로 모아줄 수 있게 된다. 

 

 

 

이제 지오메트리 셰이더를 보자.

 

#version 430 core
#define M_PI 3.14159265358979323846
#define WGS84_EQUATORIAL_RADIUS 6378137.0
#define STATIC_DURATION 360
#define LIFETIME 30
uniform mat4 pmv;


layout (points) in;
layout (triangle_strip) out;
layout (max_vertices = 100) out;

in Vertex
{	
	vec2 pos;
	int index;
} vertex[];

out GS
{
	vec4 color;		
	vec2 txcoord;
} gs;

layout (binding = 3) uniform sampler2D wdVel0;
uniform int timeSimul;

uniform vec3 c12 = vec3(0xff/255.0f, 0x00/255.0f, 0x80/255.0f);
uniform vec3 c11 = vec3(0xff/255.0f, 0x8c/255.0f, 0x00/255.0f);
uniform vec3 c10 = vec3(0x40/255.0f, 0xe0/255.0f, 0xd0/255.0f);


vec3 ColorOri3(float t)
{
	vec3 col;
	if (t>0.5) col = mix(c21, c22, (t-0.5)/0.5);
	else if (t>0) col = mix(c20, c21, (t)/0.5);
	else col = c10;

	return col;
}

vec2 getVelocity(vec2 lonlat)
{
	vec2 texCoord = vec2(lonlat.x / 360.0f, (lonlat.y + 90.0f)/180.0f);
	vec2 velocity = texture(wdVel0, texCoord).rg;//vec2(3.5, 3.5);//
	return velocity;
}


...일부 함수 생략. 앞서 기술함...

void main() 
{	
	float elev = 1000.0f;
	float convertedStrength = 2000.0f;
	//시작점
	vec2 oriWgs84 = vertex[0].pos;	
	float lon = oriWgs84.x;
	float lat = oriWgs84.y;
	float width = 40000.0f;


	//회전은 왼손법칙
	//(1,0,0)은 그리니치 0도 방향
	//(0,1,0)은 동경 90도 방향
	//(0,0,1)이 북극 방향

	//static_duration 동안 같은 진행 상태로 유지한다. 360 후에 progress가 1 증가한다.
	//즉 생멸 전체를 보려면 360 * 30 = 10800초 동안 seed가 동일하게 유지되어야 한다.
	int tt = timeSimul/STATIC_DURATION + vertex[0].index;
	int progress = tt%LIFETIME; //0~29

	vec2 posnow = vec2(0,0);	
	vec2 posWgs84 = oriWgs84;
	vec3 posCartesianPre = geodeticToCartesian(posWgs84,elev);
	vec2 velocity = vec2(0,0);
	
	for (int i=9 ; i< LIFETIME ; i++)
	{ 
		
		vec3 posCartesian = geodeticToCartesian(posWgs84,elev);
		
		if (i>=progress && i<progress+10)
		{
			float speed = length(velocity) /20.0f;


			vec3 posNormal = normalize(posCartesian - posCartesianPre);
			mat4 rotPL = rotationMatrix(posCartesian, M_PI/2);
			vec4 pl = rotPL * vec4(posNormal,1.0);

			//진행에 따른 두께
			float t = (i-max(10,progress))  /  float(min(progress-1,min(9,LIFETIME-progress-1)));//19.0f;

			//i에 영향받지 않는 두께.
			float widthT = float(min(progress-1,min(9,LIFETIME-progress-1)))/10.0f;

			gs.color = vec4(ColorOri3(speed), widthT);
			float w = pow(widthT,3) * (sin(t*M_PI)+0.001) * width * pow(speed,0.5);

			gs.txcoord = vec2(-w,w);
			gl_Position = pmv * vec4(posCartesian + pl.xyz * w,1.0);
			EmitVertex();
			gs.txcoord = vec2(w,w);
			gl_Position = pmv * vec4(posCartesian - pl.xyz * w,1.0);
			EmitVertex();

		}
		//원점 중심 업데이트 된 위치를 구한다.
		velocity = getVelocity(posWgs84);
		posnow = velocity * convertedStrength; //일단 원점에 속도를 더한다.

		//그 위치를 회전시킨다. 
		mat4 rotLon = rotationMatrix(vec3(0,0,-1), posWgs84.x * M_PI / 180.0f);
		mat4 rotLat = rotationMatrix(vec3(0,1,0), posWgs84.y* M_PI / 180.0f);		
		vec4 posRot = rotLon * rotLat * vec4(WGS84_EQUATORIAL_RADIUS + elev,
					posnow,
					1.0);
		vec3 posRotNormalized = normalize(posRot.xyz); 
		
		//회전된 벡터의 새로운 각도를 구한다.
		float lonNew = atan(posRotNormalized.y/posRotNormalized.x) / M_PI * 180.0f;
		float latNew = asin(posRotNormalized.z) / M_PI * 180.0f;
		if (posRotNormalized.x>0 && posRotNormalized.y<0) lonNew+=360.0f;
		else if (posRotNormalized.x<0) lonNew +=180.0f;
		else if (posRotNormalized.x==0.0f && posRotNormalized.y>0) lonNew = 90.0f;
		else if (posRotNormalized.x==0.0f && posRotNormalized.y<0) lonNew = 270.0f;
		posWgs84 = vec2(lonNew, latNew);
		posCartesianPre = posCartesian;

	}

	EndPrimitive();

	
} //main

앞서 설명한 geometry shader와의 차이점은, 타원형 개체의 발생과 소멸을 위해 루프를 돌리면서 값들을 계산한 부분이다. 회전 매트릭스로 회전을 시키는 등 3차원 상의 좌표 변환에 사용하거나 두께 등을 정의해주기 위해서 이런저런 자잘한 연산을 했다.

 

마지막 부분에서는 결과값을 다시 위도 경도로 변환해서  lonNew, latNew에 저장한 뒤 루프의 다음번 회차에 사용하였다. 그 뒤에 경도가 0~360의 범위를 지닐 수 있도록 조정해주었다.

 

이제 아래와 같이 벡터 필드 위에서 바람이 흐른다.

 

 

seed 값을 준비할 때 아래처럼 vertex[0].index ( = gl_VertexID)를 더해주는 까닭은, 시간상에서 각 개체들의 시작 시각들 다르게 해주기 위함이다.

int tt = timeSimul/STATIC_DURATION + vertex[0].index;

만약 vertex shader와 geometry shader에서 각각의 고유 index를 더해주지 않으면 아래처럼 모두 한꺼번에 나타났다가 모두 한꺼번에 사라지게 된다.

 

 

STATIC_DURATION * LIFETIME 만큼의 시간이 지나면 개체가 한번 발생했다가 소멸된다. 그 후에 새로운 개체가 발생되는데, 이미 시간에 따라 seed 값이 달라져 있으므로 다른 위치에 발생하게 된다.

 

이렇게 벡터 필드 위를 유영하는 불규칙한 다수의 개체들은 바람이 흐르는 것과 같은 결과를 만들어준다.

 

 

 

유튜브에 올린 영상은 다수의 시점의 벡터필드를 구성한 후 한 시점에서 앞 뒤 두 시점의 벡터필드 텍스쳐를 참조하여 그린 것이다. 시간 진행에 따라 바꿔주고, 시간 사이사이는 앞뒤시점의 데이터를 바탕으로 선형 보간을 해서 만들었다.

 

쓰고나니 실제로 그리는 부분이 좀 복잡하고 설명이 부족하다. 초급은 약간 벗어나는 내용이라 어쩔 수 없는 부분도 있다. 모든 내용을 다 설명할 수는 없기 때문이다.

기본적으로 공간 좌표의 행렬변환에 익숙하지 않은 사람들은 별도로 그 부분을 익혀야 한다. 셰이더에서 무언가를 만들어본 사람들이라면, 이런 개별 개체 중심의 셰이더 표현 방식에 익숙할 것이므로 그리 어렵지 않게 이해할 수 있을 것 같다.

 

 

계산비용 개선의 여지

 

사실 위처럼 지오메트리 셰이더에서 루프를 돌면서 다수의 점들을 발생시키는 방식은 그리 좋지는 못하다. 물론 개체가 cuda 코어 수를 넘어가는 순간 그게 그거인것 같기도 하지만... 다른 방식으로 실험해보지는 않았다. 어쨌든 위에서는 그리지 않는 부분도 루프를 돌면서 위치를 순차적으로 계산해주고 있다. 그게 좀 비효율적인 것 같다.

 

이 부분은 좀 개선해볼 수 있을 것 같은데, 개체 하나당 30개 슬롯(위의 LIFETIME만큼)의 line strip 버퍼를 마련해놓고, 계산한 만큼은 저장해놓고, 다음 프레임에서 업데이트 되는 위치만 계산해주는 방식이다. 물론 이 경우도 메모리를 30배 읽어야 하는 단점이 있는데, 메모리를 읽는게 빠를지, 루프를 쓸데없이 돌더라도 한번 읽고 나머지 모두를 계산해주는게 빠를지는 비교해봐야 알 수 있다. 생각보다 메모리 읽는 비용이 매우 비싸기 때문이다.

 

끝.