BookmarkSubscribeRSS Feed
Choose Language Hide Translation Bar
Craige_Hales

Staff

Joined:

Mar 21, 2013

Wind Visualization

Making a visualization of the wind, similar to some that have appeared online for every hurricane recently.

JSL below with the parameters turned down for better interactive performance.

The model the code is using for the wind field is designed to make a pretty visualization and does not represent reality.

The vectors for the video look like this; the visualization changes when you drag the vectors.

 

vectors.PNGThe script uses Mouse Trap to position the arrows

 

 

NTRAIL = 50; // <<< try 3 to 50 or so -- longer trail of fading dots
NPARTICLES = 1000; // try 100...10000. interacts with NBINS, need less bins with less particles
SLOWER=50; // try 20 to 200 or so, smaller divisor makes bigger steps -- interacts with smoothing that affects cornering ability
SMOOTHING = 9; // try 1..10 or so to smooth the corners; smaller is more accurate, larger is prettier
NBINS=15; // try 2 to 19 or so. determines the granularity of the spawning locations where low densities are found
WINDFIELDSIZE = 40; // granularity of field, interacts with smoothing
PARTICLESIZE = 1; // try 0..10 or so may affect rendering speed
FADEPOWER = 1.5; // power tweaks fade curve
FADEMIN=20/255; // bottom end of fade curve 1..FADEMIN - darkest color that displays
graphxsize=1920/2;
graphysize=1080/2;
//outdir="F:\weather\windfieldpics\";
//==================

// flow visualization -- make a wind field
dx = J( WINDFIELDSIZE, WINDFIELDSIZE, 0 );
dy = J( WINDFIELDSIZE, WINDFIELDSIZE, 0 );

Parallel Assign( {SLOWER=SLOWER}, dx[r( -1, 1 ), c( -1, 1 )] = Sin( r ) / SLOWER ); 
Parallel Assign( {SLOWER=SLOWER}, dy[r( -1, 1 ), c( -1, 1 )] = Sin( c ) / SLOWER ); // 

windlist = {{{0.171875, 0.0458333333333334}, {0.15125, -0.715}}, {{0.171875, -0.9075},
{0.763125, -0.348333333333333}}, {{-0.859375, -0.504166666666667}, {-0.1925, -0.9075
}}, {{-0.385, 0.889166666666667}, {-0.955625, 0.119166666666667}}, {{0.78375, -0.165
}, {0.268125, 0.128333333333333}}, {{-0.240625, 0.0458333333333334}, {-0.78375,
-0.22}}, {{0.86625, -0.0183333333333333}, {0.2475, 0.8525}}, {{-0.144375,
-0.705833333333333}, {-0.2475, -0.0641666666666667}}};
/*
{
    {{.10,.20},{.50,.30}},
    {{.20,.20},{.50,.40}},
    {{.30,.20},{.50,.50}},
    {{.40,.20},{.50,.60}},
    {{.50,.20},{.50,.70}},
    {{.60,.20},{.50,.80}},
    {{.70,.20},{.50,.90}}
};
*/
caprow=.;
capcol=.;
ReleaseFunction = Function( {x, y},
    {},
    caprow = capcol = .;
            	// put the wind vectors into the dx/dy matrices
	
    Parallel Assign(
        {SLOWER = SLOWER, windlist = windlist},
        dx[r( -1, 1 ), c( -1, 1 )] = (x = 0 ; totalscale = 0 ; For( iwind = 1, iwind <= N Items( windlist ), iwind += 1,
            cx = (windlist[iwind][1][1] + windlist[iwind][2][1]) / 2;
            cy = (windlist[iwind][1][2] + windlist[iwind][2][2]) / 2;
            dx = windlist[iwind][2][1] - windlist[iwind][1][1]; // positive toward arrow
            scale = 1 / (Sqrt( (c - cx) ^ 2 + (r - cy) ^ 2 ) + .001);
            totalscale += scale;
            x += dx * scale;
        ) ; x / totalscale / SLOWER)
    ); 

    Parallel Assign(
        {SLOWER = SLOWER, windlist = windlist},
        dy[r( -1, 1 ), c( -1, 1 )] = (y = 0 ; totalscale = 0 ; For( iwind = 1, iwind <= N Items( windlist ), iwind += 1,
            cx = (windlist[iwind][1][1] + windlist[iwind][2][1]) / 2;
            cy = (windlist[iwind][1][2] + windlist[iwind][2][2]) / 2;
            dy = windlist[iwind][2][2] - windlist[iwind][1][2]; // positive toward arrow
            scale = 1 / (Sqrt( (c - cx) ^ 2 + (r - cy) ^ 2 ) + .001);
            totalscale += scale;
            y += dy * scale;
        ) ; y / totalscale / SLOWER)
    ); 
            	
);
  
ReleaseFunction(0,0);
            
New Window( "Example",
    Graph Box(
    	X Scale( -1.1, 1.1 ),
		Y Scale( -1.1, 1.1 ),
        Pen Size( 4 );
        For( iwind = 1, iwind <= N Items( windlist ), iwind += 1,
            Arrow( windlist[iwind][1], windlist[iwind][2] )
        );
        Mousetrap(
            Function( {x, y},
                {irow, nearrow = ., icol, nearcol = ., dist, neardist = 1e9},
                If( Is Missing( caprow ),
                    For( irow = 1, irow <= N Items( windlist ), irow += 1,
                        For( icol = 1, icol <= 2, icol += 1,
                            dist = Sqrt( (windlist[irow][icol][1] - x) ^ 2 + (windlist[irow][icol][2] - y) ^ 2 );
                            If( dist < neardist,
                                nearrow = irow;
                                nearcol = icol;
                                neardist = dist;
                            );
                        )
                    );
                    If( !Is Missing( nearrow ),
                        caprow = nearrow;
                        capcol = nearcol;
                    );
                );
                If( !Is Missing( caprow ) & !Is Missing( capcol ),
                    windlist[caprow][capcol][1] = x;
                    windlist[caprow][capcol][2] = y;
                );
            )//move
        , //
            ReleaseFunction(x,y)//release
        );
    )
);

// ===============================

SPREAD=1/NBINS; // mu for random normal spawn location
density = J( NBINS, NBINS, 0 );
//density[1, 1] = 1; // first loc max needs something to work with

_DX_ = 1;
_DY_ = 2;
_COLOR_ = 3;
_TRAILIDX_ = 4;
_TRAIL_ = 5; // NTRAIL begins here
_X_ = _TRAIL_;
_Y_ = _TRAIL_ + 1;
OUTOFPLAY = 99; // well beyond the graph's edge
DFRACTION=SMOOTHING;
NFRACTION = DFRACTION-1;
tran = J( NTRAIL, 1, 0 );
Parallel Assign( {FADEPOWER=FADEPOWER,FADEMIN=FADEMIN^(1/FADEPOWER)}, tran[r( 1, FADEMIN ), c] = r ^ FADEPOWER );
_STRAIL_ = 2; // stride x y pairs
particles = J( NPARTICLES , _TRAIL_ - 1 + _STRAIL_ * NTRAIL, OUTOFPLAY ); // dx,dy,color,xy trail has alpha fade
particles[0, _DX_] = 0;
particles[0, _DY_] = 0;
particles[0, _TRAILIDX_] = NTRAIL; // used to fade in, fade out
ExpectedParticlesPerBin = NPARTICLES/nitems(density);
New Window( "Example",
	g = Graph Box(
		Frame Size( graphxsize, graphysize ),
		X Scale( -1.1, 1.1 ),
		Y Scale( -1.1, 1.1 ),
		suppressaxes,
		<<Marker Size( PARTICLESIZE ),
		<<backgroundcolor( RGB Color( 0, .1, .2 ) ), // contrast with "white", below
		minDenLocs = Loc( density == Min( density ) );
		minDenYLocs = Floor( (minDenLocs - 1) / N Cols( Density ) + 1 );
		minDenXLocs = minDenLocs - (minDenYLocs - 1) * N Cols( Density );
		yymin = .5;
		yymax = N Rows( Density ) + .5;
		minDenYLocs = 2 * ((minDenYLocs - yymin) / (yymax - yymin)) - 1;
		xxmin = .5;
		xxmax = N Cols( Density ) + .5;
		minDenXLocs = 2 * ((minDenXLocs - xxmin) / (xxmax - xxmin)) - 1;

		density[0, 0] = 0;
		For( i = 1, i <= N Rows( particles ), i += 1, 
			fadeIdx = particles[i, _TRAILIDX_];
			If(
				fadeIdx > NTRAIL, // point is in the fading sequence
					// shift old particles out
                    particles[i, (_TRAIL_ + (_STRAIL_ * 1)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL) - 1)] = particles[i, (_TRAIL_ + (_STRAIL_ * 0)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL - 1) - 1)];
					nn = (particles[i, _TRAILIDX_] -= 1); // counting time while previous point fades away
					If( nn == NTRAIL, // restart this completely faded point
						rr = Random Integer( 1, N Rows( minDenXLocs ) );
						tx = particles[i, _X_ + _STRAIL_ * (NTRAIL - 1)] = Max( -.999, Min( .999, Random Normal( minDenXLocs[rr], SPREAD ) ) ); // next to last
						ty = particles[i, _Y_ + _STRAIL_ * (NTRAIL - 1)] = Max( -.999, Min( .999, Random Normal( minDenYLocs[rr], SPREAD ) ) );
						// index to field data
						idx = 1 + Round( (N Cols( dx ) - 1) * ((tx + 1) / 2) ); // tx+1 is 0...2
						idy = 1 + Round( (N Rows( dx ) - 1) * ((ty + 1) / 2) ); // divide by 2 is 0...1 becomes 0..nr-1...round, +1 becomes 1...nr
						particles[i, _DX_] = dx[idy, idx]; // dx,dy is the field speed. init the particles DX,DY
						particles[i, _DY_] = dy[idy, idx]; // to the correct value
						// the values just stored at NTRAIL-2 will display, slightly bright, on this cycle
						// and move to NTRAIL-1 (dimmest) on next cycle
					);//
				, /* elseif */ fadeIdx >= 1, // point is in the start sequence, or started					
					If( fadeIdx > 1, // still fading in
						(particles[i, _TRAILIDX_] -= 1);// for next round
						particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx - 2)] = particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx-1)] + particles[i, _DX_::_DY_];
					,// else...sticks at 1
                        // shift old particles out
                        particles[i, (_TRAIL_ + (_STRAIL_ * 1)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL) - 1)] = particles[i, (_TRAIL_ + (_STRAIL_ * 0)) :: (_TRAIL_ + _STRAIL_ * (NTRAIL - 1) - 1)];
                        particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx - 1)] = particles[i, (_X_ :: _Y_) + _STRAIL_ * (fadeIdx-0)] + particles[i, _DX_::_DY_];
					);

					If(
						!(-1 < particles[i, _X_ + _STRAIL_ * (fadeIdx - 1)] < 1) | !(-1 < particles[i, _Y_ + _STRAIL_ * (fadeIdx - 1)] < 1) |
						Sqrt( particles[i, _DX_] ^ 2 + particles[i, _DY_] ^ 2 ) < .0001, 
						// fix up particles leaving graph
						// leaving because a 99,99 was incremented by dx/dy
						particles[i, _TRAILIDX_] = 2 * NTRAIL;// enter fade-out sequence
						particles[i, _X_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
						particles[i, _Y_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
						particles[i, _DX_] = 0;
						particles[i, _DY_] = 0;	//
					, // else 
						// index to field data
						idx = 1 + Round( (N Cols( dx ) - 1) * ((particles[i, _X_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
						idy = 1 + Round( (N Rows( dx ) - 1) * ((particles[i, _Y_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
						// new deltas based on nearest field delta
						particles[i, _DX_] = (NFRACTION * particles[i, _DX_] + dx[idy, idx]) / DFRACTION; 
						particles[i, _DY_] = (NFRACTION * particles[i, _DY_] + dy[idy, idx]) / DFRACTION; 
						// nearest density bin...increment
						yyy = 1 + Round( (N Rows( density ) - 1) * ((particles[i, _Y_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
						xxx = 1 + Round( (N Cols( density ) - 1) * ((particles[i, _X_ + _STRAIL_ * (fadeIdx - 1)] + 1) / 2) );
						q=((density[yyy, xxx] += 1))/ExpectedParticlesPerBin;
						if( q*randomuniform()>2, 
							particles[i, _TRAILIDX_] = 2 * NTRAIL;// enter fade-out sequence
                            particles[i, _X_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
                            particles[i, _Y_ + _STRAIL_ * (1 - 1)] = OUTOFPLAY;
                            particles[i, _DX_] = 0;
                            particles[i, _DY_] = 0;	//
						);
					);//
				, // else bad _TRAILIDX_
				Throw( "???" )
			);
		);
		For( i = NTRAIL, i >= 1, i -= 1,
			Transparency( tran[i] );
			Marker( Color State( "white" ), particles[0, _X_ + _STRAIL_ * (i - 1)],/*flip -*/ particles[0, _Y_ + _STRAIL_ * (i - 1)] );
		);
	)
);
wait(.01);
picno = 1e6;
while(1, g << inval<<updatewindow;wait(0);
//pic = newimage(((g[framebox(1)]<<getimage)<<getpixels)[3::1082,3::1922]);
//pic<<saveimage(outdir||char(picno+=1)||".png","png");

 );
2 Comments
Super User

 So cool and an excellent example for Parallel Assign!

Staff

Thanks! I forgot to remove the two Parallel assigns at the top after I added the interactive controller. The two that are inside the ReleaseFunction (runs when the mouse button is released after dragging the controller arrows) are making a weighted sum of the vectors for each location in a grid. That grid is used to change the direction of every particle on each update. The release function could probably be called from the MouseTrap's move function and make the interaction even better--the particles would change direction while the arrow is moving.

Article Tags
Labels