// count the number of detections in measurement process
int count_detections (void)
{
// set cases of measurement results
double mtype[4];
mtype [0] = 0.0;
mtype [1] = 0.5;
mtype [2] = mtype[1] + 0.3;
mtype [3] = mtype[2] + 0.2;
// cout << "check measurement type3 = " << mtype[3] << endl; // just to check
// select a measurement case
double mrn = uniform_random();
int type = 1;
for ( int k = 0; k < 3; k++ )
{
if ( mrn < mtype[k] )
{
type = k;
break;
}
}
return type;
}
// distance between measurement and prediction
double distance(CvMat* measurement, CvMat* prediction)
{
double distance2 = 0;
double differance = 0;
for (int u = 0; u < 2; u++)
{
differance = cvmGet(measurement,u,0) - cvmGet(prediction,u,0);
distance2 += differance * differance;
}
return sqrt(distance2);
}
// likelihood based on multivariative normal distribution (ref. 15p eqn(2.4))
double likelihood(CvMat *mean, CvMat *covariance, CvMat *sample) {
/*
struct particle
{
double weight; // weight of a particle
CvMat* loc_p = cvCreateMat(2, 1, CV_64FC1); // previously estimated position of a particle
CvMat* loc = cvCreateMat(2, 1, CV_64FC1); // currently estimated position of a particle
CvMat* velocity = cvCreateMat(2, 1, CV_64FC1); // estimated velocity of a particle
cvSub(loc, loc_p, velocity); // loc - loc_p -> velocity
};
*/
int main (int argc, char * const argv[]) {
srand(time(NULL));
int width = 400; // width of image window
int height = 400; // height of image window
CvMat* state = cvCreateMat(2, 1, CV_64FC1);
// declare particles
/* particle pb[N]; // N estimated particles
particle pp[N]; // N predicted particles
particle pu[N]; // temporary variables for updating particles
*/
CvMat* pb [N]; // estimated particles
CvMat* pp [N]; // predicted particles
CvMat* pu [N]; // temporary variables to update a particle
CvMat* v[N]; // estimated velocity of each particle
CvMat* vu[N]; // temporary varialbe to update the velocity
double w[N]; // weight of each particle
for (int n = 0; n < N; n++)
{
pb[n] = cvCreateMat(2, 1, CV_64FC1);
pp[n] = cvCreateMat(2, 1, CV_64FC1);
pu[n] = cvCreateMat(2, 1, CV_64FC1);
v[n] = cvCreateMat(2, 1, CV_64FC1);
vu[n] = cvCreateMat(2, 1, CV_64FC1);
}
// initialize particles and the state
for (int n = 0; n < N; n++)
{
w[n] = (double) 1 / (double) N; // equally weighted
for (int row=0; row < 2; row++)
{
cvmSet(pb[n], row, 0, 0.0);
cvmSet(v[n], row, 0, 15 * uniform_random());
cvmSet(state, row, 0, 0.0);
}
}
// set the system
CvMat* groundtruth = cvCreateMat(2, 1, CV_64FC1); // groundtruth of states
CvMat* measurement [3]; // measurement of states
for (int k = 0; k < 3; k++) // 3 types of measurement
{
measurement[k] = cvCreateMat(2, 1, CV_64FC1);
}
// evaluate measurements
double range = (double) width; // range to search measurements for each particle
// cout << "range of distances = " << range << endl;
int mselected;
for (int index = 0; index < count; index++)
{
double d = distance(measurement[index], pp[n]);
if ( d < range )
{
range = d;
mselected = index; // selected measurement
}
}
/// cout << "selected measurement # = " << mselected << endl;
like[n] = likelihood(measurement[mselected], measurement_noise, pp[n]);
/// cout << "likelihood of #" << n << " = " << like[n] << endl;
Lessons:
1. transition noise와 measurement noise, (그것들의 covariances) 그리고 각 입자의 초기 위치와 상태를 알맞게 설정하는 것이 관건임을 배웠다. 그것을 Tuning이라 부른다는 것도.
1-1. 코드에서 가정한 system에서는 특히 입자들의 초기 속도를 어떻게 주느냐에 따라 tracking의 성공 여부가 좌우된다.
2. 실제 상태는 등가속도 운동을 하는 비선형 시스템이나, 여기에서는 프레임 간의 운동을 등속으로 가정하여 선형 시스템으로 근사한 모델을 적용한 것이다.
2-1. 그러므로 여기에 Kalman filter를 적용하여 결과를 비교해 볼 수 있겠다. 이 경우, 3차원 Kalman gain을 계산해야 한다.
2-2. 분홍색 부분을 고쳐 비선형 모델로 만든 후 Particle filtering을 하면 결과가 더 좋지 않을까? Tuning도 더 쉬어지지 않을까?
3. 코드 좀 정리하자. 너무 지저분하다. ㅡㅡ;
4. 아니, 근데 영 헤매다가 갑자기 따라잡는 건 뭐지??? (아래 결과)
// count the number of detections in measurement process
int count_detections (void)
{
// set cases of measurement results
double mtype[4];
mtype [0] = 0.0;
mtype [1] = 0.5;
mtype [2] = mtype[1] + 0.3;
mtype [3] = mtype[2] + 0.2;
// cout << "check measurement type3 = " << mtype[3] << endl; // just to check
// select a measurement case
double mrn = uniform_random();
int type = 1;
for ( int k = 0; k < 3; k++ )
{
if ( mrn < mtype[k] )
{
type = k;
break;
}
}
return type;
}
// get measurements
int get_measurement (CvMat* measurement[], CvMat* measurement_noise, double x, double y)
{
// set measurement types
double c1 = 1.0, c2 = 4.0;
// measured point 1
cvmSet(measurement[0], 0, 0, x + (c1 * cvmGet(measurement_noise,0,0) * gaussian_random())); // x-value
cvmSet(measurement[0], 1, 0, y + (c1 * cvmGet(measurement_noise,1,1) * gaussian_random())); // y-value
// measured point 2
cvmSet(measurement[1], 0, 0, x + (c2 * cvmGet(measurement_noise,0,0) * gaussian_random())); // x-value
cvmSet(measurement[1], 1, 0, y + (c2 * cvmGet(measurement_noise,1,1) * gaussian_random())); // y-value
// measured point 3 // clutter noise
cvmSet(measurement[2], 0, 0, width*uniform_random()); // x-value
cvmSet(measurement[2], 1, 0, height*uniform_random()); // y-value
// count the number of measured points
int number = count_detections(); // number of detections
cout << "# of measured points = " << number << endl;
// get measurement
for (int index = 0; index < number; index++)
{
cout << "measurement #" << index << " : "
<< cvmGet(measurement[index],0,0) << " " << cvmGet(measurement[index],1,0) << endl;
// set the system
CvMat* state = cvCreateMat(2, 1, CV_64FC1); // state of the system to be estimated
CvMat* groundtruth = cvCreateMat(2, 1, CV_64FC1); // groundtruth of states
CvMat* measurement [3]; // measurement of states
for (int k = 0; k < 3; k++) // 3 types of measurement
{
measurement[k] = cvCreateMat(2, 1, CV_64FC1);
}
// declare particles
CvMat* pb [N]; // estimated particles
CvMat* pp [N]; // predicted particles
CvMat* pu [N]; // temporary variables to update a particle
CvMat* v[N]; // estimated velocity of each particle
CvMat* vu[N]; // temporary varialbe to update the velocity
double w[N]; // weight of each particle
for (int n = 0; n < N; n++)
{
pb[n] = cvCreateMat(2, 1, CV_64FC1);
pp[n] = cvCreateMat(2, 1, CV_64FC1);
pu[n] = cvCreateMat(2, 1, CV_64FC1);
v[n] = cvCreateMat(2, 1, CV_64FC1);
vu[n] = cvCreateMat(2, 1, CV_64FC1);
}
// initialize the state and particles
for (int n = 0; n < N; n++)
{
w[n] = (double) 1 / (double) N; // equally weighted
for (int row=0; row < 2; row++)
{
cvmSet(state, row, 0, 0.0);
cvmSet(pb[n], row, 0, 0.0);
cvmSet(v[n], row, 0, 15 * uniform_random());
}
}
// set the process noise
// covariance of Gaussian noise to control
CvMat* transition_noise = cvCreateMat(2, 2, CV_64FC1);
cvmSet(transition_noise, 0, 0, 3); //set transition_noise(0,0) to 3.0
cvmSet(transition_noise, 0, 1, 0.0);
cvmSet(transition_noise, 1, 0, 0.0);
cvmSet(transition_noise, 1, 1, 6);
// set the measurement noise
// covariance of Gaussian noise to measurement
CvMat* measurement_noise = cvCreateMat(2, 2, CV_64FC1);
cvmSet(measurement_noise, 0, 0, 2); //set measurement_noise(0,0) to 2.0
cvmSet(measurement_noise, 0, 1, 0.0);
cvmSet(measurement_noise, 1, 0, 0.0);
cvmSet(measurement_noise, 1, 1, 2);
// initialize the image window
cvZero(iplImg);
cvNamedWindow("ParticleFilter-3d", 0);
for (int t = 0; t < step; t++) // for "step" steps
{
// cvZero(iplImg);
cout << "step " << t << endl;
// get the groundtruth
double gx = 10 * t;
double gy = (-1.0 / width ) * (gx - width) * (gx - width) + height;
get_groundtruth(groundtruth, gx, gy);
// get measurements
int count = get_measurement(measurement, measurement_noise, gx, gy);
double like[N]; // likelihood between measurement and prediction
double like_sum = 0; // sum of likelihoods
for (int n = 0; n < N; n++) // for "N" particles
{
// predict
double prediction;
for (int row = 0; row < 2; row++)
{
prediction = cvmGet(pb[n],row,0) + cvmGet(v[n],row,0) + cvmGet(transition_noise,row,row) * gaussian_random();
cvmSet(pp[n], row, 0, prediction);
}
// cvLine(iplImg, cvPoint(cvRound(cvmGet(pp[n],0,0)), cvRound(cvmGet(pp[n],1,0))), cvPoint(cvRound(cvmGet(pb[n],0,0)), cvRound(cvmGet(pb[n],1,0))), CV_RGB(100,100,0), 1);
// cvCircle(iplImg, cvPoint(cvRound(cvmGet(pp[n],0,0)), cvRound(cvmGet(pp[n],1,0))), 1, CV_RGB(255, 255, 0));
// evaluate measurements
double range = (double) width; // range to search measurements for each particle
// cout << "range of distances = " << range << endl;
int mselected;
for (int index = 0; index < count; index++)
{
double d = distance(measurement[index], pp[n]);
if ( d < range )
{
range = d;
mselected = index; // selected measurement
}
}
// cout << "selected measurement # = " << mselected << endl;
like[n] = likelihood(measurement[mselected], measurement_noise, pp[n]);
// cout << "likelihood of #" << n << " = " << like[n] << endl;
like_sum += like[n];
}
// cout << "sum of likelihoods = " << like_sum << endl;
// estimate states
double state_x = 0.0;
double state_y = 0.0;
// estimate the state with predicted particles
for (int n = 0; n < N; n++) // for "N" particles
{
w[n] = like[n] / like_sum; // update normalized weights of particles
// cout << "w" << n << "= " << w[n] << " ";
state_x += w[n] * cvmGet(pp[n], 0, 0); // x-value
state_y += w[n] * cvmGet(pp[n], 1, 0); // y-value
}
cvmSet(state, 0, 0, state_x);
cvmSet(state, 1, 0, state_y);
// define integrated portions of each particles; 0 < a[0] < a[1] < a[2] = 1
a[0] = w[0];
for (int n = 1; n < N; n++)
{
a[n] = a[n - 1] + w[n];
// cout << "a" << n << "= " << a[n] << " ";
}
// cout << "a" << N << "= " << a[N] << " " << endl;
for (int n = 0; n < N; n++)
{
// select a particle from the distribution
int pselected;
for (int k = 0; k < N; k++)
{
if ( uniform_random() < a[k] )
{
pselected = k;
break;
}
}
// cout << "p " << n << " => " << pselected << " ";
Klein, G. and Murray, D. 2007. Parallel Tracking and Mapping for Small AR Workspaces
In Proceedings of the 2007 6th IEEE and ACM international Symposium on Mixed and Augmented Reality - Volume 00
(November 13 - 16, 2007). Symposium on Mixed and Augmented Reality.
IEEE Computer Society, Washington, DC, 1-10. DOI=
http://dx.doi.org/10.1109/ISMAR.2007.4538852
1. parallel threads of tracking and mapping
2. mapping from smaller keyframes: batch techniques (Bundle Adjustment)
3. Initializing the map from 5-point Algorithm
4. Initializing new points with epipolar search
5. mapping thousands of points
affine warp
warping matrix <- (1) back-projecting unit pixel displacements in the source keyframe pyramid level onto the patch's plane and then (2) projecting these into the current (target) frame
Rendering for an Interactive 360 Light Field Display
Andrew Jones Ian McDowall Hideshi Yamada Mark Bolas Paul Debevec
University of Southern California Institute for Creative Technologies Fakespace Labs Sony Corporation University of Southern California School of Cinematic Arts
(full 3D hemispherical scattering measurements ->)
1. The primary specular highlight continues all the way around the
hair, while the secondary highlight is confined to the side of the hair
toward the source.
2. A pair of large out-of-plane peaks, or glints, are present, and as
the incidence angle increases the peaks move closer to the incidence
plane, eventually merging and disappearing.
3. The scattering distribution depends on the angle of rotation of the
hair fiber about its axis. (Because hair fibers are not generally
circular in cross section.)
4. Three trasport modes are derived: surfacce reflection, transmission, and internal reflection.
2 Fibers 2.1 Hair fibers and fiber scattering
The fiber is modeled as a dielectric cylinder covered with tilted
sacles (the cuticle) and with a pigmented interior (the cortex).
The cones of the R and TRT components shift in opposite directions,
causing them to separate into two visually distinguishable highlights.
(The R highlight is white and the TRT highlight is colored.)
2.2 Scattering
The bidirectional scattering function S for a fiber (different from the bidirectional reflection distribution function f_r for a surface):
The scattering integral (the curve intensity scattered from an infinitesimal length of fiber):
The presence D in this equation (1) indicates that
a thick fiber intercepts more light, and therefore appears brighter from a distance, than a thin fiber.
3 Scattering measurements 3.1 Incidence plane
As the scattering angle increases, the secondary highlight fades out, while the primary highlight maintains more constant amplitude. Both peaks maintain approximately constant width.
The equal-angle peak
3.2 Normal plane
The hair has a 180 degree rotational symmetry and is bilaterally symmetric in cross section.
The evolution of the peaks as the fiber roatates appears similar to the internal reflection from a transparent elliptical cylinder.
3.3 3D hemispherical measurements 3.3.1 Changes in glints with angle of incidence
The azimuth at which the glints occur changes as a function of incidence angle, with the glints moving toward the incidence plane as the incidence moves from normal to grazing.
3.3.2 Hemispherical scattering
3.4 Summary
4 Theory of scattering from fibers 4.1 Scattering from cylinders
<A Simple, Efficient Method for Realistic Animation of Clouds>
Yoshinori Dobashi & Tsuyoshi Okita, Hiroshima City University Kazufumi Kaneda & Hideo Yamashita, Hiroshima University Tomoyuki Nishita, University of Tokyo
http://www.siggraph.org/s2006/main.php?f=conference&p=etech&s=morpho The line-shaped images of the digital projector are used as lighting
for the high-speed rotating model. The model's high-speed rotation is
synchronized with the lighting by a rotating polygon mirror that
reflects the line-shaped CG images of the digital projector. The type
of lighting is controlled by a touch panel.
Takashi Fukaya
NHK Science and Technical Research Laboratories
fukaya.t-ha (at) nhk.or.jp
Toshio Iwai
Media Artist
Yuko Yamanouchi
NHK Science and Technical Research Laboratories
Willian T. Reeves <Particle Systems—a Technique for Modeling a Class of Fuzzy Objects>, ACM Transactions on Graphics, Vol.2, no.2, April 1983, Page 91-108
The model definition is procedural and is controlled by random numbers.
efficiency of human design time (to obtain a highly detailed model)
ability to adjust the level of detail (to suit a specific set of viewing parameters)
fractal surfaces
It is easier to model "alive" objects changing form over a period of time.
keywords: image synthesis stochastic process Stochastics fractal surfaces procedure random numbers stochastic modeling fractal modeling
2. BASIC MODEL OF PARTICLE SYSTEMS
A particle system is a collection of many minute particles that together represent a fuzzy object. Over a period of time, particles are generated into a system, move and change from within the system, and die from the system.
initial position => the origin of a particle system initial velocity initial color <= average RGB values and the maximum deviation from them initial transparency initial size shape => a region of newly born random particles about its origin lifetime
A particle's initial color, transparency and size are determined by mean values like MeasSpeed, maximum variations like VarSpeed of below:
InitialSpeed = MeanSpeed + Rand()*VarSpeed
More complicated generation shapes based on the law of nature or on chaotic attractors have been envisioned. eg. streaked spherical shapes => motion-blur particles
2.3 Particle Dynamics 2.4 Particle Extinction
when a particle's lifetime reaches zero
when the intensity of a particle, calculated from its color and transparency, drops belowa specified threshold
when a particle moves more than a given distance in a given direction from the origin of its parent particle system
2.5 Particle Rendering (1) Explosions and fire, the two fuzzy objects we have worekd with the most, are modeled well with the assumption that each particle can be displayed as a point light source. (Other fuzzy objects, such as clouds and water, are not.) (2) Since particles do not reflect but emit light, shadows are no longer a problem. 2.6 Particle Hierarchy
The initial direction of the particles' movement was constrained by the
system's ejection angle to fall within the region bounded by the
inverted cone. As particles flew upward, the gravity parameter pulled them back down to the planet's surface, giving them a parabolic motion path. The number of particles generated per frame was based on the amount of screen area covered by the particle system. Varying the mean velocity parameter caused the explosions to be of different heights. The rate at which a particle's color changed simulated the cooling of a glowing piece of some hypothetical material.
When a motion picture camera is used to film live action at 24 frames per second, the camera shutter typically remains open for 1/50 of a second. The image captured on a frame is actually an integration of approximately half the motion that occurred between successive frames. An object moving quickly appears blurred in the individual still frames.
> Ideas of modeling objects as collections of particles Evans & Sutherland @http://www.es.com Roger Wilson at Ohio State Alvy Ray Smith and Jim Blinn, Cosmos series Alan Norton Jim Blinn, light reflection functions
> Ideas of volumetric representations Norm Badler and Joe O'Rourke, solid modeling, "bubbleman"
“A particle system is a collection of many many minute particles that together represent a fuzzy object. Over a period of time, particles are generated into a system, move and change from within the system, and die from the system.”
Alain Fournier (University of Toronto) & Don Fussell (The University of Texas at Austin)
& Loren Carpenter (Lucasfilm) <Computer Rendering of Stochastic Models>
Parallel Tracking and Mapping for Small AR Workspaces.pdf
