From de1987eba84f828fdfddeb393d7faf089eb63c5b Mon Sep 17 00:00:00 2001 From: Martin Poirier Date: Thu, 10 Jul 2008 21:04:29 +0000 Subject: [PATCH] First draft for simulated annealing optimization method (disabled because not ok yet) --- source/blender/src/autoarmature.c | 162 ++++++++++++++++++++++++++++-- 1 file changed, 154 insertions(+), 8 deletions(-) diff --git a/source/blender/src/autoarmature.c b/source/blender/src/autoarmature.c index 9af3d639894..0644b484575 100644 --- a/source/blender/src/autoarmature.c +++ b/source/blender/src/autoarmature.c @@ -47,6 +47,7 @@ #include "BLI_editVert.h" #include "BLI_ghash.h" #include "BLI_graph.h" +#include "BLI_rand.h" #include "BDR_editobject.h" @@ -591,7 +592,7 @@ static float calcMaximumDistance(ReebArcIterator *iter, float *vec0, float *vec1 if (v1_inpf > 0) { int j; - for (j = i0; j < i1; j++) + for (j = i0 + 1; j < i1 - 1; j++) { float dist; @@ -660,6 +661,8 @@ static float calcCost(ReebArcIterator *iter, RigEdge *e1, RigEdge *e2, float *ve return new_cost; } +#define MAX_COST 100 /* FIX ME */ + static void calcGradient(RigEdge *e1, RigEdge *e2, ReebArcIterator *iter, int index, int nb_joints, float *cost_cube, int *positions, float **vec_cache) { EmbedBucket *bucket = NULL; @@ -700,7 +703,7 @@ static void calcGradient(RigEdge *e1, RigEdge *e2, ReebArcIterator *iter, int in if (index + 1 < nb_joints && next_position == positions[index + 1]) { - cost_cube[index * 3 + 2] = 1; + cost_cube[index * 3 + 2] = MAX_COST; } else { @@ -708,7 +711,7 @@ static void calcGradient(RigEdge *e1, RigEdge *e2, ReebArcIterator *iter, int in if (bucket == NULL) { - cost_cube[index * 3 + 2] = 1; + cost_cube[index * 3 + 2] = MAX_COST; } else { @@ -722,7 +725,7 @@ static void calcGradient(RigEdge *e1, RigEdge *e2, ReebArcIterator *iter, int in if (index - 1 > -1 && next_position == positions[index - 1]) { - cost_cube[index * 3] = 1; + cost_cube[index * 3] = MAX_COST; } else { @@ -730,7 +733,7 @@ static void calcGradient(RigEdge *e1, RigEdge *e2, ReebArcIterator *iter, int in if (bucket == NULL) { - cost_cube[index * 3] = 1; + cost_cube[index * 3] = MAX_COST; } else { @@ -741,6 +744,73 @@ static void calcGradient(RigEdge *e1, RigEdge *e2, ReebArcIterator *iter, int in } } +static float probability(float delta_cost, float iterations) +{ + if (delta_cost < 0) + { + return 1; + } + else + { + float temperature = (1 - iterations); + return (float)exp(delta_cost) * temperature; + } +} + +static int neighbour(int nb_joints, float *cost_cube, int *moving_joint, int *moving_direction) +{ + int total = 0; + int chosen = 0; + int i; + + for (i = 0; i < nb_joints; i++) + { + if (cost_cube[i * 3] < MAX_COST) + { + total++; + } + + if (cost_cube[i * 3 + 2] < MAX_COST) + { + total++; + } + } + + if (total == 0) + { + return 0; + } + + chosen = (int)BLI_drand() * total; + + for (i = 0; i < nb_joints; i++) + { + if (cost_cube[i * 3] < MAX_COST) + { + if (chosen == 0) + { + *moving_joint = i; + *moving_direction = -1; + break; + } + chosen--; + } + + if (cost_cube[i * 3 + 2] < MAX_COST) + { + if (chosen == 0) + { + *moving_joint = i; + *moving_direction = 1; + break; + } + chosen--; + } + } + + return 1; +} + static void retargetArctoArcAggresive(RigArc *iarc) { ReebArcIterator iter; @@ -798,7 +868,7 @@ static void retargetArctoArcAggresive(RigArc *iarc) vec_cache[0] = node_start->p; vec_cache[nb_edges] = node_end->p; -#if 1 +#if 1 /* BRUTE FORCE */ while(1) { float cost = 0; @@ -958,12 +1028,88 @@ static void retargetArctoArcAggresive(RigArc *iarc) memcpy(best_positions, positions, sizeof(int) * nb_joints); } } -#else +#elif 1 /* SIMULATED ANNEALING */ + { + RigEdge *previous; + float *cost_cube; + int k, kmax; + + kmax = 100; + + /* [joint: index][position: -1, 0, +1] */ + cost_cube = MEM_callocN(sizeof(float) * 3 * nb_joints, "Cost Cube"); + + initArcIterator(&iter, earc, node_start); + + /* init vec_cache */ + for (i = 0; i < nb_joints; i++) + { + bucket = peekBucket(&iter, positions[i]); + vec_cache[i + 1] = bucket->p; + } + + /* init cost cube */ + for (previous = iarc->edges.first, edge = previous->next, i = 0; + edge; + previous = edge, edge = edge->next, i += 1) + { + calcGradient(previous, edge, &iter, i, nb_joints, cost_cube, positions, vec_cache); + } + + for (k = 0; k < kmax; k++) + { + int status; + int moving_joint = -1; + int move_direction = -1; + float delta_cost; + + printf("\r%i", k); + + status = neighbour(nb_joints, cost_cube, &moving_joint, &move_direction); + + if (status == 0) + { + break; + } + + delta_cost = cost_cube[moving_joint * 3 + (1 + move_direction)]; + + if (probability(delta_cost, (float)k / (float)kmax) > BLI_frand()) + { + /* update position */ + positions[moving_joint] += move_direction; + + /* update vector cache */ + bucket = peekBucket(&iter, positions[moving_joint]); + vec_cache[moving_joint + 1] = bucket->p; + + /* update cost cube */ + for (previous = iarc->edges.first, edge = previous->next, i = 0; + edge; + previous = edge, edge = edge->next, i += 1) + { + if (i == moving_joint - 1 || + i == moving_joint || + i == moving_joint + 1) + { + calcGradient(previous, edge, &iter, i, nb_joints, cost_cube, positions, vec_cache); + } + } + } + } + + printf("\n"); + + memcpy(best_positions, positions, sizeof(int) * nb_joints); + + MEM_freeN(cost_cube); + } +#else /* GRADIENT DESCENT*/ { RigEdge *previous; float *cost_cube; - /* [joint: index][position: -1, 0, +2] */ + /* [joint: index][position: -1, 0, +1] */ cost_cube = MEM_callocN(sizeof(float) * 3 * nb_joints, "Cost Cube"); initArcIterator(&iter, earc, node_start);