dense_tracking.cpp 67 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979
  1. /*
  2. * dense_tracking.cpp
  3. *
  4. * Created on: Mar 7, 2016
  5. * Author: Janai
  6. */
  7. #include "configuration.h"
  8. #include <iostream>
  9. #include <fstream>
  10. #include <stdio.h>
  11. #include <string>
  12. #include <cmath>
  13. #include <omp.h>
  14. #include <random>
  15. #include <set>
  16. #include <flann/flann.hpp>
  17. #include <chrono>
  18. #include <opencv2/core.hpp>
  19. #include <opencv2/highgui.hpp>
  20. #include <opencv2/imgproc.hpp>
  21. #include <boost/filesystem.hpp>
  22. #include <boost/multi_index_container.hpp>
  23. #include <boost/multi_index/indexed_by.hpp>
  24. #include <boost/multi_index/ordered_index.hpp>
  25. #include <boost/multi_index/composite_key.hpp>
  26. #include <boost/multi_index/member.hpp>
  27. #include <boost/chrono/thread_clock.hpp>
  28. #include <gsl/gsl_multifit.h>
  29. #include <gsl/gsl_randist.h>
  30. #include "epic_flow_extended/image.h"
  31. #include "epic_flow_extended/io.h"
  32. #include "epic_flow_extended/epic.h"
  33. #include "epic_flow_extended/variational_mt.h"
  34. #include "utils/utils.h"
  35. #include "utils/hypothesis.h"
  36. #include "utils/parameter_list.h"
  37. #include "penalty_functions/penalty_function.h"
  38. #include "penalty_functions/lorentzian.h"
  39. #include "penalty_functions/modified_l1_norm.h"
  40. #include "penalty_functions/quadratic_function.h"
  41. // include Hamilton-Adams demosaicing
  42. extern "C"
  43. {
  44. #ifdef DMGUNTURK
  45. #include DMGUNTURK_PATH(/dmha.h)
  46. #endif
  47. }
  48. // include flowcode (middlebury devkit)
  49. #include MIDDLEBURY_PATH(/colorcode.h)
  50. #include MIDDLEBURY_PATH(/flowIO.h)
  51. // include TRWS
  52. #include TRWS_PATH(/instances.h)
  53. #include TRWS_PATH(/MRFEnergy.h)
  54. void HADemosaicing(float *Output, const float *Input, int Width, int Height, int RedX, int RedY) {
  55. #ifdef DMGUNTURK
  56. HamiltonAdamsDemosaic(Output, Input, Width, Height, RedX, RedY); // Hamilton-Adams implemented by Pascal Getreuer
  57. #endif
  58. }
  59. using namespace std;
  60. using namespace cv;
  61. using namespace boost::multi_index;
  62. using namespace boost::chrono;
  63. void setDefaultVariational(ParameterList& params) {
  64. // general
  65. params.insert("verbose", "0", true);
  66. params.insert("threads", "1", true);
  67. params.insert("scale", "1.0f", true);
  68. params.insert("slow_flow_S", "2", true);
  69. // energy function
  70. params.insert("slow_flow_alpha", "4.0f");
  71. params.insert("slow_flow_gamma", "6.0f", true);
  72. params.insert("slow_flow_delta", "1.0f", true);
  73. // image pyramid
  74. params.insert("slow_flow_layers", "1", true);
  75. params.insert("slow_flow_p_scale", "0.9f", true);
  76. // optimization
  77. params.insert("slow_flow_niter_alter", "10", true);
  78. params.insert("slow_flow_niter_outer", "10", true);
  79. params.insert("slow_flow_thres_outer", "1e-5", true);
  80. params.insert("slow_flow_niter_inner", "1", true);
  81. params.insert("slow_flow_thres_inner", "1e-5", true);
  82. params.insert("slow_flow_niter_solver", "30", true);
  83. params.insert("slow_flow_sor_omega", "1.9f", true);
  84. // occlusion reasoning
  85. params.insert("slow_flow_occlusion_reasoning", "1", true);
  86. params.insert("slow_flow_occlusion_penalty", "0.1", true);
  87. params.insert("slow_flow_occlusion_alpha", "0.1", true);
  88. params.insert("slow_flow_output_occlusions", "1", true);
  89. // regularization
  90. params.insert("slow_flow_robust_color", "1", true);
  91. params.insert("slow_flow_robust_color_eps", "0.001", true);
  92. params.insert("slow_flow_robust_color_truncation", "0.5", true);
  93. params.insert("slow_flow_robust_reg", "1", true);
  94. params.insert("slow_flow_robust_reg_eps", "0.001", true);
  95. params.insert("slow_flow_robust_reg_truncation", "0.5", true);
  96. }
  97. void setDefault(ParameterList& params) {
  98. params.insert("verbose", "0", true);
  99. // dense tracking parameters
  100. params.insert("scale", "1"); // scale input images
  101. params.insert("acc_skip_pixel", "1"); // skip number of pixel for final resolution)
  102. params.insert("acc_occlusion", "0"); // set to 1 to use occlusions from jet optimization
  103. // data driven proposal generation
  104. params.insert("acc_consistency_threshold", "1.0"); // threshold for forward backward consistency
  105. params.insert("acc_discard_inconsistent", "1"); // set to 1 to discard inconsistent trajectories
  106. params.insert("acc_epic_interpolation", "1"); // set to 1 to use epic flow interpolation
  107. params.insert("acc_epic_skip", "2"); // skip pixel for epic flow interpolation
  108. // energy weights and penalites
  109. params.insert("acc_jet_consistency", "1.0"); // weight of flow data term
  110. params.insert("acc_brightness_constancy", "0.1"); // weight of brightness constancy
  111. params.insert("acc_gradient_constancy", "1.0"); // weight of gradient constancy
  112. params.insert("acc_occlusion_penalty", "500.0"); // penalty for occluded frames
  113. params.insert("acc_beta", "10.0"); // spatial smoothness flow
  114. params.insert("acc_satial_occ", "10.0"); // spatial smoothness occlusion
  115. params.insert("acc_temporal_occ", "10.0"); // temporal smoothness occlusion
  116. params.insert("acc_cv", "0.0"); // temporal flow smoothness
  117. params.insert("acc_traj_sim_method", "1"); // spatial trajectory smoothness compares 0: adjacent flow, 1: accumulated flow, 2: final flow
  118. params.insert("acc_traj_sim_thres", "0.1"); // threshold for distance between two hypotheses
  119. // occlusion initialization
  120. params.insert("acc_occlusion_threshold", "5.0"); // threshold for comparison of high speed flow and trajectory
  121. params.insert("acc_occlusion_fb_threshold", "5.0"); // threshold for forward backward consistency check
  122. // discrete optimization
  123. params.insert("acc_alternate", "5"); // alternate between trws and propagation
  124. params.insert("acc_approach", "0"); // 0: TRW-S, 1: BP
  125. params.insert("acc_trws_eps", "1e-5"); // stopping criterion for lower bound
  126. params.insert("acc_trws_max_iter", "10"); // maximum number of iterations
  127. // neighbor propagation
  128. params.insert("acc_neigh_hyp", "5"); // number of proposals propagated to neighbors
  129. params.insert("acc_neigh_hyp_radius", "100.0"); // radius for proposal propagation
  130. params.insert("acc_neigh_skip1", "2"); // skipped pixel for nearest neighbor search tree
  131. params.insert("acc_neigh_skip2", "4"); // skipped pixel for nearest neighbor search tree
  132. params.insert("acc_hyp_neigh_tryouts", "20"); // maximum number of tryouts for sampling
  133. // penalty
  134. params.insert("acc_penalty_fct_data", "1"); // penality fct data term 0: quadratic 1: modified l1 norm 2: lorentzian
  135. params.insert("acc_penalty_fct_data_eps", "0.001");
  136. params.insert("acc_penalty_fct_reg", "1"); // penality fct regularization term 0: quadratic 1: modified l1 norm 2: lorentzian
  137. params.insert("acc_penalty_fct_reg_eps", "0.001");
  138. }
  139. inline bool insideImg(double x, double y, int width, int height) {
  140. return (y >= 0 && y < height && x >= 0 && x < width);
  141. }
  142. bool compareHypotheses(hypothesis* h1, hypothesis* h2) {
  143. return (h1 != NULL && (h2 == NULL || h1->score() < h2->score())); // ascending order considering the score and null hypothesis to the end
  144. }
  145. float addJC(hypothesis* h, const Mat* obs, double acc_jc, double acc_cv, PenaltyFunction* phi_d, ParameterList& params, const Mat* occlusion_masks = NULL) {
  146. Point2d p = h->p;
  147. int height = obs[0].rows;
  148. int width = obs[0].cols;
  149. double jenergy = 0;
  150. double cvenergy = 0;
  151. int contribution = 0;
  152. for (uint32_t j = 0; j < params.Jets; j++) {
  153. double u_j = h->u(j);
  154. double v_j = h->v(j); // x, y
  155. double u_jm1 = 0;
  156. double v_jm1 = 0;
  157. if (j > 0) {
  158. u_jm1 = h->u((j - 1));
  159. v_jm1 = h->v((j - 1)); // x, y
  160. }
  161. // exclude unknown flow
  162. if (u_j > UNKNOWN_FLOW_THRESH || v_j > UNKNOWN_FLOW_THRESH)
  163. break;
  164. // ############################################### compare flow to jet estimation
  165. if (insideImg(p.x + u_jm1, p.y + v_jm1, width, height)) {
  166. if (h->occluded(j) == 1 || h->occluded(j+1) == 1) // flow from j to j+1
  167. continue;
  168. double I_x = bilinearInterp<double>(p.x + u_jm1, p.y + v_jm1, obs[j], 1);
  169. double I_y = bilinearInterp<double>(p.x + u_jm1, p.y + v_jm1, obs[j], 0); // x, y
  170. jenergy += 0.5 * phi_d->apply(((u_j - u_jm1 - I_x) * (u_j - u_jm1 - I_x) + (v_j - v_jm1 - I_y) * (v_j - v_jm1 - I_y)));
  171. contribution++;
  172. }
  173. double u_jp1 = 0;
  174. double v_jp1 = 0;
  175. if (j + 1 < params.Jets) {
  176. u_jp1 = h->u((j + 1));
  177. v_jp1 = h->v((j + 1)); // x, y
  178. }
  179. // ############################################### assume constant velocity
  180. double u_sq = 2 * u_j - u_jm1 - u_jp1;
  181. double v_sq = 2 * v_j - v_jm1 - v_jp1;
  182. u_sq *= u_sq;
  183. v_sq *= v_sq;
  184. cvenergy += sqrt(u_sq + v_sq);
  185. }
  186. if(contribution > 0) jenergy /= contribution;
  187. return acc_jc * jenergy + acc_cv * cvenergy;
  188. }
  189. double dt_warp_time = 0, dt_med_time = 0, dt_sum_time = 0;
  190. /*
  191. * ################################### CONSIDER ALL POSSIBLE PAIRWISE TERMS #########################################
  192. */
  193. float addBCGC(hypothesis* h, color_image_t const* const * obs, color_image_t const* const * dx, color_image_t const* const * dy, double acc_bc, double acc_gc, int skip, ParameterList& params, const Mat* occlusion_masks = NULL) {
  194. Point2d p = h->p;
  195. int r = 0.5f * (skip + 1);
  196. int height = obs[0]->height;
  197. int width = obs[0]->width;
  198. int stride = obs[0]->stride;
  199. double wenergy = 0;
  200. double neighs = 0;
  201. time_t dt_start, dt_end;
  202. // ############################# Mean brightness and gradient constancy #############################
  203. for (int off_x = (p.x - r); off_x <= (p.x + r); off_x++) {
  204. for (int off_y = (p.y - r); off_y <= (p.y + r); off_y++) {
  205. if (off_x < 0 || off_x >= width || off_y < 0 || off_y >= height)
  206. continue;
  207. // warp images
  208. uint32_t visible = 0;
  209. vector<vector<double> > I(3, vector<double>(params.Jets + 1, 0));
  210. vector<vector<double> > Ix(3, vector<double>(params.Jets + 1, 0));
  211. vector<vector<double> > Iy(3, vector<double>(params.Jets + 1, 0));
  212. time(&dt_start);
  213. for (uint32_t j = 0; j < params.Jets + 1; j++) {
  214. double x_j = off_x;
  215. double y_j = off_y;
  216. if (j == 0) {
  217. int idx = stride * y_j + x_j;
  218. I[0][j] = obs[j]->c3[idx];
  219. I[1][j] = obs[j]->c2[idx];
  220. I[2][j] = obs[j]->c1[idx];
  221. Ix[0][j] = dx[j]->c3[idx];
  222. Ix[1][j] = dx[j]->c2[idx];
  223. Ix[2][j] = dx[j]->c1[idx];
  224. Iy[0][j] = dy[j]->c3[idx];
  225. Iy[1][j] = dy[j]->c2[idx];
  226. Iy[2][j] = dy[j]->c1[idx];
  227. visible++;
  228. } else {
  229. x_j += h->u(j - 1);
  230. y_j += h->v(j - 1);
  231. // stop in the case the object will stay occluded!
  232. if (insideImg(x_j, y_j, width, height) && (occlusion_masks == NULL || occlusion_masks[j].at<uchar>(y_j, x_j) != 0)) {
  233. I[0][j] = bilinearInterp(x_j, y_j, obs[j]->c3, height, width, stride);
  234. I[1][j] = bilinearInterp(x_j, y_j, obs[j]->c2, height, width, stride);
  235. I[2][j] = bilinearInterp(x_j, y_j, obs[j]->c1, height, width, stride);
  236. Ix[0][j] = bilinearInterp(x_j, y_j, dx[j]->c3, height, width, stride);
  237. Ix[1][j] = bilinearInterp(x_j, y_j, dx[j]->c2, height, width, stride);
  238. Ix[2][j] = bilinearInterp(x_j, y_j, dx[j]->c1, height, width, stride);
  239. Iy[0][j] = bilinearInterp(x_j, y_j, dy[j]->c3, height, width, stride);
  240. Iy[1][j] = bilinearInterp(x_j, y_j, dy[j]->c2, height, width, stride);
  241. Iy[2][j] = bilinearInterp(x_j, y_j, dy[j]->c1, height, width, stride);
  242. visible++;
  243. }
  244. }
  245. }
  246. time(&dt_end);
  247. dt_warp_time += difftime(dt_end, dt_start);
  248. int contribution = 0;
  249. double e_p = 0;
  250. time(&dt_start);
  251. // compute data terms
  252. for (uint32_t i = 0; i < visible; i++) {
  253. for (uint32_t j = (i + 1); j < visible; j++) {
  254. double x_i = off_x;
  255. double y_i = off_y;
  256. if (i > 0) {
  257. x_i += h->u(i - 1);
  258. y_i += h->v(i - 1);
  259. }
  260. double x_j = off_x + h->u(j - 1);
  261. double y_j = off_y + h->v(j - 1); // x, y
  262. if (insideImg(x_i, y_i, width, height) && insideImg(x_j, y_j, width, height)) {
  263. if (h->occluded(i) == 1 || h->occluded(j) == 1)
  264. continue; // skip occluded jets
  265. e_p += acc_bc * 0.3334 * (fabs(I[0][i] - I[0][j]) + fabs(I[1][i] - I[1][j]) + fabs(I[2][i] - I[2][j]));
  266. e_p += acc_gc * 0.3334 * (fabs(Ix[0][i] - Ix[0][j]) + fabs(Ix[1][i] - Ix[1][j]) + fabs(Ix[2][i] - Ix[2][j]) + fabs(Iy[0][i] - Iy[0][j]) + fabs(Iy[1][i] - Iy[1][j]) + fabs(Iy[2][i] - Iy[2][j]));
  267. contribution++;
  268. }
  269. }
  270. }
  271. time(&dt_end);
  272. dt_sum_time += difftime(dt_end, dt_start);
  273. if(contribution > 0) e_p /= contribution;
  274. wenergy += e_p;
  275. neighs++;
  276. }
  277. }
  278. if(neighs > 0) wenergy /= neighs;
  279. return wenergy;
  280. }
  281. float addOC(hypothesis* h, double acc_occ, double acc_temporal_occ, ParameterList& params) {
  282. int occlusions = 0;
  283. int change = 0;
  284. for (uint32_t i = 0; i < params.Jets + 1; i++) {
  285. // ############################# prefer smaller number of occlusions #############################
  286. occlusions += h->occluded(i); // skip occluded jets
  287. // ############################# expect temporal consistancy #############################
  288. if (i < params.Jets && h->occluded(i) != h->occluded(i + 1))
  289. change++;
  290. }
  291. return acc_occ * occlusions + acc_temporal_occ * change;
  292. }
  293. void computeSmoothnessWeight(const color_image_t *im, image_t* lum, float coef, const convolution_t *deriv, float avg_1, float avg_2, float avg_3, float std_1, float std_2, float std_3, bool hbit) {
  294. int i;
  295. image_t *lum_x = image_new(im->width, im->height), *lum_y = image_new(im->width, im->height);
  296. // compute luminance
  297. v4sf *im1p = (v4sf*) im->c1, *im2p = (v4sf*) im->c2, *im3p = (v4sf*) im->c3, *lump = (v4sf*) lum->data;
  298. for (i = 0; i < im->height * im->stride / 4; i++) {
  299. if (hbit)
  300. *lump = (0.299f * ((*im1p) * std_1 + avg_1) + 0.587f * ((*im2p) * std_2 + avg_2) + 0.114f * ((*im3p) * std_3 + avg_3)) / 65535.0f; // channels normalized to 0..1
  301. else
  302. *lump = (0.299f * ((*im1p) * std_1 + avg_1) + 0.587f * ((*im2p) * std_2 + avg_2) + 0.114f * ((*im3p) * std_3 + avg_3)) / 255.0f; // channels normalized to 0..1
  303. lump += 1;
  304. im1p += 1;
  305. im2p += 1;
  306. im3p += 1;
  307. }
  308. // compute derivatives with five-point tencil
  309. convolve_horiz(lum_x, lum, deriv);
  310. convolve_vert(lum_y, lum, deriv);
  311. // compute lum norm
  312. lump = (v4sf*) lum->data;
  313. v4sf *lumxp = (v4sf*) lum_x->data, *lumyp = (v4sf*) lum_y->data;
  314. for (i = 0; i < lum->height * lum->stride / 4; i++) {
  315. *lump = -coef * __builtin_ia32_sqrtps((*lumxp) * (*lumxp) + (*lumyp) * (*lumyp));
  316. lump[0][0] = 0.5f * expf((float) lump[0][0]);
  317. lump[0][1] = 0.5f * expf((float) lump[0][1]);
  318. lump[0][2] = 0.5f * expf((float) lump[0][2]);
  319. lump[0][3] = 0.5f * expf((float) lump[0][3]);
  320. lump += 1;
  321. lumxp += 1;
  322. lumyp += 1;
  323. }
  324. image_delete(lum_x);
  325. image_delete(lum_y);
  326. }
  327. /* show usage information */
  328. void usage(){
  329. printf("usage:\n");
  330. printf(" ./dense_tracking [cfg] -select [estimation for one specific final pair] -resume\n");
  331. printf("\n");
  332. }
  333. int main(int argc, char **argv) {
  334. if (argc < 2) {
  335. usage();
  336. return -1;
  337. }
  338. /*
  339. * getting config file location
  340. */
  341. uint32_t selected = 0;
  342. uint32_t selected_end = 0;
  343. /*
  344. * read in parameters and print config
  345. */
  346. string file = string(argv[1]);
  347. if (boost::filesystem::exists(file)) {
  348. printf("using parameters %s\n", file.c_str());
  349. } else {
  350. usage();
  351. return -1;
  352. }
  353. ParameterList params;
  354. setDefault(params);
  355. params.read(file);
  356. bool resume_frame = false;
  357. int max_fps = params.parameter<int>("max_fps", "0"); // fps of original sequence
  358. int threads = params.parameter<int>("threads", "1");
  359. #define isarg(key) !strcmp(a,key)
  360. if (argc > 1) {
  361. int current_arg = 1;
  362. while (current_arg < argc) {
  363. const char* a = argv[current_arg++];
  364. if (a[0] != '-') {
  365. continue;
  366. }
  367. if ( isarg("-h") || isarg("-help"))
  368. usage();
  369. else if (isarg("-output"))
  370. params.output = string(argv[current_arg++]);
  371. else if (isarg("-threads"))
  372. threads = atoi(argv[current_arg++]);
  373. else if( isarg("-resume") )
  374. resume_frame = true;
  375. else if( isarg("-select") ) {
  376. selected = atoi(argv[current_arg++]);
  377. selected_end = selected + 1;
  378. } else {
  379. fprintf(stderr, "unknown argument %s", a);
  380. usage();
  381. exit(1);
  382. }
  383. }
  384. } else {
  385. usage();
  386. return -1;
  387. }
  388. // check path
  389. for (uint32_t i = 0; i < params.jet_estimation.size(); i++)
  390. if (params.jet_estimation[i].back() != '/') params.jet_estimation[i] += "/";
  391. bool sintel = params.parameter<bool>("sintel", "0");
  392. bool subframes = params.parameter<bool>("subframes", "0"); // are subframes specified
  393. int skip_pixel = params.parameter<int>("acc_skip_pixel", "0");
  394. int ref_fps_F = params.parameter<int>("ref_fps_F", "1");
  395. uint32_t rates = params.jet_estimation.size();
  396. vector<float> weight_jet_estimation(rates, 0);
  397. for (uint32_t i = 0; i < rates; i++) {
  398. weight_jet_estimation[i] = i;
  399. if(params.jet_weight.size() > i)
  400. weight_jet_estimation[i] = params.jet_weight[i];
  401. }
  402. int min_fps_idx = params.parameter<int>("acc_min_fps", "0");
  403. /*
  404. * ############################### get step sizes ###############################
  405. */
  406. if (rates != params.jet_S.size()) {
  407. bool error = false;
  408. params.jet_S.resize(rates);
  409. for(uint32_t r = 0; r < rates; r++) {
  410. if(access((params.jet_estimation[r] + "/config.cfg").c_str(), F_OK) != -1) {
  411. ParameterList tmp((params.jet_estimation[r] + "/config.cfg"));
  412. if(tmp.exists("slow_flow_S"))
  413. params.jet_S[r] = tmp.parameter<int>("slow_flow_S");
  414. else {
  415. cerr << "Error reading S from " << params.jet_estimation[r] << "/config.cfg" << endl;
  416. error = true;
  417. }
  418. } else {
  419. cerr << "Error reading " << params.jet_estimation[r] << "/config.cfg" << endl;
  420. error = true;
  421. }
  422. }
  423. if(error) {
  424. cerr << "Frame rate or window size for Jet estimations missing!" << endl;
  425. exit(-1);
  426. }
  427. }
  428. // first slow flow estimation is used as reference
  429. int steps = params.jet_S[min_fps_idx] - 1;
  430. /*
  431. * get frame rates
  432. */
  433. if (rates <= 0) {
  434. cerr << "No Jet estimation specified!" << endl;
  435. exit(-1);
  436. }
  437. if (rates != params.jet_fps.size()) {
  438. bool error = false;
  439. params.jet_fps.resize(rates);
  440. for(uint32_t r = 0; r < rates; r++) {
  441. if(access((params.jet_estimation[r] + "/config.cfg").c_str(), F_OK) != -1) {
  442. ParameterList tmp((params.jet_estimation[r] + "/config.cfg"));
  443. if(tmp.exists("jet_fps")) {
  444. params.jet_fps[r] = tmp.parameter<int>("jet_fps");
  445. } else {
  446. cerr << "Error reading jet_fps from " << params.jet_estimation[r] << "/config.cfg" << endl;
  447. error = true;
  448. }
  449. } else {
  450. cerr << "Error reading " << params.jet_estimation[r] << "/config.cfg" << endl;
  451. error = true;
  452. }
  453. }
  454. if(error) {
  455. cerr << "Frame rate or window size for Jet estimations missing!" << endl;
  456. exit(-1);
  457. }
  458. }
  459. // ############################### specify reference framerate by maximum flow ########################################################
  460. params.Jets = params.jet_fps[min_fps_idx] / (1.0f * params.parameter<int>("ref_fps") * steps);
  461. /*
  462. * decompose sequence in subsets and adjust number of frames if necessary
  463. */
  464. uint32_t Jets = params.Jets;
  465. int skip = (1.0f * max_fps) / params.jet_fps[min_fps_idx]; // number of frames to skip for jet fps
  466. uint32_t gtFrames = params.Jets * skip;
  467. vector<int> compression_params;
  468. compression_params.push_back(CV_IMWRITE_PNG_COMPRESSION);
  469. compression_params.push_back(0);
  470. compression_params.push_back(CV_IMWRITE_JPEG_QUALITY);
  471. compression_params.push_back(100);
  472. // MAKE SURE FOLDER IS NOT OVERWRITTEN
  473. if(!resume_frame) {
  474. string newPath = params.output;
  475. if(newPath[newPath.length() - 1] == '/') newPath.erase(newPath.length() - 1);
  476. int num = 1;
  477. while(boost::filesystem::exists(newPath)) {
  478. cerr << newPath << " already exists!" << endl;
  479. stringstream tmp;
  480. tmp << newPath << "_" << num++;
  481. newPath = tmp.str();
  482. }
  483. params.output = newPath;
  484. }
  485. if(params.output[params.output.length() - 1] != '/') params.output = params.output + "/";
  486. boost::filesystem::create_directories(params.output); // accumulate folder
  487. ofstream infos;
  488. infos.open((params.output + "/config.cfg").c_str());
  489. infos << "# Slow Flow Accumulation\n";
  490. infos << params.cfgString(true);
  491. infos.close();
  492. cout << (params.output + "config.cfg").c_str() << endl;
  493. double traj_sim_thres = params.parameter<double>("acc_traj_sim_thres");
  494. int traj_sim_method = params.parameter<int>("acc_traj_sim_method");
  495. bool use_oracle = params.parameter<bool>("acc_oracle");
  496. bool use_occlusion = params.parameter<bool>("acc_occlusion");
  497. float acc_jc = params.parameter<float>("acc_jet_consistency");
  498. float acc_bc = params.parameter<float>("acc_brightness_constancy");
  499. float acc_gc = params.parameter<float>("acc_gradient_constancy");
  500. float acc_occ = params.parameter<float>("acc_occlusion_penalty");
  501. double acc_beta = params.parameter<double>("acc_beta"); // weight for flow spatial smoothness
  502. double acc_spatial_occ = params.parameter<double>("acc_spatial_occ"); // weight for occlusion spatial smoothness
  503. double acc_temporal_occ = params.parameter<double>("acc_temporal_occ"); // weight for occlusion spatial smoothness
  504. double acc_cv = params.parameter<double>("acc_cv"); // weight for constant velocity
  505. double outlier_beta = params.parameter<double>("acc_outlier_beta"); // energy for occlusion penalty
  506. float occlusion_threshold = params.parameter<float>("acc_occlusion_threshold"); // hypotheses from neighbors
  507. float occlusion_fb_threshold = params.parameter<float>("acc_occlusion_fb_threshold"); // hypotheses from neighbors
  508. // propagate hypotheses to neighbors
  509. int alternate = params.parameter<int>("acc_alternate"); // alternate between trws and perturbation
  510. uint32_t perturb_keep = params.parameter<int>("acc_perturb_keep"); // top x hypotheses will be kept
  511. bool discard_inconsistent = params.parameter<bool>("acc_discard_inconsistent");
  512. bool use_jet_occlusions = params.parameter<bool>("acc_use_jet_occlusions");
  513. uint32_t hyp_neigh = params.parameter<int>("acc_neigh_hyp"); // hypotheses from neighbors
  514. float hyp_neigh_radius = params.parameter<int>("acc_neigh_hyp_radius"); // radius to draw neighbors from
  515. bool draw_nn_from_radius = (hyp_neigh_radius > 0);
  516. float hyp_neigh_draws = params.parameter<int>("acc_neigh_draws"); // radius to draw neighbors from
  517. int hyp_neigh_tryouts = params.parameter<int>("acc_hyp_neigh_tryouts"); // maximum number of tryouts for sampling
  518. int nn_skip1 = params.parameter<int>("acc_neigh_skip1", "2");
  519. int nn_skip2 = params.parameter<int>("acc_neigh_skip2", "4");
  520. // discrete optimization robust penalty function
  521. int penalty_fct_data = params.parameter<int>("acc_penalty_fct_data"); // 0:QUADRATIC, 1:MOD_L1, 2:LORETZIAN
  522. double penalty_fct_data_eps = params.parameter<double>("acc_penalty_fct_data_eps");
  523. int penalty_fct_reg = params.parameter<int>("acc_penalty_fct_reg"); // 0:QUADRATIC, 1:MOD_L1, 2:LORETZIAN
  524. double penalty_fct_reg_eps = params.parameter<double>("acc_penalty_fct_reg_eps");
  525. unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
  526. if (params.exists("seed"))
  527. seed = params.parameter<int>("seed");
  528. default_random_engine generator(seed);
  529. // EPIC FLOW
  530. float epic_skip = params.parameter<float>("acc_epic_skip", "2");
  531. epic_params_t epic_params;
  532. epic_params_default(&epic_params);
  533. epic_params.pref_nn = 25;
  534. epic_params.nn = 160;
  535. epic_params.coef_kernel = 1.1f;
  536. // discrete optimization set options
  537. MRFEnergy<TypeGeneral>::Options options;
  538. options.m_eps = params.parameter<double>("acc_trws_eps");
  539. options.m_iterMax = params.parameter<int>("acc_trws_max_iter");
  540. options.m_printIter = 1;
  541. options.m_printMinIter = 0;
  542. // nearest neighbor search
  543. flann::SearchParams flann_params;
  544. flann_params.checks = 128; //
  545. flann_params.max_neighbors = -1; // unlimited
  546. flann_params.sorted = true;
  547. flann::SearchParams flann_params_p;
  548. flann_params_p.checks = 128; //
  549. flann_params_p.max_neighbors = -1; // unlimited
  550. flann_params_p.sorted = true;
  551. // discrete optimization set robust function
  552. PenaltyFunction *phi_d = NULL, *phi_s = NULL;
  553. switch (penalty_fct_data) {
  554. case 0:
  555. phi_d = new QuadraticFunction();
  556. break;
  557. case 1:
  558. phi_d = new ModifiedL1Norm(penalty_fct_data_eps);
  559. break;
  560. default:
  561. phi_d = new Lorentzian(penalty_fct_data_eps);
  562. break;
  563. }
  564. switch (penalty_fct_reg) {
  565. case 0:
  566. phi_s = new QuadraticFunction();
  567. break;
  568. case 1:
  569. phi_s = new ModifiedL1Norm(penalty_fct_reg_eps);
  570. break;
  571. default:
  572. phi_s = new Lorentzian(penalty_fct_reg_eps);
  573. break;
  574. }
  575. stringstream acc_folder;
  576. acc_folder << params.output;
  577. if(use_oracle)
  578. cout << endl << "Oracle uses same penalty as first jet estimation!!!!" << endl;
  579. boost::filesystem::create_directories(acc_folder.str());
  580. boost::filesystem::create_directories(acc_folder.str() + "gt_occlusions/"); // accumulate folder
  581. boost::filesystem::create_directories(acc_folder.str() + "occlusions/"); // accumulate folder
  582. boost::filesystem::create_directories(acc_folder.str() + "tmp/"); // accumulate folder
  583. params.insert("accumulation", acc_folder.str() + "frame_%i.flo", true);
  584. params.print();
  585. // in case of sintel data we would like to be able to distinguish frame number from 24 fps and 1008 fps
  586. if (sintel && !subframes)
  587. params.sequence_start = params.sequence_start * 1000; //
  588. if(selected_end == 0)
  589. selected_end = ref_fps_F;
  590. // iterate while changing start
  591. #pragma omp parallel for num_threads(threads) schedule(static,1)
  592. for(uint32_t start_jet = selected; start_jet < selected_end; start_jet++) {
  593. // alternate between continuous and discrete optimization
  594. stringstream numVariablesStream, factorsStream;
  595. double avg_unary_time = 0, avg_pairw_time = 0, avg_optimization_time = 0;
  596. // thread safe
  597. ParameterList thread_params;
  598. #pragma omp critical
  599. {
  600. thread_params = ParameterList(params);
  601. }
  602. // update start jet
  603. thread_params.sequence_start = thread_params.sequence_start + start_jet * Jets * steps * skip;
  604. int start_format = (thread_params.file.find_last_of('/') + 1);
  605. int end_format = thread_params.file.length() - start_format;
  606. string sequence_path = thread_params.file.substr(0, start_format);
  607. string format = "frame_%i.png";
  608. string flow_format = thread_params.parameter<string>("flow_format", "frame_%i");
  609. int len_format = flow_format.find_last_of('.');
  610. flow_format = flow_format.substr(0, len_format);
  611. format = thread_params.file.substr(start_format, end_format);
  612. if (sequence_path[sequence_path.length() - 1] != '/')
  613. sequence_path = sequence_path + "/";
  614. char final_file[1024];
  615. if (!sintel)
  616. sprintf(final_file, (acc_folder.str() + "/" + flow_format + ".flo").c_str(), thread_params.sequence_start);
  617. else
  618. sprintf(final_file, (acc_folder.str() + "/s" + flow_format + ".flo").c_str(), thread_params.sequence_start, 0);
  619. if (access(final_file, F_OK) != -1) {
  620. cout << "Flow file " << final_file << " already exists!" << endl;
  621. continue;
  622. }
  623. cout << "Flow file " << final_file << " does not exists!" << endl;
  624. /*
  625. * ################### read in image sequence ###################
  626. */
  627. vector<int> red_loc = thread_params.splitParameter<int>("raw_red_loc", "0,0");
  628. // read in sequence, forward and backward flow, and segmentation
  629. Mat *sequence = new Mat[Jets + 1]; // high speed video sequence
  630. Mat reference;
  631. color_image_t **data = new color_image_t*[Jets + 1]; // high speed video sequence
  632. color_image_t **data_dx = new color_image_t*[Jets + 1]; // derivative of high speed video sequence
  633. color_image_t **data_dy = new color_image_t*[Jets + 1]; // derivative of high speed video sequence
  634. Mat* gt = NULL;
  635. image_t*** gt_sfr = NULL;
  636. Mat* gt_occlusions = NULL;
  637. Mat* occlusions = new Mat[Jets];
  638. Mat *forward_flow = new Mat[Jets];
  639. Mat *backward_flow = new Mat[Jets];
  640. /*
  641. * ################### read in image sequence ###################
  642. */
  643. float norm = 1;
  644. for (uint32_t f = 0; f < (Jets + 1); f++) {
  645. char img_file[1024];
  646. if (!sintel) {
  647. sprintf(img_file, (sequence_path + format).c_str(), thread_params.sequence_start + f * steps * skip);
  648. } else {
  649. int sintel_frame = thread_params.sequence_start / 1000;
  650. int hfr_frame = f * steps * skip + (thread_params.sequence_start % 1000);
  651. while (hfr_frame < 0) {
  652. sintel_frame--;
  653. hfr_frame = 42 + hfr_frame;
  654. }
  655. while (hfr_frame > 41) {
  656. sintel_frame++;
  657. hfr_frame = hfr_frame - 42;
  658. }
  659. sprintf(img_file, (sequence_path + format).c_str(), sintel_frame, hfr_frame);
  660. }
  661. cout << "Reading " << img_file << "..." << endl;
  662. sequence[f] = imread(string(img_file), CV_LOAD_IMAGE_UNCHANGED); // load images
  663. if (sequence[f].type() == 2 || sequence[f].type() == 18)
  664. norm = 1.0f / 255; // for 16 bit images
  665. // convert to floating point
  666. sequence[f].convertTo(sequence[f], CV_32FC(sequence[f].channels()));
  667. /*
  668. * DEMOSAICING
  669. */
  670. if (thread_params.exists("raw") && thread_params.parameter<bool>("raw")) {
  671. Mat tmp = sequence[f].clone();
  672. color_image_t* tmp_in = color_image_new(sequence[f].cols, sequence[f].rows);
  673. color_image_t* tmp_out = color_image_new(sequence[f].cols, sequence[f].rows);
  674. switch (thread_params.parameter<int>("raw_demosaicing", "0")) {
  675. case 0: // use bilinear demosaicing
  676. sequence[f] = Mat::zeros(tmp.rows, tmp.cols, CV_32FC3);
  677. // bayer2rgb(tmp, sequence[f], green_start, blue_start); // red green
  678. bayer2rgbGR(tmp, sequence[f], red_loc[0], red_loc[1]); // red green
  679. break;
  680. case 1: // use hamilton adams demosaicing
  681. mat2colorImg<float>(sequence[f], tmp_in);
  682. HADemosaicing(tmp_out->c1, tmp_in->c1, tmp_in->width, tmp_in->height, red_loc[0], red_loc[1]); // Hamilton-Adams implemented by Pascal Getreuer
  683. sequence[f] = Mat::zeros(sequence[f].rows, sequence[f].cols, CV_32FC3);
  684. colorImg2colorMat<Vec3f>(tmp_out, sequence[f]);
  685. break;
  686. case 2: // use opencv demosaicing
  687. tmp.convertTo(tmp, CV_8UC1);
  688. sequence[f] = Mat::zeros(tmp.rows, tmp.cols, CV_8UC3);
  689. int code = CV_BayerBG2RGB;
  690. if (red_loc[1] == 0) // y
  691. if (red_loc[0] == 0) // x
  692. code = CV_BayerBG2RGB;
  693. else
  694. code = CV_BayerGB2RGB;
  695. else if (red_loc[0] == 0) // x
  696. code = CV_BayerGR2RGB;
  697. else
  698. code = CV_BayerRG2RGB;
  699. cv::cvtColor(tmp, sequence[f], code); // components from second row, second column !!!!!!!!!!!!!!!!!
  700. sequence[f].convertTo(sequence[f], CV_32FC(sequence[f].channels()));
  701. break;
  702. }
  703. color_image_delete(tmp_in);
  704. color_image_delete(tmp_out);
  705. } else {
  706. // covert to RGB
  707. cv::cvtColor(sequence[f], sequence[f], CV_BGR2RGB);
  708. }
  709. if (thread_params.parameter<bool>("grayscale", "0"))
  710. cvtColor(sequence[f], sequence[f], CV_RGB2GRAY);
  711. // use only a part of the images
  712. if (thread_params.extent.x > 0 || thread_params.extent.y > 0) {
  713. sequence[f] = sequence[f].rowRange(Range(thread_params.center.y - thread_params.extent.y / 2, thread_params.center.y + thread_params.extent.y / 2));
  714. sequence[f] = sequence[f].colRange(Range(thread_params.center.x - thread_params.extent.x / 2, thread_params.center.x + thread_params.extent.x / 2));
  715. }
  716. // rescale image with gaussian blur to avoid anti-aliasing
  717. double img_scale = thread_params.parameter<double>("scale", "1.0");
  718. if (img_scale != 1) {
  719. GaussianBlur(sequence[f], sequence[f], Size(0, 0), 1 / sqrt(2 * img_scale), 1 / sqrt(2 * img_scale), BORDER_REPLICATE);
  720. resize(sequence[f], sequence[f], Size(0, 0), img_scale, img_scale, INTER_LINEAR);
  721. }
  722. // print to file
  723. char file[1024];
  724. sprintf(file, (acc_folder.str() + "sequence/frame_%i.png").c_str(), thread_params.sequence_start - steps * skip + f * skip);
  725. Mat output_img;
  726. if (thread_params.parameter<bool>("16bit", "0")) {
  727. sequence[f].convertTo(output_img, CV_16UC(sequence[f].channels()));
  728. } else {
  729. sequence[f].convertTo(output_img, CV_8UC(sequence[f].channels()), norm);
  730. }
  731. if (thread_params.parameter<bool>("grayscale", "0"))
  732. cv::cvtColor(output_img, output_img, CV_GRAY2BGR); // OpenCV uses BGR
  733. else
  734. cv::cvtColor(output_img, output_img, CV_RGB2BGR); // OpenCV uses BGR
  735. if (thread_params.verbosity(WRITE_FILES)) {
  736. imwrite(file, output_img, compression_params);
  737. }
  738. data[f] = color_image_new(sequence[f].cols, sequence[f].rows);
  739. if (thread_params.parameter<bool>("grayscale", "0"))
  740. mat2colorImg<float>(sequence[f], data[f]);
  741. else
  742. colorMat2colorImg<Vec3f>(sequence[f], data[f]);
  743. }
  744. // normalize data terms
  745. normalize(data, Jets + 1, thread_params);
  746. // compute derivatives
  747. for (uint32_t f = 0; f < (uint32_t) (Jets + 1); f++) {
  748. data_dx[f] = color_image_new(sequence[f].cols, sequence[f].rows);
  749. data_dy[f] = color_image_new(sequence[f].cols, sequence[f].rows);
  750. float deriv_filter[3] = { 0.0f, -8.0f / 12.0f, 1.0f / 12.0f };
  751. convolution_t *deriv = convolution_new(2, deriv_filter, 0);
  752. color_image_convolve_hv(data_dx[f], data[f], deriv, NULL);
  753. color_image_convolve_hv(data_dy[f], data[f], NULL, deriv);
  754. convolution_delete(deriv);
  755. }
  756. // store reference image
  757. sequence[0].convertTo(reference, CV_8UC(sequence[0].channels()), norm);
  758. double img_scale = 1.0f / (skip_pixel + 1);
  759. if (img_scale != 1) {
  760. GaussianBlur(reference, reference, Size(0, 0), 1 / sqrt(2 * img_scale), 1 / sqrt(2 * img_scale), BORDER_REPLICATE);
  761. resize(reference, reference, Size(0, 0), img_scale, img_scale, INTER_LINEAR);
  762. }
  763. // reference image and edges for epic flow
  764. float_image forward_edges;
  765. color_image_t *imlab = NULL;
  766. image_t *epic_wx = NULL, *epic_wy = NULL;
  767. if(thread_params.parameter<bool>("acc_epic_interpolation", "1")) {
  768. char image_f[1000], edges_f[1000], edges_cmd[1000];
  769. sprintf(image_f, (acc_folder.str() + "tmp/frame_epic_%i.png").c_str(), thread_params.sequence_start);
  770. sprintf(edges_f, (acc_folder.str() + "tmp/edges_%i.dat").c_str(), thread_params.sequence_start);
  771. Mat tmp = reference.clone();
  772. // get lab image
  773. color_image_t* tmp_img = color_image_new(tmp.cols, tmp.rows);
  774. if(tmp.channels() == 1) {
  775. mat2colorImg<uchar>(tmp, tmp_img);
  776. } else
  777. colorMat2colorImg<Vec3b>(tmp, tmp_img);
  778. imlab = rgb_to_lab(tmp_img);
  779. color_image_delete(tmp_img);
  780. // get edges
  781. cv::cvtColor(tmp, tmp, CV_RGB2BGR); // OpenCV uses BGR
  782. imwrite(image_f, tmp, compression_params);
  783. sprintf(edges_cmd, "matlab -nodesktop -nojvm -r \"addpath(\'%s/matlab/\'); detect_edges(\'%s\',\'%s\'); exit\"", SOURCE_PATH.c_str(), image_f, edges_f);
  784. system(edges_cmd);
  785. forward_edges = read_edges(edges_f, tmp.cols, tmp.rows);
  786. }
  787. // ######### compute smoothness weighting
  788. float avg_1 = thread_params.parameter<double>("img_norm_avg_1", "0"), avg_2 = thread_params.parameter<double>("img_norm_avg_2", "0"), avg_3 = thread_params.parameter<double>("img_norm_avg_3", "0"), std_1 = thread_params.parameter<double>(
  789. "img_norm_std_1", "1"), std_2 = thread_params.parameter<double>("img_norm_std_2", "1"), std_3 = thread_params.parameter<double>("img_norm_std_3", "1");
  790. // use preset local smooth weights or compute
  791. convolution_t *deriv;
  792. float deriv_filter[3] = { 0.0f, -8.0f / 12.0f, 1.0f / 12.0f };
  793. deriv = convolution_new(2, deriv_filter, 0);
  794. image_t *smooth_weight = image_new(data[0]->width, data[0]->height);
  795. computeSmoothnessWeight(data[0], smooth_weight, 5.0, deriv, avg_1, avg_2, avg_3, std_1, std_2, std_3, thread_params.parameter<bool>("16bit", "0"));
  796. convolution_delete(deriv);
  797. /*
  798. * ################### read in ground truth ###################
  799. */
  800. if (use_oracle && thread_params.file_gt_list.size() > 0) {
  801. gt = new Mat[gtFrames];
  802. gt_sfr = new image_t**[1];
  803. for (uint32_t f = 0; f < gtFrames; f++) {
  804. char gtF[1024];
  805. if (!sintel)
  806. sprintf(gtF, thread_params.file_gt.c_str(), thread_params.sequence_start + f);
  807. else {
  808. int sintel_frame = thread_params.sequence_start / 1000;
  809. int hfr_frame = f + (thread_params.sequence_start % 1000);
  810. while (hfr_frame < 0) {
  811. sintel_frame--;
  812. hfr_frame = 42 + hfr_frame;
  813. }
  814. while (hfr_frame > 41) {
  815. sintel_frame++;
  816. hfr_frame = hfr_frame - 42;
  817. }
  818. sprintf(gtF, thread_params.file_gt.c_str(), sintel_frame, hfr_frame);
  819. }
  820. if (access(gtF, F_OK) == -1) {
  821. continue;
  822. cerr << "Error reading " << gtF << "!" << endl;
  823. }
  824. gt[f] = readGTMiddlebury(gtF);
  825. // rescale gt flow
  826. float rescale = (1.0f * sequence[0].cols) / gt[f].cols;
  827. resize(gt[f], gt[f], Size(0, 0), rescale, rescale, INTER_LINEAR);
  828. gt[f] = rescale * gt[f];
  829. }
  830. }
  831. /*
  832. * ################### read in ground truth occlusions ###################
  833. */
  834. if (use_oracle && thread_params.occlusions_list.size() > 0) {
  835. gt_occlusions = new Mat[gtFrames];
  836. for (u_int32_t f = 0; f < (uint32_t) (gtFrames); f++) {
  837. char oocF[1024];
  838. if (!sintel)
  839. sprintf(oocF, thread_params.occlusions_list[0].c_str(), thread_params.sequence_start + f);
  840. else {
  841. int sintel_frame = thread_params.sequence_start / 1000;
  842. int hfr_frame = f + (thread_params.sequence_start % 1000);
  843. while (hfr_frame < 0) {
  844. sintel_frame--;
  845. hfr_frame = 42 + hfr_frame;
  846. }
  847. while (hfr_frame > 41) {
  848. sintel_frame++;
  849. hfr_frame = hfr_frame - 42;
  850. }
  851. sprintf(oocF, thread_params.occlusions_list[0].c_str(), sintel_frame, hfr_frame);
  852. }
  853. // check if file exists
  854. if (access(oocF, F_OK) != -1) {
  855. gt_occlusions[f] = imread(string(oocF));
  856. float rescale = (1.0f * sequence[0].cols) / gt_occlusions[f].cols;
  857. resize(gt_occlusions[f], gt_occlusions[f], Size(0, 0), rescale, rescale, INTER_CUBIC);
  858. // use only a part of the images
  859. if (thread_params.extent.x > 0 || thread_params.extent.y > 0) {
  860. gt_occlusions[f] = gt_occlusions[f].rowRange(Range(thread_params.center.y - thread_params.extent.y / 2, thread_params.center.y + thread_params.extent.y / 2));
  861. gt_occlusions[f] = gt_occlusions[f].colRange(Range(thread_params.center.x - thread_params.extent.x / 2, thread_params.center.x + thread_params.extent.x / 2));
  862. }
  863. memset(oocF, 0, 1024);
  864. sprintf(oocF, "%s/gt_occlusions/occ_%05i.png", acc_folder.str().c_str(), thread_params.sequence_start + f);
  865. // write occlusion to file
  866. imwrite(oocF, gt_occlusions[f]);
  867. if (gt_occlusions[f].channels() > 1)
  868. cvtColor(gt_occlusions[f], gt_occlusions[f], CV_BGR2GRAY);
  869. gt_occlusions[f].convertTo(gt_occlusions[f], CV_8UC1);
  870. } else {
  871. cerr << "Error reading " << oocF << "!" << endl;
  872. }
  873. }
  874. }
  875. if (!use_occlusion) {
  876. delete[] occlusions;
  877. occlusions = NULL;
  878. }
  879. /*
  880. * ########################################### data driven proposal generation ###########################################
  881. */
  882. int owidth(sequence[0].cols), oheight(sequence[0].rows);
  883. int xy_incr = skip_pixel + 1;
  884. int xy_start = 0.5f * skip_pixel;
  885. int height = floor((1.0f * oheight) / xy_incr);
  886. int width = floor((1.0f * owidth) / xy_incr);
  887. int size = width * height;
  888. vector<vector<hypothesis*> > fb_hypotheses(size);
  889. // generate hypotheses from each jet estimation
  890. Mat consistent = Mat::zeros(height, width, CV_32SC1);
  891. for (uint32_t r = 0; r < rates; r++) {
  892. cout << "Processing " << thread_params.jet_estimation[r] << endl;
  893. int r_steps = thread_params.jet_S[r] - 1;
  894. float ratio = (1.0f * params.jet_fps[r]) / params.jet_fps[min_fps_idx];
  895. uint32_t r_Jets = (ratio * Jets);
  896. int r_skip = (1.0f * max_fps) / thread_params.jet_fps[r]; // number of frames to skip for jet fps
  897. Mat *r_forward_flow = new Mat[r_Jets];
  898. Mat *r_backward_flow = new Mat[r_Jets];
  899. /*
  900. * ################### read in flow fields ###################
  901. */
  902. for (uint32_t f = 0; f < r_Jets; f++) {
  903. char fFlowF[1024];
  904. char bFlowF[1024];
  905. sprintf(fFlowF, ("%s" + flow_format + ".flo").c_str(), thread_params.jet_estimation[r].c_str(), (thread_params.sequence_start + f * r_steps * r_skip));
  906. sprintf(bFlowF, ("%s" + flow_format + "_back.flo").c_str(), thread_params.jet_estimation[r].c_str(), (thread_params.sequence_start + f * r_steps * r_skip + r_steps * r_skip));
  907. if (!boost::filesystem::exists(fFlowF)) {
  908. cerr << fFlowF << " does not exist!" << endl;
  909. break;
  910. }
  911. if (!boost::filesystem::exists(bFlowF)) {
  912. cerr << bFlowF << " does not exist!" << endl;
  913. break;
  914. }
  915. r_forward_flow[f] = readGTMiddlebury(fFlowF);
  916. r_backward_flow[f] = readGTMiddlebury(bFlowF);
  917. // crop images
  918. if (thread_params.center.x > 0) {
  919. // NOTE: ROWRANGE IS INDUCING ERRORS IN ACCUMULATION!!!!!!
  920. r_forward_flow[f] = crop(r_forward_flow[f], thread_params.center, thread_params.extent);
  921. r_backward_flow[f] = crop(r_backward_flow[f], thread_params.center, thread_params.extent);
  922. }
  923. // rescale image
  924. float rescale = (1.0f * sequence[0].cols) / r_forward_flow[f].cols;
  925. resize(r_forward_flow[f], r_forward_flow[f], Size(0, 0), rescale, rescale, INTER_LINEAR);
  926. resize(r_backward_flow[f], r_backward_flow[f], Size(0, 0), rescale, rescale, INTER_LINEAR);
  927. r_forward_flow[f] *= rescale;
  928. r_backward_flow[f] *= rescale;
  929. if((int) r == min_fps_idx) {
  930. forward_flow[f] = r_forward_flow[f];
  931. backward_flow[f] = r_backward_flow[f];
  932. }
  933. }
  934. /*
  935. * ################### read in occlusions ###################
  936. */
  937. Mat* r_occlusions = NULL;
  938. if(use_jet_occlusions) {
  939. r_occlusions = new Mat[r_Jets];
  940. for (uint32_t f = 0; f < (r_Jets); f++) {
  941. char seqF[1024];
  942. sprintf(seqF, "%s/occlusion/frame_%i.pbm", thread_params.jet_estimation[r].c_str(), (thread_params.sequence_start + f * r_steps * r_skip));
  943. if (!boost::filesystem::exists(seqF)) {
  944. cerr << seqF << " does not exist!" << endl;
  945. break;
  946. }
  947. r_occlusions[f] = imread(string(seqF), CV_LOAD_IMAGE_UNCHANGED); // load images
  948. // crop images
  949. if (thread_params.center.x > 0)
  950. r_occlusions[f] = crop(r_occlusions[f], thread_params.center, thread_params.extent);
  951. // rescale image
  952. float rescale = (1.0f * sequence[0].cols) / r_occlusions[f].cols;
  953. resize(r_occlusions[f], r_occlusions[f], Size(0, 0), rescale, rescale, INTER_CUBIC);
  954. medianBlur ( r_occlusions[f], r_occlusions[f], 3 );
  955. if (thread_params.verbosity(WRITE_FILES)) {
  956. char file[1024];
  957. sprintf(file, (acc_folder.str() + "occlusion_%i.png").c_str(), thread_params.sequence_start + steps * skip + f * skip);
  958. imwrite(file, r_occlusions[f], compression_params);
  959. }
  960. // create mask out of the occlusion
  961. r_occlusions[f].convertTo(r_occlusions[f], CV_8UC1);
  962. r_occlusions[f] = (255 - r_occlusions[f]);
  963. if (thread_params.verbosity(VER_IN_GT)) {
  964. namedWindow("Occlusion"); // Create a window for display.
  965. imshow("Occlusion", r_occlusions[f]);
  966. waitKey(0);
  967. }
  968. if(min_fps_idx && (int) r == min_fps_idx) {
  969. occlusions[f] = r_occlusions[f];
  970. }
  971. }
  972. }
  973. /*
  974. * ################### accumulate consitent trajectories ###################
  975. */
  976. cout << "Generating hypotheses from consistent accumulations..." << endl;
  977. int created_hypotheses = 0;
  978. int rejected_hypotheses = 0;
  979. Mat r_consistent = Mat::zeros(height, width, CV_32SC1);
  980. int r_consistent_num = 0;
  981. Mat* acc_cons_flow = new Mat[r_Jets];
  982. Mat tracked;
  983. double threshold = thread_params.parameter<double>("acc_consistency_threshold");
  984. if(use_jet_occlusions)
  985. tracked = accumulateConsistentBatches(acc_cons_flow, r_forward_flow, r_backward_flow, r_occlusions, r_Jets, threshold, skip_pixel, discard_inconsistent, 0);
  986. else
  987. tracked = accumulateConsistentBatches(acc_cons_flow, r_forward_flow, r_backward_flow, NULL, r_Jets, threshold, skip_pixel, discard_inconsistent, 0);
  988. for (int y = 0; y < height; y++) {
  989. for (int x = 0; x < width; x++) {
  990. // only consider consistent trajectories
  991. if(tracked.at<int>(y,x) == (int) r_Jets) {
  992. consistent.at<int>(y,x) = 1;
  993. r_consistent.at<int>(y,x) = 1;
  994. r_consistent_num++;
  995. // get original pixel location
  996. int oy = y * xy_incr + xy_start;
  997. int ox = x * xy_incr + xy_start;
  998. // generate hypothesis
  999. double* f_x = new double[r_Jets];
  1000. double* f_y = new double[r_Jets];
  1001. for (int f = 0; f < tracked.at<int>(y,x); f++) {
  1002. f_y[f] = acc_cons_flow[f].at<Vec2d>(y, x)[0];
  1003. f_x[f] = acc_cons_flow[f].at<Vec2d>(y, x)[1];
  1004. }
  1005. hypothesis* hyp = new hypothesis(r_Jets, 0, tracked.at<int>(y,x), f_x, f_y, ox, oy);
  1006. fb_hypotheses[y * width + x].push_back(new hypothesis(hyp));
  1007. fb_hypotheses[y * width + x].back()->setLocation(ox, oy);
  1008. fb_hypotheses[y * width + x].back()->jet_est = r;
  1009. fb_hypotheses[y * width + x].back()->adaptFPS(Jets);
  1010. // compute energy
  1011. fb_hypotheses[y * width + x].back()->setOcclusions(forward_flow, backward_flow, occlusion_threshold, occlusion_fb_threshold);
  1012. fb_hypotheses[y * width + x].back()->energy = addJC(fb_hypotheses[y * width + x].back(), forward_flow, acc_jc, acc_cv, phi_d, thread_params, occlusions)
  1013. + addBCGC(fb_hypotheses[y * width + x].back(), data, data_dx, data_dy, acc_bc, acc_gc, skip_pixel, thread_params, occlusions)
  1014. + addOC(fb_hypotheses[y * width + x].back(), acc_occ, acc_temporal_occ, thread_params)
  1015. + weight_jet_estimation[r];
  1016. created_hypotheses++;
  1017. delete hyp;
  1018. } else
  1019. rejected_hypotheses++;
  1020. }
  1021. }
  1022. // obtain consistent maps
  1023. removeSmallSegments(r_consistent, 0.1, 100);
  1024. /*
  1025. * ################### epic flow interpolation in occluded regions ###################
  1026. */
  1027. if(thread_params.parameter<bool>("acc_epic_interpolation", "1")) {
  1028. vector<hypothesis*> epic_hypotheses(height * width, NULL);
  1029. // epic interpolation using consistent flow
  1030. for (uint32_t j = 0; j < r_Jets; j++) {
  1031. float_image matches = empty_image(float, 4, r_consistent_num);
  1032. int m = 0;
  1033. for (int y = floor(0.5f * epic_skip); y < height; y+=epic_skip) {
  1034. for (int x = floor(0.5f * epic_skip); x < width; x+=epic_skip) {
  1035. if(r_consistent.at<int>(y,x) == 1) {
  1036. matches.pixels[4*m] = x;
  1037. matches.pixels[4*m + 1] = y;
  1038. matches.pixels[4*m + 2] = x + acc_cons_flow[j].at<Vec2d>(y,x)[1] / (skip_pixel + 1);
  1039. matches.pixels[4*m + 3] = y + acc_cons_flow[j].at<Vec2d>(y,x)[0] / (skip_pixel + 1);
  1040. m++;
  1041. }
  1042. }
  1043. }
  1044. matches.ty = m;
  1045. cout << "Using " << m << " Matches!" << endl;
  1046. // prepare variables
  1047. image_t *wx = image_new(imlab->width, imlab->height), *wy = image_new(imlab->width, imlab->height);
  1048. epic(wx, wy, imlab, &matches, &forward_edges, &epic_params, 1);
  1049. if(thread_params.parameter<bool>("acc_epic_interpolation", "1")) {
  1050. for (int y = 0; y < height; y++) {
  1051. for (int x = 0; x < width; x++) {
  1052. int oy = y * xy_incr + xy_start;
  1053. int ox = x * xy_incr + xy_start;
  1054. // init epic hypotheses
  1055. if(epic_hypotheses[y * width + x] == NULL) {
  1056. epic_hypotheses[y * width + x] = new hypothesis(r_Jets, 0, r_Jets, new double[r_Jets], new double[r_Jets], ox, oy);
  1057. epic_hypotheses[y * width + x]->jet_est = r;
  1058. }
  1059. // set flow
  1060. epic_hypotheses[y * width + x]->flow_x[j] = wx->data[y * wx->stride + x] * (skip_pixel + 1);
  1061. epic_hypotheses[y * width + x]->flow_y[j] = wy->data[y * wy->stride + x] * (skip_pixel + 1);
  1062. // finish hypothesis
  1063. if(j == r_Jets - 1) {
  1064. epic_hypotheses[y * width + x]->adaptFPS(Jets);
  1065. epic_hypotheses[y * width + x]->setOcclusions(forward_flow, backward_flow, occlusion_threshold, occlusion_fb_threshold);
  1066. epic_hypotheses[y * width + x]->energy = addJC(epic_hypotheses[y * width + x], forward_flow, acc_jc, acc_cv, phi_d, thread_params, occlusions)
  1067. + addBCGC(epic_hypotheses[y * width + x], data, data_dx, data_dy, acc_bc, acc_gc, skip_pixel, thread_params, occlusions)
  1068. + addOC(epic_hypotheses[y * width + x], acc_occ, acc_temporal_occ, thread_params)
  1069. + weight_jet_estimation[r];
  1070. fb_hypotheses[y * width + x].push_back(epic_hypotheses[y * width + x]);
  1071. }
  1072. }
  1073. }
  1074. }
  1075. // write flow image to file
  1076. if(thread_params.verbosity(WRITE_FILES)) {
  1077. Mat floImg = flowColorImg(wx, wy, params.verbosity(VER_CMD));
  1078. if(floImg.data && !acc_folder.str().empty()) {
  1079. stringstream flowF;
  1080. flowF << acc_folder.str() << "tmp/epic_" << params.jet_fps[r] << "fps_" << thread_params.sequence_start << "_" << j << ".png";
  1081. imwrite((flowF.str()), floImg);
  1082. }
  1083. }
  1084. free(matches.pixels);
  1085. if((int) r == min_fps_idx && j == r_Jets - 1) {
  1086. epic_wx = wx;
  1087. epic_wy = wy;
  1088. image_mul_scalar(epic_wx, (skip_pixel + 1));
  1089. image_mul_scalar(epic_wy, (skip_pixel + 1));
  1090. } else {
  1091. image_delete(wx);
  1092. image_delete(wy);
  1093. }
  1094. }
  1095. }
  1096. if (thread_params.verbosity(VER_CMD))
  1097. cout << created_hypotheses << " trajectory hypotheses generated! (" << rejected_hypotheses << " rejected)" << endl;
  1098. // cleanup
  1099. delete[] acc_cons_flow;
  1100. delete[] r_forward_flow;
  1101. delete[] r_backward_flow;
  1102. if(r_occlusions != NULL) delete[] r_occlusions;
  1103. }
  1104. /*
  1105. * ###################### dense tracking formulation ######################
  1106. */
  1107. time_t unary_start, unary_end;
  1108. time_t opt_start, opt_end;
  1109. // reset data term time
  1110. dt_warp_time = 0;
  1111. dt_med_time = 0;
  1112. dt_sum_time = 0;
  1113. // initialize MRF object
  1114. uint32_t numVariables = height * width; // maximum number of variables
  1115. uint32_t factors = 0;
  1116. Mat selected_hyp = Mat::zeros(height, width, CV_32SC1);
  1117. /*
  1118. * ###################### discret optimization: flow reasoning ######################
  1119. */
  1120. for (int p_it = 0; p_it < alternate; p_it++) {
  1121. default_random_engine generator(seed + p_it);
  1122. if (p_it > 0) {
  1123. // keep only top x
  1124. for (int y = 0; y < height; y++) {
  1125. for (int x = 0; x < width; x++) {
  1126. int idx = y * width + x;
  1127. vector<hypothesis*> temp = fb_hypotheses[idx];
  1128. fb_hypotheses[idx].clear();
  1129. // add best hypotheses from last iteration
  1130. int last = selected_hyp.at<int>(y, x);
  1131. if (last >= 0 && (uint32_t) last < temp.size()) {
  1132. fb_hypotheses[idx].push_back(temp[last]);
  1133. temp[last] = NULL;
  1134. }
  1135. // sort in ascending order
  1136. sort(temp.begin(), temp.end(), compareHypotheses);
  1137. // add to get top x
  1138. uint32_t h = 0;
  1139. for (; h < temp.size(); h++) {
  1140. if(temp[h] == NULL)
  1141. continue;
  1142. if(fb_hypotheses[idx].size() <= perturb_keep)
  1143. fb_hypotheses[idx].push_back(temp[h]);
  1144. else
  1145. // clean up the rest
  1146. delete temp[h];
  1147. }
  1148. }
  1149. }
  1150. } else {
  1151. // SORT ACCORDING DATATERM!
  1152. for (int y = 0; y < height; y++) {
  1153. for (int x = 0; x < width; x++) {
  1154. int idx = y * width + x;
  1155. // sort all hypothesis besides outlier hypothesis
  1156. if (fb_hypotheses[idx].size() > 1) {
  1157. sort(fb_hypotheses[idx].begin(), fb_hypotheses[idx].end(), compareHypotheses);
  1158. }
  1159. }
  1160. }
  1161. }
  1162. /*
  1163. * ########################################### hypotheses from neighbors ###########################################
  1164. */
  1165. {
  1166. cout << "neighbors" << endl;
  1167. vector<vector<hypothesis*> > tmp = fb_hypotheses;;
  1168. clock_t start_nn = clock();
  1169. // build nearest neighbor tree
  1170. int counter2 = 0;
  1171. int counter4 = 0;
  1172. for (int y = 1; y < height; y++) {
  1173. for (int x = 1; x < width; x++) {
  1174. // only use consistent estimates for first draw
  1175. if(consistent.at<int>(y,x) == 1 || p_it > 0) {
  1176. if((y - 1) % nn_skip1 == 0 && (x - 1) % nn_skip1 == 0) counter2++;
  1177. if((y - 2) % nn_skip2 == 0 && (x - 2) % nn_skip2 == 0) counter4++;
  1178. }
  1179. }
  1180. }
  1181. vector<flann::Matrix<float>> dataset(2);
  1182. dataset[0] = flann::Matrix<float>(new float[counter2 * 2], counter2, 2);
  1183. dataset[1] = flann::Matrix<float>(new float[counter4 * 2], counter4, 2);
  1184. counter2 = 0;
  1185. counter4 = 0;
  1186. for (int y = 1; y < height; y++) {
  1187. for (int x = 1; x < width; x++) {
  1188. // only use consistent estimates for first draw
  1189. if(consistent.at<int>(y,x) == 1 || p_it > 0) {
  1190. if((y - 1) % nn_skip1 == 0 && (x - 1) % nn_skip1 == 0) {
  1191. dataset[0].ptr()[counter2 * 2] = x;
  1192. dataset[0].ptr()[counter2 * 2 + 1] = y;
  1193. counter2++;
  1194. }
  1195. if((y - 2) % nn_skip2 == 0 && (x - 2) % nn_skip2 == 0) {
  1196. dataset[1].ptr()[counter4 * 2] = x;
  1197. dataset[1].ptr()[counter4 * 2 + 1] = y;
  1198. counter4++;
  1199. }
  1200. }
  1201. }
  1202. }
  1203. vector<flann::Index<flann::L2<float> >*> nonempty;
  1204. nonempty.push_back(new flann::Index<flann::L2<float> >(flann::KDTreeSingleIndexParams()));
  1205. nonempty.push_back(new flann::Index<flann::L2<float> >(flann::KDTreeSingleIndexParams()));
  1206. nonempty[0]->buildIndex(dataset[0]);
  1207. nonempty[1]->buildIndex(dataset[1]);
  1208. cout << "First NN tree " << counter2 << endl << "Second NN tree " << counter4 << endl;
  1209. for (int y = 0; y < height; y++) {
  1210. for (int x = 0; x < width; x++) {
  1211. int i = y * width + x;
  1212. int oy = y * xy_incr + xy_start;
  1213. int ox = x * xy_incr + xy_start;
  1214. for(uint32_t t = 0; t < nonempty.size(); t++) {
  1215. // compare to existent trajectories (nearest neighbors)
  1216. flann::Matrix<float> query(new float[2], 1, 2);
  1217. std::vector<std::vector<int> > indices;
  1218. std::vector<std::vector<float> > dists;
  1219. std::vector<double> prob;
  1220. std::vector<double> cum;
  1221. query.ptr()[0] = x;
  1222. query.ptr()[1] = y;
  1223. int found_nn = 0;
  1224. if(draw_nn_from_radius) {
  1225. found_nn = nonempty[t]->radiusSearch(query, indices, dists, (t + 1) * (hyp_neigh_radius), flann_params_p);
  1226. if(found_nn < 50) {
  1227. found_nn = nonempty[t]->knnSearch(query, indices, dists, 50, flann_params_p);
  1228. }
  1229. } else {
  1230. found_nn = nonempty[t]->knnSearch(query, indices, dists, hyp_neigh_draws, flann_params_p);
  1231. }
  1232. // remove x,y
  1233. if(fb_hypotheses[i].size() > 0) {
  1234. std::vector<std::vector<int> > tmp_indices = indices;
  1235. std::vector<std::vector<float> > tmp_dists = dists;
  1236. indices[0].clear();
  1237. dists[0].clear();
  1238. int tmp_found_nn = 0;
  1239. for(int i = 0; i < found_nn; i++) {
  1240. if(dataset[t].ptr()[2 * tmp_indices[0][i]] != x && dataset[t].ptr()[2 * tmp_indices[0][i] + 2] != y) {
  1241. indices[0].push_back(tmp_indices[0][i]);
  1242. dists[0].push_back(tmp_dists[0][i]);
  1243. tmp_found_nn++;
  1244. }
  1245. }
  1246. found_nn = tmp_found_nn;
  1247. }
  1248. std::uniform_real_distribution<double> uniform1(0, 1);
  1249. uniform_int_distribution<int> uniform2(0, found_nn - 1);
  1250. int tryouts = 0;
  1251. // draw half of the neighbors
  1252. while(tryouts < hyp_neigh_tryouts && (fb_hypotheses[i].size() - tmp[i].size()) < (t + 1) * hyp_neigh) {
  1253. tryouts++;
  1254. int ridx = uniform2(generator);
  1255. int nx = x, ny = y, ni = ny * width + nx;
  1256. nx = dataset[t].ptr()[indices[0][ridx] * 2];
  1257. ny = dataset[t].ptr()[indices[0][ridx] * 2 + 1];
  1258. ni = ny * width + nx;
  1259. hypothesis* top_n_h = NULL;
  1260. // select best hypothesis at nx,ny!
  1261. // hypotheses are sort in ascending order (the last selected hypothesis is the first, the others are sorted using the data term!)
  1262. top_n_h = new hypothesis(tmp[ni][0]);
  1263. top_n_h->setLocation(ox, oy);
  1264. top_n_h->setOcclusions(forward_flow, backward_flow, occlusion_threshold, occlusion_fb_threshold);
  1265. top_n_h->energy = addJC(top_n_h, forward_flow, acc_jc, acc_cv, phi_d, thread_params, occlusions) + addBCGC(top_n_h, data, data_dx, data_dy, acc_bc, acc_gc, skip_pixel, thread_params, occlusions)
  1266. + addOC(top_n_h, acc_occ, acc_temporal_occ, thread_params)
  1267. + weight_jet_estimation[top_n_h->jet_est];
  1268. if(top_n_h != NULL) {
  1269. // check if a similar hypothesis does not already exists
  1270. bool discard = false;
  1271. for (uint32_t h = 0; h < fb_hypotheses[i].size() && !discard; h++) {
  1272. // only store dissimilar or longer trajectories
  1273. discard = (fb_hypotheses[i][h]->compare((*top_n_h), traj_sim_thres, traj_sim_method) >= 0);
  1274. }
  1275. if(!discard) {
  1276. fb_hypotheses[i].push_back(top_n_h);
  1277. } else
  1278. delete top_n_h;
  1279. }
  1280. }
  1281. if(fb_hypotheses[i].size() == 0)
  1282. cout << found_nn << endl;
  1283. }
  1284. }
  1285. }
  1286. delete dataset[0].ptr();
  1287. delete dataset[1].ptr();
  1288. delete nonempty[0];
  1289. delete nonempty[1];
  1290. clock_t end_nn = clock();
  1291. cout << "NN propagation took " << (double) (end_nn-start_nn) / CLOCKS_PER_SEC << "secs" << endl;
  1292. }
  1293. /*
  1294. * ########################################### non maximum suppression ###########################################
  1295. */
  1296. int all_h = 0;
  1297. int remaining_h = 0;
  1298. for (int y = 0; y < height; y++) {
  1299. for (int x = 0; x < width; x++) {
  1300. int idx = y * width + x;
  1301. all_h += fb_hypotheses[idx].size();
  1302. if (fb_hypotheses[idx].size() > 1) {
  1303. vector<hypothesis*> nms_hypotheses;
  1304. // resort after adding neighboring hypotheses!
  1305. // hypotheses are sort in ascending order (the last selected hypothesis is the first, the others are sorted using the data term!)
  1306. if(p_it > 0)
  1307. sort(fb_hypotheses[idx].begin() + 1, fb_hypotheses[idx].end(), compareHypotheses);
  1308. else
  1309. sort(fb_hypotheses[idx].begin(), fb_hypotheses[idx].end(), compareHypotheses);
  1310. nms_hypotheses.push_back(fb_hypotheses[idx][0]);
  1311. // iterate over all other hypotheses
  1312. for (uint32_t nh = 1; nh < fb_hypotheses[idx].size(); nh++) {
  1313. bool discard = false;
  1314. // compare to already added hypotheses
  1315. for (uint32_t h = 0; h < nms_hypotheses.size(); h++) {
  1316. if (fb_hypotheses[idx][nh]->distance(*nms_hypotheses[h], traj_sim_method) < traj_sim_thres)
  1317. discard = true;
  1318. }
  1319. // add or delete
  1320. if (!discard) {
  1321. nms_hypotheses.push_back(fb_hypotheses[idx][nh]);
  1322. } else {
  1323. delete fb_hypotheses[idx][nh];
  1324. fb_hypotheses[idx][nh] = NULL;
  1325. break;
  1326. }
  1327. }
  1328. fb_hypotheses[idx].clear();
  1329. fb_hypotheses[idx] = nms_hypotheses;
  1330. }
  1331. remaining_h += fb_hypotheses[idx].size();
  1332. }
  1333. }
  1334. if (thread_params.verbosity(VER_CMD))
  1335. cout << endl << "From " << all_h << " " << remaining_h << " hypotheses are remaining after non maxima suppression!" << endl;
  1336. /*
  1337. * ########################################### select minimal hypothesis ###########################################
  1338. */
  1339. for (int y = 0; y < height; y++) {
  1340. for (int x = 0; x < width; x++) {
  1341. int i = y * width + x;
  1342. if (fb_hypotheses[i].empty())
  1343. continue;
  1344. vector<hypothesis*> sorted_hyp = fb_hypotheses[i];
  1345. sort(sorted_hyp.begin(), sorted_hyp.end(), compareHypotheses);
  1346. }
  1347. }
  1348. // initialize MRF object
  1349. numVariables = height * width; // maximum number of variables
  1350. factors = 0;
  1351. // create MRF object
  1352. MRFEnergy<TypeGeneral>* mrf = new MRFEnergy<TypeGeneral>(TypeGeneral::GlobalSize());
  1353. vector<MRFEnergy<TypeGeneral>::NodeId> nodes(height * width);
  1354. TypeGeneral::REAL energy, lowerBound;
  1355. cout << "Generate MRF object..." << endl;
  1356. // data term: add unary nodes
  1357. time(&unary_start);
  1358. cout << " adding unary potentials ..." << flush;
  1359. // everything visible in first frame
  1360. for (int y = 0; y < height; y++) {
  1361. for (int x = 0; x < width; x++) {
  1362. int idx = y * width + x;
  1363. uint32_t labels = fb_hypotheses[idx].size();
  1364. if(labels == 0) {
  1365. throw std::out_of_range("One pixel without hypotheses! Please add outlier hypotheses!");
  1366. }
  1367. double* e_mex = new double[labels];
  1368. // add data terms
  1369. for (uint32_t h = 0; h < labels; h++) {
  1370. e_mex[h] = fb_hypotheses[idx][h]->energy;
  1371. }
  1372. if (thread_params.verbosity(VER_CMD)) {
  1373. if (x == 120 && y == 90) {
  1374. cout << "Unary: P(" << x << ", " << y << ")" << endl;
  1375. for (uint32_t h = 0; h < labels; h++)
  1376. cout << " " << h;
  1377. cout << endl;
  1378. for (uint32_t h = 0; h < labels; h++)
  1379. cout << " " << e_mex[h];
  1380. cout << endl << endl;
  1381. }
  1382. }
  1383. nodes[idx] = mrf->AddNode(TypeGeneral::LocalSize(labels), TypeGeneral::NodeData(e_mex));
  1384. delete[] e_mex;
  1385. factors++;
  1386. }
  1387. }
  1388. time(&unary_end);
  1389. cout << " took " << difftime(unary_end, unary_start) << endl;
  1390. #pragma omp critical
  1391. {
  1392. avg_unary_time += difftime(unary_end, unary_start);
  1393. }
  1394. time(&unary_start);
  1395. cout << " adding pairwise potentials ..." << flush;
  1396. // adding pairwise
  1397. for (int y = 0; y < height; y++) {
  1398. for (int x = 0; x < width; x++) {
  1399. int idx1 = y * width + x;
  1400. int oidx1 = (y * xy_incr + xy_start) * owidth + x * xy_incr + xy_start;
  1401. uint32_t labels1 = fb_hypotheses[idx1].size();
  1402. if (labels1 == 0)
  1403. continue;
  1404. vector<int> oneighbors;
  1405. vector<int> neighbors;
  1406. if (x + 1 < width) {
  1407. neighbors.push_back(y * width + x + 1);
  1408. oneighbors.push_back((y * xy_incr + xy_start) * owidth + (x + 1) * xy_incr + xy_start);
  1409. }
  1410. if (y + 1 < height) {
  1411. neighbors.push_back((y + 1) * width + x);
  1412. oneighbors.push_back(((y + 1) * xy_incr + xy_start) * owidth + x * xy_incr + xy_start);
  1413. }
  1414. vector<int>::iterator oidx2 = oneighbors.begin();
  1415. for (vector<int>::iterator idx2 = neighbors.begin(); idx2 != neighbors.end(); idx2++, oidx2++) {
  1416. uint32_t labels2 = fb_hypotheses[*idx2].size();
  1417. if (labels2 == 0)
  1418. continue;
  1419. // create energy matrix
  1420. TypeGeneral::REAL *P = new TypeGeneral::REAL[labels1 * labels2];
  1421. for (uint32_t h1 = 0; h1 < labels1; h1++) {
  1422. for (uint32_t h2 = 0; h2 < labels2; h2++) {
  1423. float dist = outlier_beta;
  1424. float smooth_occ = 0;
  1425. hypothesis* hyp_p1 = fb_hypotheses[idx1][h1];
  1426. hypothesis* hyp_p2 = fb_hypotheses[*idx2][h2];
  1427. // similarity of hypotheses
  1428. dist = hyp_p1->distance(*hyp_p2, traj_sim_method);
  1429. // spatial smoothness for occlusion
  1430. for (uint32_t j = 0; j < Jets + 1; j++)
  1431. if (hyp_p1->occluded(j) != hyp_p2->occluded(j))
  1432. smooth_occ++;
  1433. P[h2 * labels1 + h1] = (smooth_weight->data[oidx1] + smooth_weight->data[*oidx2]) * (acc_beta * dist + acc_spatial_occ * smooth_occ);
  1434. }
  1435. }
  1436. if (thread_params.verbosity(VER_CMD)) {
  1437. if (x == 120 && y == 90) {
  1438. // ##### DEBUG EXTENDED OUTPUT
  1439. if (thread_params.verbosity(VER_CMD)) {
  1440. cout << "Pairwise: SP(" << x << ", " << y << "), SP(" << *idx2 << ")" << endl;
  1441. cout << " ";
  1442. for (uint32_t h_2 = 0; h_2 < labels2; h_2++) {
  1443. cout << "h" << h_2 << " ";
  1444. }
  1445. cout << endl;
  1446. for (uint32_t h_1 = 0; h_1 < labels1; h_1++) {
  1447. cout << " h" << h_1 << " ";
  1448. for (uint32_t h_2 = 0; h_2 < labels2; h_2++) {
  1449. cout << P[h_2 * labels1 + h_1] << " ";
  1450. }
  1451. cout << endl;
  1452. }
  1453. cout << endl;
  1454. }
  1455. }
  1456. }
  1457. mrf->AddEdge(nodes[idx1], nodes[*idx2], TypeGeneral::EdgeData(TypeGeneral::GENERAL, P));
  1458. delete[] P;
  1459. factors++;
  1460. }
  1461. }
  1462. }
  1463. time(&unary_end);
  1464. cout << " took " << difftime(unary_end, unary_start) << endl;
  1465. #pragma omp critical
  1466. {
  1467. avg_pairw_time += difftime(unary_end, unary_start);
  1468. }
  1469. cout << "Variables:\t" << numVariables << endl;
  1470. cout << "Factors:\t" << factors << endl;
  1471. // run discrete optimization TRW-S
  1472. cout << "Run discrete optimization..." << endl;
  1473. time(&opt_start);
  1474. mrf->SetAutomaticOrdering();
  1475. if (params.parameter<int>("acc_approach") == 0) {
  1476. mrf->Minimize_TRW_S(options, lowerBound, energy);
  1477. time(&opt_end);
  1478. printf("TRW-S finished. Time: %i\n", (int) difftime(opt_end, opt_start));
  1479. } else {
  1480. mrf->Minimize_BP(options, energy);
  1481. lowerBound = std::numeric_limits<double>::signaling_NaN();
  1482. time(&opt_end);
  1483. printf("BP finished. Time: %i\n", (int) difftime(opt_end, opt_start));
  1484. }
  1485. #pragma omp critical
  1486. {
  1487. avg_optimization_time += difftime(opt_end, opt_start);
  1488. }
  1489. // get solutions
  1490. Mat oracle_selection = Mat::zeros(height, width, CV_8UC3);
  1491. Mat oracle_present = Mat::zeros(height, width, CV_8UC3);
  1492. Mat occlusion_map = Mat::zeros(height, width, CV_32FC1);
  1493. Mat flow = Mat::zeros(height, width, CV_64FC2);
  1494. for (int y = 0; y < height; y++) {
  1495. for (int x = 0; x < width; x++) {
  1496. int idx = y * width + x;
  1497. int h = -1;
  1498. if (fb_hypotheses[idx].size() > 0)
  1499. h = (int) (mrf->GetSolution(nodes[idx]));
  1500. // remember best hypothesis
  1501. selected_hyp.at<int>(y, x) = h;
  1502. // get flow of best hypothesis
  1503. flow.at<Vec2d>(y, x)[1] = fb_hypotheses[idx][h]->u(Jets - 1) / xy_incr;
  1504. flow.at<Vec2d>(y, x)[0] = fb_hypotheses[idx][h]->v(Jets - 1) / xy_incr;
  1505. // get occlusion map
  1506. occlusion_map.at<float>(y, x) = fb_hypotheses[idx][h]->occluded(0);
  1507. for (uint32_t f = 0; f < Jets; f++)
  1508. occlusion_map.at<float>(y, x) = max((int) occlusion_map.at<float>(y, x), fb_hypotheses[idx][h]->occluded(f + 1));
  1509. }
  1510. }
  1511. #pragma omp critical
  1512. {
  1513. numVariablesStream << "\t" << numVariables;
  1514. factorsStream << "\t" << factors;
  1515. }
  1516. delete mrf;
  1517. /*
  1518. * ########################################### visualization and store accumulated flow ###########################################
  1519. */
  1520. if (p_it == (alternate - 1)) {
  1521. Mat flow_img = flowColorImg(flow, thread_params.verbosity(VER_CMD));
  1522. // Check for invalid input
  1523. if (!flow_img.data) {
  1524. cout << "No flow for frame " << thread_params.sequence_start << std::endl;
  1525. continue;
  1526. }
  1527. // write final estimation
  1528. char occF[1024];
  1529. sprintf(occF, (acc_folder.str() + "/occlusions/frame_%i.pbm").c_str(), thread_params.sequence_start);
  1530. occlusion_map.convertTo(occlusion_map, CV_8UC1, 255);
  1531. imwrite(occF, occlusion_map);
  1532. char forward_flow_file[1024];
  1533. if (!sintel)
  1534. sprintf(forward_flow_file, (acc_folder.str() + "/" + flow_format).c_str(), thread_params.sequence_start);
  1535. else
  1536. sprintf(forward_flow_file, (acc_folder.str() + "/" + flow_format).c_str(), thread_params.sequence_start, 0);
  1537. imwrite((string(forward_flow_file) + "_vis.png"), flow_img);
  1538. writeFlowMiddlebury(flow, (string(forward_flow_file) + ".flo"));
  1539. }
  1540. }
  1541. cout << endl << "Writing results to: " << acc_folder.str() << endl << endl;
  1542. //----------------------------------
  1543. // ########################################### final clean up ###########################################
  1544. //----------------------------------
  1545. for (uint32_t i = 0; i < fb_hypotheses.size(); i++) {
  1546. for (uint32_t h = 0; h < fb_hypotheses[i].size(); h++) {
  1547. delete fb_hypotheses[i][h];
  1548. }
  1549. }
  1550. fb_hypotheses.clear();
  1551. fb_hypotheses = vector<vector<hypothesis*> >();
  1552. for (uint32_t f = 0; f < (Jets + 1); f++) {
  1553. color_image_delete(data[f]);
  1554. color_image_delete(data_dx[f]);
  1555. color_image_delete(data_dy[f]);
  1556. }
  1557. delete[] data;
  1558. delete[] data_dx;
  1559. delete[] data_dy;
  1560. image_delete(smooth_weight);
  1561. if(imlab != NULL) {
  1562. color_image_delete(imlab);
  1563. free(forward_edges.pixels);
  1564. }
  1565. if(epic_wx != NULL) image_delete(epic_wx);
  1566. if(epic_wy != NULL) image_delete(epic_wy);
  1567. delete[] sequence;
  1568. if (gt_sfr != NULL) {
  1569. image_delete(gt_sfr[0][0]);
  1570. image_delete(gt_sfr[0][1]);
  1571. delete[] gt_sfr[0];
  1572. delete[] gt_sfr;
  1573. }
  1574. if (gt != NULL) delete[] gt;
  1575. if (gt_occlusions != NULL) delete[] gt_occlusions;
  1576. if (occlusions != NULL) delete[] occlusions;
  1577. delete[] forward_flow;
  1578. delete[] backward_flow;
  1579. /*
  1580. * write results and infos to file
  1581. */
  1582. if (!params.output.empty()) {
  1583. ofstream infos;
  1584. infos.open((acc_folder.str() + "/result.info").c_str());
  1585. infos << "# Discrete optimization file\n\n";
  1586. infos << "Warping took " << dt_warp_time << "s.\n";
  1587. infos << "Median took " << dt_med_time << "s.\n";
  1588. infos << "Data term computation took " << dt_sum_time << "s.\n";
  1589. infos << "Adding unary potentials took " << avg_unary_time << "s.\n";
  1590. infos << "Adding pairwise potentials took " << avg_pairw_time << "s.\n";
  1591. infos << "Run discrete optimization took " << avg_optimization_time << "s.\n\n";
  1592. infos << "Discrete Optimization:\n";
  1593. infos << "\tVariables:\t" << numVariablesStream.str() << "\n";
  1594. infos << "\tFactors:\t" << factorsStream.str() << "\n\n";
  1595. infos.close();
  1596. }
  1597. }
  1598. if (phi_d != NULL)
  1599. delete phi_d;
  1600. if (phi_d != NULL)
  1601. delete phi_s;
  1602. cout << "Done!" << endl;
  1603. return 0;
  1604. }