SemSegContextTree3D.cpp 65 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276
  1. #include "SemSegContextTree3D.h"
  2. #include "core/basics/FileName.h"
  3. #include "core/basics/numerictools.h"
  4. #include "core/basics/quadruplet.h"
  5. #include "core/basics/StringTools.h"
  6. #include "core/basics/Timer.h"
  7. #include "core/basics/vectorio.h"
  8. #include "core/image/Filter.h"
  9. #include "core/image/FilterT.h"
  10. #include "core/image/Morph.h"
  11. #include "core/imagedisplay/ImageDisplay.h"
  12. #include "vislearning/baselib/cc.h"
  13. #include "vislearning/baselib/Globals.h"
  14. #include "vislearning/baselib/ICETools.h"
  15. #include "vislearning/cbaselib/CachedExample.h"
  16. #include "vislearning/cbaselib/PascalResults.h"
  17. #include "segmentation/RSGraphBased.h"
  18. #include "segmentation/RSMeanShift.h"
  19. #include "segmentation/RSSlic.h"
  20. #include <omp.h>
  21. #include <iostream>
  22. #define VERBOSE
  23. #undef DEBUG
  24. #undef VISUALIZE
  25. #undef WRITEREGIONS
  26. using namespace OBJREC;
  27. using namespace std;
  28. using namespace NICE;
  29. //###################### CONSTRUCTORS #########################//
  30. SemSegContextTree3D::SemSegContextTree3D () : SemanticSegmentation ()
  31. {
  32. this->lfcw = NULL;
  33. this->firstiteration = true;
  34. this->run3Dseg = false;
  35. this->maxSamples = 2000;
  36. this->minFeats = 50;
  37. this->maxDepth = 10;
  38. this->windowSize = 15;
  39. this->contextMultiplier = 3;
  40. this->featsPerSplit = 200;
  41. this->useShannonEntropy = true;
  42. this->nbTrees = 10;
  43. this->randomTests = 10;
  44. this->useAltTristimulus = false;
  45. this->useGradient = true;
  46. this->useWeijer = false;
  47. this->useAdditionalLayer = false;
  48. this->useHoiemFeatures = false;
  49. this->useCategorization = false;
  50. this->cndir = "";
  51. this->fasthik = NULL;
  52. this->saveLoadData = false;
  53. this->fileLocation = "tmp.txt";
  54. this->pixelWiseLabeling = true;
  55. this->segmentation = NULL;
  56. this->useFeat0 = true;
  57. this->useFeat1 = false;
  58. this->useFeat2 = true;
  59. this->useFeat3 = true;
  60. this->useFeat4 = false;
  61. this->useFeat5 = false;
  62. }
  63. SemSegContextTree3D::SemSegContextTree3D (
  64. const Config *conf,
  65. const MultiDataset *md )
  66. : SemanticSegmentation ( conf, & ( md->getClassNames ( "train" ) ) )
  67. {
  68. this->conf = conf;
  69. string section = "SSContextTree";
  70. string featsec = "Features";
  71. this->lfcw = NULL;
  72. this->firstiteration = true;
  73. this->run3Dseg = conf->gB ( section, "run_3dseg", false );
  74. this->maxSamples = conf->gI ( section, "max_samples", 2000 );
  75. this->minFeats = conf->gI ( section, "min_feats", 50 );
  76. this->maxDepth = conf->gI ( section, "max_depth", 10 );
  77. this->windowSize = conf->gI ( section, "window_size", 15 );
  78. this->contextMultiplier = conf->gI ( section, "context_multiplier", 3 );
  79. this->featsPerSplit = conf->gI ( section, "feats_per_split", 200 );
  80. this->useShannonEntropy = conf->gB ( section, "use_shannon_entropy", true );
  81. this->nbTrees = conf->gI ( section, "amount_trees", 10 );
  82. this->randomTests = conf->gI ( section, "random_tests", 10 );
  83. this->useAltTristimulus = conf->gB ( featsec, "use_alt_trist", false );
  84. this->useGradient = conf->gB ( featsec, "use_gradient", true );
  85. this->useWeijer = conf->gB ( featsec, "use_weijer", true );
  86. this->useAdditionalLayer = conf->gB ( featsec, "use_additional_layer", false );
  87. this->useHoiemFeatures = conf->gB ( featsec, "use_hoiem_features", false );
  88. this->useCategorization = conf->gB ( section, "use_categorization", false );
  89. this->cndir = conf->gS ( "SSContextTree", "cndir", "" );
  90. this->saveLoadData = conf->gB ( "debug", "save_load_data", false );
  91. this->fileLocation = conf->gS ( "debug", "datafile", "tmp.txt" );
  92. this->pixelWiseLabeling = conf->gB ( section, "pixelWiseLabeling", false );
  93. if ( useCategorization && cndir == "" )
  94. this->fasthik = new GPHIKClassifierNICE ( conf );
  95. else
  96. this->fasthik = NULL;
  97. if ( useWeijer )
  98. this->lfcw = new LocalFeatureColorWeijer ( conf );
  99. this->classnames = md->getClassNames ( "train" );
  100. // feature types
  101. this->useFeat0 = conf->gB ( section, "use_feat_0", true); // pixel pair features
  102. this->useFeat1 = conf->gB ( section, "use_feat_1", false); // region feature
  103. this->useFeat2 = conf->gB ( section, "use_feat_2", true); // integral features
  104. this->useFeat3 = conf->gB ( section, "use_feat_3", true); // integral contex features
  105. this->useFeat4 = conf->gB ( section, "use_feat_4", false); // pixel pair context features
  106. this->useFeat5 = conf->gB ( section, "use_feat_5", false); // ray features
  107. string segmentationtype = conf->gS ( section, "segmentation_type", "slic" );
  108. if ( segmentationtype == "meanshift" )
  109. this->segmentation = new RSMeanShift ( conf );
  110. else if ( segmentationtype == "felzenszwalb" )
  111. this->segmentation = new RSGraphBased ( conf );
  112. else if ( segmentationtype == "slic" )
  113. this->segmentation = new RSSlic ( conf );
  114. else if ( segmentationtype == "none" )
  115. {
  116. this->segmentation = NULL;
  117. this->pixelWiseLabeling = true;
  118. this->useFeat1 = false;
  119. }
  120. else
  121. throw ( "no valid segmenation_type\n please choose between none, meanshift, slic and felzenszwalb\n" );
  122. if ( !useGradient || imagetype == IMAGETYPE_RGB)
  123. this->useFeat5 = false;
  124. if ( useFeat0 )
  125. this->featTypes.push_back(0);
  126. if ( useFeat1 )
  127. this->featTypes.push_back(1);
  128. if ( useFeat2 )
  129. this->featTypes.push_back(2);
  130. if ( useFeat3 )
  131. this->featTypes.push_back(3);
  132. if ( useFeat4 )
  133. this->featTypes.push_back(4);
  134. if ( useFeat5 )
  135. this->featTypes.push_back(5);
  136. this->initOperations();
  137. }
  138. //###################### DESTRUCTORS ##########################//
  139. SemSegContextTree3D::~SemSegContextTree3D()
  140. {
  141. }
  142. //#################### MEMBER FUNCTIONS #######################//
  143. void SemSegContextTree3D::initOperations()
  144. {
  145. string featsec = "Features";
  146. // operation prototypes
  147. vector<Operation3D*> tops0, tops1, tops2, tops3, tops4, tops5;
  148. if ( conf->gB ( featsec, "int", true ) )
  149. {
  150. tops2.push_back ( new IntegralOps3D() );
  151. Operation3D* o = new IntegralOps3D();
  152. o->setContext(true);
  153. tops3.push_back ( o );
  154. }
  155. if ( conf->gB ( featsec, "bi_int_cent", true ) )
  156. {
  157. tops2.push_back ( new BiIntegralCenteredOps3D() );
  158. Operation3D* o = new BiIntegralCenteredOps3D();
  159. o->setContext(true);
  160. tops3.push_back ( o );
  161. }
  162. if ( conf->gB ( featsec, "int_cent", true ) )
  163. {
  164. tops2.push_back ( new IntegralCenteredOps3D() );
  165. Operation3D* o = new IntegralCenteredOps3D();
  166. o->setContext(true);
  167. tops3.push_back ( o );
  168. }
  169. if ( conf->gB ( featsec, "haar_horz", true ) )
  170. {
  171. tops2.push_back ( new HaarHorizontal3D() );
  172. Operation3D* o = new HaarHorizontal3D();
  173. o->setContext(true);
  174. tops3.push_back ( o );
  175. }
  176. if ( conf->gB ( featsec, "haar_vert", true ) )
  177. {
  178. tops2.push_back ( new HaarVertical3D() );
  179. Operation3D* o = new HaarVertical3D();
  180. o->setContext(true);
  181. tops3.push_back ( o );
  182. }
  183. if ( conf->gB ( featsec, "haar_stack", true ) )
  184. {
  185. tops2.push_back ( new HaarStacked3D() );
  186. Operation3D* o = new HaarStacked3D();
  187. o->setContext(true);
  188. tops3.push_back ( o );
  189. }
  190. if ( conf->gB ( featsec, "haar_diagxy", true ) )
  191. {
  192. tops2.push_back ( new HaarDiagXY3D() );
  193. Operation3D* o = new HaarDiagXY3D();
  194. o->setContext(true);
  195. tops3.push_back ( o );
  196. }
  197. if ( conf->gB ( featsec, "haar_diagxz", true ) )
  198. {
  199. tops2.push_back ( new HaarDiagXZ3D() );
  200. Operation3D* o = new HaarDiagXZ3D();
  201. o->setContext(true);
  202. tops3.push_back ( o );
  203. }
  204. if ( conf->gB ( featsec, "haar_diagyz", true ) )
  205. {
  206. tops2.push_back ( new HaarDiagYZ3D() );
  207. Operation3D* o = new HaarDiagYZ3D();
  208. o->setContext(true);
  209. tops3.push_back ( o );
  210. }
  211. if ( conf->gB ( featsec, "haar3_horz", true ) )
  212. {
  213. tops2.push_back ( new Haar3Horiz3D() );
  214. Operation3D* o = new Haar3Horiz3D();
  215. o->setContext(true);
  216. tops3.push_back ( o );
  217. }
  218. if ( conf->gB ( featsec, "haar3_vert", true ) )
  219. {
  220. tops2.push_back ( new Haar3Vert3D() );
  221. Operation3D* o = new Haar3Vert3D();
  222. o->setContext(true);
  223. tops3.push_back ( o );
  224. }
  225. if ( conf->gB ( featsec, "haar3_stack", true ) )
  226. {
  227. tops2.push_back ( new Haar3Stack3D() );
  228. Operation3D* o = new Haar3Stack3D();
  229. o->setContext(true);
  230. tops3.push_back ( o );
  231. }
  232. if ( conf->gB ( featsec, "minus", true ) )
  233. {
  234. tops0.push_back ( new Minus3D() );
  235. Operation3D* o = new Minus3D();
  236. o->setContext(true);
  237. tops4.push_back ( o );
  238. }
  239. if ( conf->gB ( featsec, "minus_abs", true ) )
  240. {
  241. tops0.push_back ( new MinusAbs3D() );
  242. Operation3D* o = new MinusAbs3D();
  243. o->setContext(true);
  244. tops4.push_back ( o );
  245. }
  246. if ( conf->gB ( featsec, "addition", true ) )
  247. {
  248. tops0.push_back ( new Addition3D() );
  249. Operation3D* o = new Addition3D();
  250. o->setContext(true);
  251. tops4.push_back ( o );
  252. }
  253. if ( conf->gB ( featsec, "only1", true ) )
  254. {
  255. tops0.push_back ( new Only13D() );
  256. Operation3D* o = new Only13D();
  257. o->setContext(true);
  258. tops4.push_back ( o );
  259. }
  260. if ( conf->gB ( featsec, "rel_x", true ) )
  261. tops0.push_back ( new RelativeXPosition3D() );
  262. if ( conf->gB ( featsec, "rel_y", true ) )
  263. tops0.push_back ( new RelativeYPosition3D() );
  264. if ( conf->gB ( featsec, "rel_z", true ) )
  265. tops0.push_back ( new RelativeZPosition3D() );
  266. if ( conf->gB ( featsec, "ray_diff", false ) )
  267. tops5.push_back ( new RayDiff3D() );
  268. if ( conf->gB ( featsec, "ray_dist", false ) )
  269. tops5.push_back ( new RayDist3D() );
  270. if ( conf->gB ( featsec, "ray_orient", false ) )
  271. tops5.push_back ( new RayOrient3D() );
  272. if ( conf->gB ( featsec, "ray_norm", false ) )
  273. tops5.push_back ( new RayNorm3D() );
  274. this->ops.push_back ( tops0 );
  275. this->ops.push_back ( tops1 );
  276. this->ops.push_back ( tops2 );
  277. this->ops.push_back ( tops3 );
  278. this->ops.push_back ( tops4 );
  279. this->ops.push_back ( tops5 );
  280. }
  281. double SemSegContextTree3D::getBestSplit (
  282. std::vector<NICE::MultiChannelImage3DT<double> > &feats,
  283. std::vector<NICE::MultiChannelImage3DT<unsigned short int> > &nodeIndices,
  284. const std::vector<NICE::MultiChannelImageT<int> > &labels,
  285. int node,
  286. Operation3D *&splitop,
  287. double &splitval,
  288. const int &tree,
  289. vector<vector<vector<double> > > &regionProbs )
  290. {
  291. Timer t;
  292. t.start();
  293. int imgCount = 0;
  294. try
  295. {
  296. imgCount = ( int ) feats.size();
  297. }
  298. catch ( Exception )
  299. {
  300. cerr << "no features computed?" << endl;
  301. }
  302. double bestig = -numeric_limits< double >::max();
  303. splitop = NULL;
  304. splitval = -1.0;
  305. vector<quadruplet<int,int,int,int> > selFeats;
  306. map<int, int> e;
  307. int featcounter = forest[tree][node].featcounter;
  308. if ( featcounter < minFeats )
  309. {
  310. return 0.0;
  311. }
  312. vector<double> fraction ( a.size(), 0.0 );
  313. for ( uint i = 0; i < fraction.size(); i++ )
  314. {
  315. if ( forbidden_classes.find ( labelmapback[i] ) != forbidden_classes.end() )
  316. fraction[i] = 0;
  317. else
  318. fraction[i] = ( ( double ) maxSamples ) / ( ( double ) featcounter * a[i] * a.size() );
  319. }
  320. featcounter = 0;
  321. for ( int iCounter = 0; iCounter < imgCount; iCounter++ )
  322. {
  323. int xsize = ( int ) nodeIndices[iCounter].width();
  324. int ysize = ( int ) nodeIndices[iCounter].height();
  325. int zsize = ( int ) nodeIndices[iCounter].depth();
  326. for ( int x = 0; x < xsize; x++ )
  327. for ( int y = 0; y < ysize; y++ )
  328. for ( int z = 0; z < zsize; z++ )
  329. {
  330. if ( nodeIndices[iCounter].get ( x, y, z, tree ) == node )
  331. {
  332. int cn = labels[iCounter].get ( x, y, ( uint ) z );
  333. double randD = ( double ) rand() / ( double ) RAND_MAX;
  334. if ( labelmap.find ( cn ) == labelmap.end() )
  335. continue;
  336. if ( randD < fraction[labelmap[cn]] )
  337. {
  338. quadruplet<int,int,int,int> quad( iCounter, x, y, z );
  339. featcounter++;
  340. selFeats.push_back ( quad );
  341. e[cn]++;
  342. }
  343. }
  344. }
  345. }
  346. // global entropy
  347. double globent = 0.0;
  348. for ( map<int, int>::iterator mapit = e.begin() ; mapit != e.end(); mapit++ )
  349. {
  350. double p = ( double ) ( *mapit ).second / ( double ) featcounter;
  351. globent += p * log2 ( p );
  352. }
  353. globent = -globent;
  354. if ( globent < 0.5 )
  355. return 0.0;
  356. // pointers to all randomly chosen features
  357. std::vector<Operation3D*> featsel;
  358. for ( int i = 0; i < featsPerSplit; i++ )
  359. {
  360. int x1, x2, y1, y2, z1, z2, ft;
  361. do
  362. {
  363. ft = ( int ) ( rand() % featTypes.size() );
  364. ft = featTypes[ft];
  365. }
  366. while ( channelsPerType[ft].size() == 0 );
  367. int tmpws = windowSize;
  368. if ( ft == 2 || ft == 4 )
  369. {
  370. //use larger window size for context features
  371. tmpws *= contextMultiplier;
  372. }
  373. // use region feature only with reasonable pre-segmentation
  374. // if ( ft == 1 && depth < 8 )
  375. // {
  376. // ft = 0;
  377. // }
  378. /* random window positions */
  379. double z_ratio = conf->gB ( "SSContextTree", "z_ratio", 1.0 );
  380. int tmp_z = ( int ) floor( (tmpws * z_ratio) + 0.5 );
  381. double y_ratio = conf->gB ( "SSContextTree", "y_ratio", 1.0 );
  382. int tmp_y = ( int ) floor( (tmpws * y_ratio) + 0.5 );
  383. x1 = ( int ) ( rand() % tmpws ) - tmpws / 2 ;
  384. x2 = ( int ) ( rand() % tmpws ) - tmpws / 2 ;
  385. y1 = ( int ) ( rand() % tmp_y ) - tmp_y / 2 ;
  386. y2 = ( int ) ( rand() % tmp_y ) - tmp_y / 2 ;
  387. z1 = ( int ) ( rand() % tmp_z ) - tmp_z / 2 ;
  388. z2 = ( int ) ( rand() % tmp_z ) - tmp_z / 2 ;
  389. // use z1/z2 as directions (angles) in ray features
  390. if ( ft == 5 )
  391. {
  392. z1 = ( int ) ( rand() % 8 );
  393. z2 = ( int ) ( rand() % 8 );
  394. }
  395. // if (conf->gB ( "SSContextTree", "z_negative_only", false ))
  396. // {
  397. // z1 = -abs(z1);
  398. // z2 = -abs(z2);
  399. // }
  400. /* random feature maps (channels) */
  401. int f1, f2;
  402. f1 = ( int ) ( rand() % channelsPerType[ft].size() );
  403. if ( (rand() % 2) == 0 )
  404. f2 = ( int ) ( rand() % channelsPerType[ft].size() );
  405. else
  406. f2 = f1;
  407. f1 = channelsPerType[ft][f1];
  408. f2 = channelsPerType[ft][f2];
  409. if ( ft == 1 )
  410. {
  411. int classes = ( int ) regionProbs[0][0].size();
  412. f2 = ( int ) ( rand() % classes );
  413. }
  414. /* random extraction method (operation) */
  415. int o = ( int ) ( rand() % ops[ft].size() );
  416. Operation3D *op = ops[ft][o]->clone();
  417. op->set ( x1, y1, z1, x2, y2, z2, f1, f2, ft );
  418. if ( ft == 3 || ft == 4 )
  419. op->setContext ( true );
  420. else
  421. op->setContext ( false );
  422. featsel.push_back ( op );
  423. }
  424. // do actual split tests
  425. for ( int f = 0; f < featsPerSplit; f++ )
  426. {
  427. double l_bestig = -numeric_limits< double >::max();
  428. double l_splitval = -1.0;
  429. vector<double> vals;
  430. double maxval = -numeric_limits<double>::max();
  431. double minval = numeric_limits<double>::max();
  432. int counter = 0;
  433. for ( vector<quadruplet<int,int,int,int> >::const_iterator it = selFeats.begin();
  434. it != selFeats.end(); it++ )
  435. {
  436. Features feat;
  437. feat.feats = &feats[ ( *it ).first ];
  438. feat.rProbs = &regionProbs[ ( *it ).first ];
  439. assert ( forest.size() > ( uint ) tree );
  440. assert ( forest[tree][0].dist.size() > 0 );
  441. double val = 0.0;
  442. val = featsel[f]->getVal ( feat, ( *it ).second, ( *it ).third, ( *it ).fourth );
  443. if ( !isfinite ( val ) )
  444. {
  445. #ifdef DEBUG
  446. cerr << "feat " << feat.feats->width() << " " << feat.feats->height() << " " << feat.feats->depth() << endl;
  447. cerr << "non finite value " << val << " for " << featsel[f]->writeInfos() << endl << (*it).second << " " << (*it).third << " " << (*it).fourth << endl;
  448. #endif
  449. val = 0.0;
  450. }
  451. vals.push_back ( val );
  452. maxval = std::max ( val, maxval );
  453. minval = std::min ( val, minval );
  454. }
  455. if ( minval == maxval )
  456. continue;
  457. // split values
  458. for ( int run = 0 ; run < randomTests; run++ )
  459. {
  460. // choose threshold randomly
  461. double sval = 0.0;
  462. sval = ( (double) rand() / (double) RAND_MAX*(maxval-minval) ) + minval;
  463. map<int, int> eL, eR;
  464. int counterL = 0, counterR = 0;
  465. counter = 0;
  466. for ( vector<quadruplet<int,int,int,int> >::const_iterator it2 = selFeats.begin();
  467. it2 != selFeats.end(); it2++, counter++ )
  468. {
  469. int cn = labels[ ( *it2 ).first ].get ( ( *it2 ).second, ( *it2 ).third, ( *it2 ).fourth );
  470. //cout << "vals[counter2] " << vals[counter2] << " val: " << val << endl;
  471. if ( vals[counter] < sval )
  472. {
  473. //left entropie:
  474. eL[cn] = eL[cn] + 1;
  475. counterL++;
  476. }
  477. else
  478. {
  479. //right entropie:
  480. eR[cn] = eR[cn] + 1;
  481. counterR++;
  482. }
  483. }
  484. double leftent = 0.0;
  485. for ( map<int, int>::iterator mapit = eL.begin() ; mapit != eL.end(); mapit++ )
  486. {
  487. double p = ( double ) ( *mapit ).second / ( double ) counterL;
  488. leftent -= p * log2 ( p );
  489. }
  490. double rightent = 0.0;
  491. for ( map<int, int>::iterator mapit = eR.begin() ; mapit != eR.end(); mapit++ )
  492. {
  493. double p = ( double ) ( *mapit ).second / ( double ) counterR;
  494. rightent -= p * log2 ( p );
  495. }
  496. //cout << "rightent: " << rightent << " leftent: " << leftent << endl;
  497. double pl = ( double ) counterL / ( double ) ( counterL + counterR );
  498. //information gain
  499. double ig = globent - ( 1.0 - pl ) * rightent - pl * leftent;
  500. //double ig = globent - rightent - leftent;
  501. if ( useShannonEntropy )
  502. {
  503. double esplit = - ( pl * log ( pl ) + ( 1 - pl ) * log ( 1 - pl ) );
  504. ig = 2 * ig / ( globent + esplit );
  505. }
  506. if ( ig > l_bestig )
  507. {
  508. l_bestig = ig;
  509. l_splitval = sval;
  510. }
  511. }
  512. if ( l_bestig > bestig )
  513. {
  514. bestig = l_bestig;
  515. splitop = featsel[f];
  516. splitval = l_splitval;
  517. }
  518. }
  519. #ifdef DEBUG
  520. cout << "globent: " << globent << " bestig " << bestig << " splitval: " << splitval << endl;
  521. #endif
  522. return bestig;
  523. }
  524. inline double SemSegContextTree3D::getMeanProb (
  525. const int &x,
  526. const int &y,
  527. const int &z,
  528. const int &channel,
  529. const MultiChannelImage3DT<unsigned short int> &nodeIndices )
  530. {
  531. double val = 0.0;
  532. for ( int tree = 0; tree < nbTrees; tree++ )
  533. {
  534. val += forest[tree][nodeIndices.get ( x,y,z,tree ) ].dist[channel];
  535. }
  536. return val / ( double ) nbTrees;
  537. }
  538. void SemSegContextTree3D::computeRayFeatImage (
  539. NICE::MultiChannelImage3DT<double> &feats,
  540. int firstChannel )
  541. {
  542. int xsize = feats.width();
  543. int ysize = feats.height();
  544. int zsize = feats.depth();
  545. const int amountDirs = 8;
  546. // compute ray feature maps from canny image
  547. for ( int z = 0; z < zsize; z++)
  548. {
  549. // canny image from raw channel
  550. NICE::Image med (xsize,ysize);
  551. NICE::median ( feats.getChannel( z, 0 ), &med, 2);
  552. NICE::Image* can = NICE::canny( med, 5, 25);
  553. for ( int dir = 0; dir < amountDirs; dir++)
  554. {
  555. NICE::Matrix dist(xsize,ysize,0);
  556. NICE::Matrix norm(xsize,ysize,0);
  557. NICE::Matrix orient(xsize,ysize,0);
  558. for (int y = 0; y < ysize; y++)
  559. for ( int x = 0; x < xsize; x++)
  560. {
  561. int xo = 0, yo = 0; // offsets
  562. int theta = 0;
  563. switch (dir)
  564. {
  565. case 0: theta = 0; yo = -1; break;
  566. case 1: theta = 45; xo = 1; yo = -1; break;
  567. case 2: theta = 90; xo = 1; x = (xsize-1)-x; break;
  568. case 3: theta = 135; xo = 1; yo = -1; break;
  569. case 4: theta = 180; yo = 1; y = (ysize-1)-y; break;
  570. case 5: theta = 225; xo = -1; yo = 1; y = (ysize-1)-y; break;
  571. case 6: theta = 270; xo = -1; break;
  572. case 7: theta = 315; xo = -1; yo = -1; break;
  573. default: return;
  574. }
  575. if (can->getPixelQuick(x,y) != 0
  576. || x+xo < 0
  577. || x+xo >= xsize
  578. || y+yo < 0
  579. || y+yo >= ysize )
  580. {
  581. double gx = feats.get(x, y, z, 1);
  582. double gy = feats.get(x, y, z, 2);
  583. //double go = atan2 (gy, gx);
  584. norm(x, y) = sqrt(gx*gx+gy*gy);
  585. orient(x, y) = ( gx*cos(theta)+gy*sin(theta) ) / norm(x,y);
  586. dist(x, y) = 0;
  587. }
  588. else
  589. {
  590. orient(x, y) = orient(x+xo,y+yo);
  591. norm(x, y) = norm(x+xo,y+yo);
  592. dist(x, y) = dist(x+xo,y+yo) + 1;
  593. }
  594. }
  595. for (int y = 0; y < ysize; y++)
  596. for (int x = 0; x < xsize; x++)
  597. {
  598. // distance feature maps
  599. feats.set( x, y, z, dist(x,y), firstChannel + dir );
  600. // norm feature maps
  601. feats.set( x, y, z, norm(x,y), firstChannel + amountDirs + dir );
  602. // orientation feature maps
  603. feats.set( x, y, z, norm(x,y), firstChannel + (amountDirs*2) + dir );
  604. }
  605. }
  606. delete can;
  607. }
  608. }
  609. void SemSegContextTree3D::updateProbabilityMaps (
  610. const NICE::MultiChannelImage3DT<unsigned short int> &nodeIndices,
  611. NICE::MultiChannelImage3DT<double> &feats,
  612. int firstChannel )
  613. {
  614. int xsize = feats.width();
  615. int ysize = feats.height();
  616. int zsize = feats.depth();
  617. int classes = ( int ) forest[0][0].dist.size();
  618. // integral images for context channels (probability maps for each class)
  619. #pragma omp parallel for
  620. for ( int c = 0; c < classes; c++ )
  621. {
  622. for ( int z = 0; z < zsize; z++ )
  623. {
  624. for ( int y = 0; y < ysize; y++ )
  625. {
  626. for ( int x = 0; x < xsize; x++ )
  627. {
  628. double val = getMeanProb ( x, y, z, c, nodeIndices );
  629. if (useFeat3)
  630. feats ( x, y, z, firstChannel + c ) = val;
  631. if (useFeat4)
  632. feats ( x, y, z, firstChannel + classes + c ) = val;
  633. }
  634. }
  635. // Gaussian filter on probability maps
  636. // NICE::ImageT<double> img = feats.getChannelT( z, firstChannel+c );
  637. // NICE::ImageT<double> gF(xsize,ysize);
  638. // NICE::FilterT<double,double,double> filt;
  639. // filt.filterGaussSigmaApproximate( img, 2, &gF );
  640. // for ( int y = 0; y < ysize; y++ )
  641. // for ( int x = 0; x < xsize; x++ )
  642. // feats.set(x, y, z, gF.getPixelQuick(x,y), firstChannel+c);
  643. }
  644. feats.calcIntegral ( firstChannel + c );
  645. }
  646. }
  647. inline double computeWeight ( const int &d, const int &dim )
  648. {
  649. if (d == 0)
  650. return 0.0;
  651. else
  652. return 1.0 / ( pow ( 2, ( double ) ( dim - d + 1 ) ) );
  653. }
  654. void SemSegContextTree3D::train ( const MultiDataset *md )
  655. {
  656. const LabeledSet trainSet = * ( *md ) ["train"];
  657. const LabeledSet *trainp = &trainSet;
  658. if ( saveLoadData )
  659. {
  660. if ( FileMgt::fileExists ( fileLocation ) )
  661. read ( fileLocation );
  662. else
  663. {
  664. train ( trainp );
  665. write ( fileLocation );
  666. }
  667. }
  668. else
  669. {
  670. train ( trainp );
  671. }
  672. }
  673. void SemSegContextTree3D::train ( const LabeledSet * trainp )
  674. {
  675. int shortsize = numeric_limits<short>::max();
  676. Timer timer;
  677. timer.start();
  678. vector<int> zsizeVec;
  679. getDepthVector ( trainp, zsizeVec, run3Dseg );
  680. //FIXME: memory usage
  681. vector<MultiChannelImage3DT<double> > allfeats;
  682. vector<MultiChannelImage3DT<unsigned short int> > nodeIndices;
  683. vector<MultiChannelImageT<int> > labels;
  684. vector<SparseVector*> globalCategorFeats;
  685. vector<map<int,int> > classesPerImage;
  686. vector<vector<int> > rSize;
  687. vector<int> amountRegionpI;
  688. std::string forbidden_classes_s = conf->gS ( "analysis", "forbidden_classes", "" );
  689. classnames.getSelection ( forbidden_classes_s, forbidden_classes );
  690. int imgCounter = 0;
  691. int amountPixels = 0;
  692. // How many channels of non-integral type do we have?
  693. if ( imagetype == IMAGETYPE_RGB )
  694. rawChannels = 3;
  695. else
  696. rawChannels = 1;
  697. if ( useGradient )
  698. rawChannels *= 3;
  699. if ( useWeijer )
  700. rawChannels += 11;
  701. if ( useHoiemFeatures )
  702. rawChannels += 8;
  703. if ( useAdditionalLayer )
  704. rawChannels += 1;
  705. ///////////////////////////// read input data /////////////////////////////////
  706. ///////////////////////////////////////////////////////////////////////////////
  707. int depthCount = 0;
  708. vector< string > filelist;
  709. NICE::MultiChannelImageT<uchar> pixelLabels;
  710. for (LabeledSet::const_iterator it = trainp->begin(); it != trainp->end(); it++)
  711. {
  712. for (std::vector<ImageInfo *>::const_iterator jt = it->second.begin();
  713. jt != it->second.end(); jt++)
  714. {
  715. int classno = it->first;
  716. ImageInfo & info = *(*jt);
  717. std::string file = info.img();
  718. filelist.push_back ( file );
  719. depthCount++;
  720. const LocalizationResult *locResult = info.localization();
  721. // getting groundtruth
  722. NICE::Image pL;
  723. pL.resize ( locResult->xsize, locResult->ysize );
  724. pL.set ( 0 );
  725. locResult->calcLabeledImage ( pL, ( *classNames ).getBackgroundClass() );
  726. pixelLabels.addChannel ( pL );
  727. if ( locResult->size() <= 0 )
  728. {
  729. fprintf ( stderr, "WARNING: NO ground truth polygons found for %s !\n",
  730. file.c_str() );
  731. continue;
  732. }
  733. fprintf ( stderr, "SSContext: Collecting pixel examples from localization info: %s\n", file.c_str() );
  734. int depthBoundary = 0;
  735. if ( run3Dseg )
  736. {
  737. depthBoundary = zsizeVec[imgCounter];
  738. }
  739. if ( depthCount < depthBoundary ) continue;
  740. // all image slices collected -> make a 3d image
  741. NICE::MultiChannelImage3DT<double> imgData;
  742. make3DImage ( filelist, imgData );
  743. int xsize = imgData.width();
  744. int ysize = imgData.height();
  745. int zsize = imgData.depth();
  746. amountPixels += xsize * ysize * zsize;
  747. MultiChannelImageT<int> tmpMat ( xsize, ysize, ( uint ) zsize );
  748. labels.push_back ( tmpMat );
  749. nodeIndices.push_back ( MultiChannelImage3DT<unsigned short int> ( xsize, ysize, zsize, nbTrees ) );
  750. nodeIndices[imgCounter].setAll ( 0 );
  751. // MultiChannelImage3DT<double> feats;
  752. // allfeats.push_back ( feats );
  753. int amountRegions;
  754. // convert color to L*a*b, add selected feature channels
  755. addFeatureMaps ( imgData, filelist, amountRegions );
  756. allfeats.push_back(imgData);
  757. if ( useFeat1 )
  758. {
  759. amountRegionpI.push_back ( amountRegions );
  760. rSize.push_back ( vector<int> ( amountRegions, 0 ) );
  761. }
  762. if ( useCategorization )
  763. {
  764. globalCategorFeats.push_back ( new SparseVector() );
  765. classesPerImage.push_back ( map<int,int>() );
  766. }
  767. for ( int x = 0; x < xsize; x++ )
  768. {
  769. for ( int y = 0; y < ysize; y++ )
  770. {
  771. for ( int z = 0; z < zsize; z++ )
  772. {
  773. if ( useFeat1 )
  774. rSize[imgCounter][allfeats[imgCounter] ( x, y, z, rawChannels ) ]++;
  775. if ( run3Dseg )
  776. classno = pixelLabels ( x, y, ( uint ) z );
  777. else
  778. classno = pL.getPixelQuick ( x,y );
  779. labels[imgCounter].set ( x, y, classno, ( uint ) z );
  780. if ( forbidden_classes.find ( classno ) != forbidden_classes.end() )
  781. continue;
  782. labelcounter[classno]++;
  783. if ( useCategorization )
  784. classesPerImage[imgCounter][classno] = 1;
  785. }
  786. }
  787. }
  788. filelist.clear();
  789. pixelLabels.reInit ( 0,0,0 );
  790. depthCount = 0;
  791. imgCounter++;
  792. }
  793. }
  794. int classes = 0;
  795. for ( map<int, int>::const_iterator mapit = labelcounter.begin();
  796. mapit != labelcounter.end(); mapit++ )
  797. {
  798. labelmap[mapit->first] = classes;
  799. labelmapback[classes] = mapit->first;
  800. classes++;
  801. }
  802. ////////////////////////// channel type configuration /////////////////////////
  803. ///////////////////////////////////////////////////////////////////////////////
  804. // Type 0: single pixel & pixel-comparison features on gray value channels
  805. for ( int i = 0; i < rawChannels; i++ )
  806. channelType.push_back ( 0 );
  807. // Type 1: region channel with unsupervised segmentation
  808. int shift = 0;
  809. if ( useFeat1 )
  810. {
  811. channelType.push_back ( 1 );
  812. shift++;
  813. }
  814. // Type 2: rectangular and Haar-like features on gray value integral channels
  815. if ( useFeat2 )
  816. for ( int i = 0; i < rawChannels; i++ )
  817. channelType.push_back ( 2 );
  818. // Type 3: type 2 features on context channels
  819. if ( useFeat3 )
  820. for ( int i = 0; i < classes; i++ )
  821. channelType.push_back ( 3 );
  822. // Type 4: type 0 features on context channels
  823. if ( useFeat4 )
  824. for ( int i = 0; i < classes; i++ )
  825. channelType.push_back ( 4 );
  826. // Type 5: ray features for shape modeling on canny-map
  827. if ( useFeat5 )
  828. for ( int i = 0; i < 24; i++ )
  829. channelType.push_back ( 5 );
  830. // 'amountTypes' sets upper bound for usable feature types
  831. int amountTypes = 6;
  832. channelsPerType = vector<vector<int> > ( amountTypes, vector<int>() );
  833. for ( int i = 0; i < ( int ) channelType.size(); i++ )
  834. {
  835. channelsPerType[channelType[i]].push_back ( i );
  836. }
  837. ///////////////////////////////////////////////////////////////////////////////
  838. ///////////////////////////////////////////////////////////////////////////////
  839. vector<vector<vector<double> > > regionProbs;
  840. if ( useFeat1 )
  841. {
  842. for ( int i = 0; i < imgCounter; i++ )
  843. {
  844. regionProbs.push_back ( vector<vector<double> > ( amountRegionpI[i], vector<double> ( classes, 0.0 ) ) );
  845. }
  846. }
  847. //balancing
  848. a = vector<double> ( classes, 0.0 );
  849. int featcounter = 0;
  850. for ( int iCounter = 0; iCounter < imgCounter; iCounter++ )
  851. {
  852. int xsize = ( int ) nodeIndices[iCounter].width();
  853. int ysize = ( int ) nodeIndices[iCounter].height();
  854. int zsize = ( int ) nodeIndices[iCounter].depth();
  855. for ( int x = 0; x < xsize; x++ )
  856. {
  857. for ( int y = 0; y < ysize; y++ )
  858. {
  859. for ( int z = 0; z < zsize; z++ )
  860. {
  861. featcounter++;
  862. int cn = labels[iCounter] ( x, y, ( uint ) z );
  863. if ( labelmap.find ( cn ) == labelmap.end() )
  864. continue;
  865. a[labelmap[cn]] ++;
  866. }
  867. }
  868. }
  869. }
  870. for ( int i = 0; i < ( int ) a.size(); i++ )
  871. {
  872. a[i] /= ( double ) featcounter;
  873. }
  874. #ifdef VERBOSE
  875. cout << "\nDistribution:" << endl;
  876. for ( int i = 0; i < ( int ) a.size(); i++ )
  877. cout << "class " << i << ": " << a[i] << endl;
  878. #endif
  879. depth = 0;
  880. uniquenumber = 0;
  881. //initialize random forest
  882. for ( int t = 0; t < nbTrees; t++ )
  883. {
  884. vector<TreeNode> singletree;
  885. singletree.push_back ( TreeNode() );
  886. singletree[0].dist = vector<double> ( classes, 0.0 );
  887. singletree[0].depth = depth;
  888. singletree[0].featcounter = amountPixels;
  889. singletree[0].nodeNumber = uniquenumber;
  890. uniquenumber++;
  891. forest.push_back ( singletree );
  892. }
  893. vector<int> startnode ( nbTrees, 0 );
  894. bool noNewSplit = false;
  895. timer.stop();
  896. cout << "\nTime for Pre-Processing: " << timer.getLastAbsolute() << " seconds\n" << endl;
  897. //////////////////////////// train the classifier ///////////////////////////
  898. /////////////////////////////////////////////////////////////////////////////
  899. timer.start();
  900. while ( !noNewSplit && (depth < maxDepth) )
  901. {
  902. depth++;
  903. #ifdef DEBUG
  904. cout << "depth: " << depth << endl;
  905. #endif
  906. noNewSplit = true;
  907. vector<MultiChannelImage3DT<unsigned short int> > lastNodeIndices = nodeIndices;
  908. vector<vector<vector<double> > > lastRegionProbs = regionProbs;
  909. if ( useFeat1 )
  910. for ( int i = 0; i < imgCounter; i++ )
  911. {
  912. int numRegions = (int) regionProbs[i].size();
  913. for ( int r = 0; r < numRegions; r++ )
  914. for ( int c = 0; c < classes; c++ )
  915. regionProbs[i][r][c] = 0.0;
  916. }
  917. // initialize & update context channels
  918. for ( int i = 0; i < imgCounter; i++)
  919. if ( useFeat3 || useFeat4 )
  920. this->updateProbabilityMaps ( nodeIndices[i], allfeats[i], rawChannels + shift );
  921. #ifdef VERBOSE
  922. Timer timerDepth;
  923. timerDepth.start();
  924. #endif
  925. double weight = computeWeight ( depth, maxDepth )
  926. - computeWeight ( depth - 1, maxDepth );
  927. #pragma omp parallel for
  928. // for each tree
  929. for ( int tree = 0; tree < nbTrees; tree++ )
  930. {
  931. const int t = ( int ) forest[tree].size();
  932. const int s = startnode[tree];
  933. startnode[tree] = t;
  934. double bestig;
  935. // for each node
  936. for ( int node = s; node < t; node++ )
  937. {
  938. if ( !forest[tree][node].isleaf && forest[tree][node].left < 0 )
  939. {
  940. // find best split
  941. Operation3D *splitfeat = NULL;
  942. double splitval;
  943. bestig = getBestSplit ( allfeats, lastNodeIndices, labels, node,
  944. splitfeat, splitval, tree, lastRegionProbs );
  945. forest[tree][node].feat = splitfeat;
  946. forest[tree][node].decision = splitval;
  947. // split the node
  948. if ( splitfeat != NULL )
  949. {
  950. noNewSplit = false;
  951. int left;
  952. #pragma omp critical
  953. {
  954. left = forest[tree].size();
  955. forest[tree].push_back ( TreeNode() );
  956. forest[tree].push_back ( TreeNode() );
  957. }
  958. int right = left + 1;
  959. forest[tree][node].left = left;
  960. forest[tree][node].right = right;
  961. forest[tree][left].init( depth, classes, uniquenumber);
  962. int leftu = uniquenumber;
  963. uniquenumber++;
  964. forest[tree][right].init( depth, classes, uniquenumber);
  965. int rightu = uniquenumber;
  966. uniquenumber++;
  967. #pragma omp parallel for
  968. for ( int i = 0; i < imgCounter; i++ )
  969. {
  970. int xsize = nodeIndices[i].width();
  971. int ysize = nodeIndices[i].height();
  972. int zsize = nodeIndices[i].depth();
  973. for ( int x = 0; x < xsize; x++ )
  974. {
  975. for ( int y = 0; y < ysize; y++ )
  976. {
  977. for ( int z = 0; z < zsize; z++ )
  978. {
  979. if ( nodeIndices[i].get ( x, y, z, tree ) == node )
  980. {
  981. // get feature value
  982. Features feat;
  983. feat.feats = &allfeats[i];
  984. feat.rProbs = &lastRegionProbs[i];
  985. double val = 0.0;
  986. val = splitfeat->getVal ( feat, x, y, z );
  987. if ( !isfinite ( val ) ) val = 0.0;
  988. #pragma omp critical
  989. {
  990. int curLabel = labels[i] ( x, y, ( uint ) z );
  991. // traverse to left child
  992. if ( val < splitval )
  993. {
  994. nodeIndices[i].set ( x, y, z, left, tree );
  995. if ( labelmap.find ( curLabel ) != labelmap.end() )
  996. forest[tree][left].dist[labelmap[curLabel]]++;
  997. forest[tree][left].featcounter++;
  998. if ( useCategorization && leftu < shortsize )
  999. ( *globalCategorFeats[i] ) [leftu]+=weight;
  1000. }
  1001. // traverse to right child
  1002. else
  1003. {
  1004. nodeIndices[i].set ( x, y, z, right, tree );
  1005. if ( labelmap.find ( curLabel ) != labelmap.end() )
  1006. forest[tree][right].dist[labelmap[curLabel]]++;
  1007. forest[tree][right].featcounter++;
  1008. if ( useCategorization && rightu < shortsize )
  1009. ( *globalCategorFeats[i] ) [rightu]+=weight;
  1010. }
  1011. }
  1012. }
  1013. }
  1014. }
  1015. }
  1016. }
  1017. // normalize distributions in child leaves
  1018. double lcounter = 0.0, rcounter = 0.0;
  1019. for ( int c = 0; c < (int)forest[tree][left].dist.size(); c++ )
  1020. {
  1021. if ( forbidden_classes.find ( labelmapback[c] ) != forbidden_classes.end() )
  1022. {
  1023. forest[tree][left].dist[c] = 0;
  1024. forest[tree][right].dist[c] = 0;
  1025. }
  1026. else
  1027. {
  1028. forest[tree][left].dist[c] /= a[c];
  1029. lcounter += forest[tree][left].dist[c];
  1030. forest[tree][right].dist[c] /= a[c];
  1031. rcounter += forest[tree][right].dist[c];
  1032. }
  1033. }
  1034. assert ( lcounter > 0 && rcounter > 0 );
  1035. // if ( lcounter <= 0 || rcounter <= 0 )
  1036. // {
  1037. // cout << "lcounter : " << lcounter << " rcounter: " << rcounter << endl;
  1038. // cout << "splitval: " << splitval << " splittype: " << splitfeat->writeInfos() << endl;
  1039. // cout << "bestig: " << bestig << endl;
  1040. // for ( int i = 0; i < imgCounter; i++ )
  1041. // {
  1042. // int xsize = nodeIndices[i].width();
  1043. // int ysize = nodeIndices[i].height();
  1044. // int zsize = nodeIndices[i].depth();
  1045. // int counter = 0;
  1046. // for ( int x = 0; x < xsize; x++ )
  1047. // {
  1048. // for ( int y = 0; y < ysize; y++ )
  1049. // {
  1050. // for ( int z = 0; z < zsize; z++ )
  1051. // {
  1052. // if ( lastNodeIndices[i].get ( x, y, tree ) == node )
  1053. // {
  1054. // if ( ++counter > 30 )
  1055. // break;
  1056. // Features feat;
  1057. // feat.feats = &allfeats[i];
  1058. // feat.rProbs = &lastRegionProbs[i];
  1059. // double val = splitfeat->getVal ( feat, x, y, z );
  1060. // if ( !isfinite ( val ) ) val = 0.0;
  1061. // cout << "splitval: " << splitval << " val: " << val << endl;
  1062. // }
  1063. // }
  1064. // }
  1065. // }
  1066. // }
  1067. // assert ( lcounter > 0 && rcounter > 0 );
  1068. // }
  1069. for ( int c = 0; c < classes; c++ )
  1070. {
  1071. forest[tree][left].dist[c] /= lcounter;
  1072. forest[tree][right].dist[c] /= rcounter;
  1073. }
  1074. }
  1075. else
  1076. {
  1077. forest[tree][node].isleaf = true;
  1078. }
  1079. }
  1080. }
  1081. }
  1082. if ( useFeat1 )
  1083. {
  1084. for ( int i = 0; i < imgCounter; i++ )
  1085. {
  1086. int xsize = nodeIndices[i].width();
  1087. int ysize = nodeIndices[i].height();
  1088. int zsize = nodeIndices[i].depth();
  1089. #pragma omp parallel for
  1090. // set region probability distribution
  1091. for ( int x = 0; x < xsize; x++ )
  1092. {
  1093. for ( int y = 0; y < ysize; y++ )
  1094. {
  1095. for ( int z = 0; z < zsize; z++ )
  1096. {
  1097. for ( int tree = 0; tree < nbTrees; tree++ )
  1098. {
  1099. int node = nodeIndices[i].get ( x, y, z, tree );
  1100. for ( int c = 0; c < classes; c++ )
  1101. {
  1102. int r = (int) ( allfeats[i] ( x, y, z, rawChannels ) );
  1103. regionProbs[i][r][c] += forest[tree][node].dist[c];
  1104. }
  1105. }
  1106. }
  1107. }
  1108. }
  1109. // normalize distribution
  1110. int numRegions = (int) regionProbs[i].size();
  1111. for ( int r = 0; r < numRegions; r++ )
  1112. {
  1113. for ( int c = 0; c < classes; c++ )
  1114. {
  1115. regionProbs[i][r][c] /= ( double ) ( rSize[i][r] );
  1116. }
  1117. }
  1118. }
  1119. }
  1120. if ( firstiteration ) firstiteration = false;
  1121. #ifdef VERBOSE
  1122. timerDepth.stop();
  1123. cout << "Depth " << depth << ": " << timerDepth.getLastAbsolute() << " seconds" <<endl;
  1124. #endif
  1125. lastNodeIndices.clear();
  1126. lastRegionProbs.clear();
  1127. }
  1128. timer.stop();
  1129. cout << "Time for Learning: " << timer.getLastAbsolute() << " seconds\n" << endl;
  1130. //////////////////////// classification using HIK ///////////////////////////
  1131. /////////////////////////////////////////////////////////////////////////////
  1132. if ( useCategorization && fasthik != NULL )
  1133. {
  1134. timer.start();
  1135. uniquenumber = std::min ( shortsize, uniquenumber );
  1136. for ( uint i = 0; i < globalCategorFeats.size(); i++ )
  1137. {
  1138. globalCategorFeats[i]->setDim ( uniquenumber );
  1139. globalCategorFeats[i]->normalize();
  1140. }
  1141. map<int,Vector> ys;
  1142. int cCounter = 0;
  1143. for ( map<int,int>::const_iterator it = labelmap.begin();
  1144. it != labelmap.end(); it++, cCounter++ )
  1145. {
  1146. ys[cCounter] = Vector ( globalCategorFeats.size() );
  1147. for ( int i = 0; i < imgCounter; i++ )
  1148. {
  1149. if ( classesPerImage[i].find ( it->first ) != classesPerImage[i].end() )
  1150. {
  1151. ys[cCounter][i] = 1;
  1152. }
  1153. else
  1154. {
  1155. ys[cCounter][i] = -1;
  1156. }
  1157. }
  1158. }
  1159. fasthik->train( reinterpret_cast<vector<const NICE::SparseVector *>&>(globalCategorFeats), ys);
  1160. timer.stop();
  1161. cerr << "Time for Categorization: " << timer.getLastAbsolute() << " seconds\n" << endl;
  1162. }
  1163. #ifdef VERBOSE
  1164. cout << "\nFEATURE USAGE" << endl;
  1165. cout << "#############\n" << endl;
  1166. // amount of used features per feature type
  1167. std::map<int, int> featTypeCounter;
  1168. for ( int tree = 0; tree < nbTrees; tree++ )
  1169. {
  1170. int t = ( int ) forest[tree].size();
  1171. for ( int node = 0; node < t; node++ )
  1172. {
  1173. if ( !forest[tree][node].isleaf && forest[tree][node].left != -1 )
  1174. {
  1175. featTypeCounter[ forest[tree][node].feat->getFeatType() ] += 1;
  1176. }
  1177. }
  1178. }
  1179. cout << "Types:" << endl;
  1180. for ( map<int, int>::const_iterator it = featTypeCounter.begin(); it != featTypeCounter.end(); it++ )
  1181. cout << it->first << ": " << it->second << endl;
  1182. cout << "\nOperations - All:" << endl;
  1183. // used operations
  1184. vector<int> opOverview ( NBOPERATIONS, 0 );
  1185. // relative use of context vs raw features per tree level
  1186. vector<vector<double> > contextOverview ( maxDepth, vector<double> ( 2, 0.0 ) );
  1187. for ( int tree = 0; tree < nbTrees; tree++ )
  1188. {
  1189. int t = ( int ) forest[tree].size();
  1190. for ( int node = 0; node < t; node++ )
  1191. {
  1192. #ifdef DEBUG
  1193. printf ( "tree[%i]: left: %i, right: %i", node, forest[tree][node].left, forest[tree][node].right );
  1194. #endif
  1195. if ( !forest[tree][node].isleaf && forest[tree][node].left != -1 )
  1196. {
  1197. cout << forest[tree][node].feat->writeInfos() << endl;
  1198. opOverview[ forest[tree][node].feat->getOps() ]++;
  1199. contextOverview[forest[tree][node].depth][ ( int ) forest[tree][node].feat->getContext() ]++;
  1200. }
  1201. #ifdef DEBUG
  1202. for ( int d = 0; d < ( int ) forest[tree][node].dist.size(); d++ )
  1203. {
  1204. cout << " " << forest[tree][node].dist[d];
  1205. }
  1206. cout << endl;
  1207. #endif
  1208. }
  1209. }
  1210. // amount of used features per operation type
  1211. cout << "\nOperations - Summary:" << endl;
  1212. for ( int t = 0; t < ( int ) opOverview.size(); t++ )
  1213. {
  1214. cout << "Ops " << t << ": " << opOverview[ t ] << endl;
  1215. }
  1216. // ratio of used context features per depth level
  1217. cout << "\nContext-Ratio:" << endl;
  1218. for ( int d = 0; d < maxDepth; d++ )
  1219. {
  1220. double sum = contextOverview[d][0] + contextOverview[d][1];
  1221. if ( sum == 0 )
  1222. sum = 1;
  1223. contextOverview[d][0] /= sum;
  1224. contextOverview[d][1] /= sum;
  1225. cout << "Depth [" << d+1 << "] Normal: " << contextOverview[d][0] << " Context: " << contextOverview[d][1] << endl;
  1226. }
  1227. #endif
  1228. }
  1229. void SemSegContextTree3D::addFeatureMaps (
  1230. NICE::MultiChannelImage3DT<double> &imgData,
  1231. const vector<string> &filelist,
  1232. int &amountRegions )
  1233. {
  1234. int xsize = imgData.width();
  1235. int ysize = imgData.height();
  1236. int zsize = imgData.depth();
  1237. amountRegions = 0;
  1238. // RGB to Lab
  1239. if ( imagetype == IMAGETYPE_RGB )
  1240. {
  1241. for ( int z = 0; z < zsize; z++ )
  1242. for ( int y = 0; y < ysize; y++ )
  1243. for ( int x = 0; x < xsize; x++ )
  1244. {
  1245. double R, G, B, X, Y, Z, L, a, b;
  1246. R = ( double )imgData.get( x, y, z, 0 ) / 255.0;
  1247. G = ( double )imgData.get( x, y, z, 1 ) / 255.0;
  1248. B = ( double )imgData.get( x, y, z, 2 ) / 255.0;
  1249. if ( useAltTristimulus )
  1250. {
  1251. ColorConversion::ccRGBtoXYZ( R, G, B, &X, &Y, &Z, 4 );
  1252. ColorConversion::ccXYZtoCIE_Lab( X, Y, Z, &L, &a, &b, 4 );
  1253. }
  1254. else
  1255. {
  1256. ColorConversion::ccRGBtoXYZ( R, G, B, &X, &Y, &Z, 0 );
  1257. ColorConversion::ccXYZtoCIE_Lab( X, Y, Z, &L, &a, &b, 0 );
  1258. }
  1259. imgData.set( x, y, z, L, 0 );
  1260. imgData.set( x, y, z, a, 1 );
  1261. imgData.set( x, y, z, b, 2 );
  1262. }
  1263. }
  1264. // Gradient layers
  1265. if ( useGradient )
  1266. {
  1267. int currentsize = imgData.channels();
  1268. imgData.addChannel ( 2*currentsize );
  1269. for ( int z = 0; z < zsize; z++ )
  1270. for ( int c = 0; c < currentsize; c++ )
  1271. {
  1272. ImageT<double> tmp = imgData.getChannelT(z, c);
  1273. ImageT<double> sobX( xsize, ysize );
  1274. ImageT<double> sobY( xsize, ysize );
  1275. NICE::FilterT<double, double, double>::sobelX ( tmp, sobX );
  1276. NICE::FilterT<double, double, double>::sobelY ( tmp, sobY );
  1277. for ( int y = 0; y < ysize; y++ )
  1278. for ( int x = 0; x < xsize; x++ )
  1279. {
  1280. imgData.set( x, y, z, sobX.getPixelQuick(x,y), c+currentsize );
  1281. imgData.set( x, y, z, sobY.getPixelQuick(x,y), c+(currentsize*2) );
  1282. }
  1283. }
  1284. }
  1285. // Weijer color names
  1286. if ( useWeijer )
  1287. {
  1288. if ( imagetype == IMAGETYPE_RGB )
  1289. {
  1290. int currentsize = imgData.channels();
  1291. imgData.addChannel ( 11 );
  1292. for ( int z = 0; z < zsize; z++ )
  1293. {
  1294. NICE::ColorImage img = imgData.getColor ( z );
  1295. NICE::MultiChannelImageT<double> cfeats;
  1296. lfcw->getFeats ( img, cfeats );
  1297. for ( int c = 0; c < cfeats.channels(); c++)
  1298. for ( int y = 0; y < ysize; y++ )
  1299. for ( int x = 0; x < xsize; x++ )
  1300. imgData.set(x, y, z, cfeats.get(x,y,(uint)c), c+currentsize);
  1301. }
  1302. }
  1303. else
  1304. {
  1305. cerr << "Can't compute weijer features of a grayscale image." << endl;
  1306. }
  1307. }
  1308. // arbitrary additional layer as image
  1309. if ( useAdditionalLayer )
  1310. {
  1311. int currentsize = imgData.channels();
  1312. imgData.addChannel ( 1 );
  1313. for ( int z = 0; z < zsize; z++ )
  1314. {
  1315. vector<string> list;
  1316. StringTools::split ( filelist[z], '/', list );
  1317. string layerPath = StringTools::trim ( filelist[z], list.back() ) + "addlayer/" + list.back();
  1318. NICE::Image layer ( layerPath );
  1319. for ( int y = 0; y < ysize; y++ )
  1320. for ( int x = 0; x < xsize; x++ )
  1321. imgData.set(x, y, z, layer.getPixelQuick(x,y), currentsize);
  1322. }
  1323. }
  1324. // read the geometric cues produced by Hoiem et al.
  1325. if ( useHoiemFeatures )
  1326. {
  1327. string hoiemDirectory = conf->gS ( "Features", "hoiem_directory" );
  1328. // we could also give the following set as a config option
  1329. string hoiemClasses_s = "sky 000 090-045 090-090 090-135 090 090-por 090-sol";
  1330. vector<string> hoiemClasses;
  1331. StringTools::split ( hoiemClasses_s, ' ', hoiemClasses );
  1332. int currentsize = imgData.channels();
  1333. imgData.addChannel ( hoiemClasses.size() );
  1334. for ( int z = 0; z < zsize; z++ )
  1335. {
  1336. FileName fn ( filelist[z] );
  1337. fn.removeExtension();
  1338. FileName fnBase = fn.extractFileName();
  1339. for ( vector<string>::const_iterator i = hoiemClasses.begin(); i != hoiemClasses.end(); i++, currentsize++ )
  1340. {
  1341. string hoiemClass = *i;
  1342. FileName fnConfidenceImage ( hoiemDirectory + fnBase.str() + "_c_" + hoiemClass + ".png" );
  1343. if ( ! fnConfidenceImage.fileExists() )
  1344. {
  1345. fthrow ( Exception, "Unable to read the Hoiem geometric confidence image: " << fnConfidenceImage.str() << " (original image is " << filelist[z] << ")" );
  1346. }
  1347. else
  1348. {
  1349. Image confidenceImage ( fnConfidenceImage.str() );
  1350. if ( confidenceImage.width() != xsize || confidenceImage.height() != ysize )
  1351. {
  1352. fthrow ( Exception, "The size of the geometric confidence image does not match with the original image size: " << fnConfidenceImage.str() );
  1353. }
  1354. // copy standard image to double image
  1355. for ( int y = 0 ; y < confidenceImage.height(); y++ )
  1356. for ( int x = 0 ; x < confidenceImage.width(); x++ )
  1357. imgData ( x, y, z, currentsize ) = ( double ) confidenceImage ( x, y );
  1358. currentsize++;
  1359. }
  1360. }
  1361. }
  1362. }
  1363. // region feature (unsupervised segmentation)
  1364. int shift = 0;
  1365. if ( useFeat1 )
  1366. {
  1367. shift = 1;
  1368. MultiChannelImageT<int> regions;
  1369. regions.reInit( xsize, ysize, zsize );
  1370. amountRegions = segmentation->segRegions ( imgData, regions, imagetype );
  1371. int currentsize = imgData.channels();
  1372. imgData.addChannel ( 1 );
  1373. for ( int z = 0; z < ( int ) regions.channels(); z++ )
  1374. for ( int y = 0; y < regions.height(); y++ )
  1375. for ( int x = 0; x < regions.width(); x++ )
  1376. imgData.set ( x, y, z, regions ( x, y, ( uint ) z ), currentsize );
  1377. }
  1378. // intergal images of raw channels
  1379. if ( useFeat2 )
  1380. {
  1381. imgData.addChannel ( rawChannels );
  1382. #pragma omp parallel for
  1383. for ( int i = 0; i < rawChannels; i++ )
  1384. {
  1385. int corg = i;
  1386. int cint = i + rawChannels + shift;
  1387. for ( int z = 0; z < zsize; z++ )
  1388. for ( int y = 0; y < ysize; y++ )
  1389. for ( int x = 0; x < xsize; x++ )
  1390. imgData ( x, y, z, cint ) = imgData ( x, y, z, corg );
  1391. imgData.calcIntegral ( cint );
  1392. }
  1393. }
  1394. int classes = classNames->numClasses();
  1395. if ( useFeat3 )
  1396. imgData.addChannel ( classes );
  1397. if ( useFeat4 )
  1398. imgData.addChannel ( classes );
  1399. if ( useFeat5 )
  1400. {
  1401. imgData.addChannel ( 24 );
  1402. this->computeRayFeatImage( imgData, imgData.channels()-24);
  1403. }
  1404. }
  1405. void SemSegContextTree3D::classify (
  1406. const std::vector<std::string> & filelist,
  1407. NICE::MultiChannelImageT<double> & segresult,
  1408. NICE::MultiChannelImage3DT<double> & probabilities )
  1409. {
  1410. ///////////////////////// build MCI3DT from files ///////////////////////////
  1411. /////////////////////////////////////////////////////////////////////////////
  1412. NICE::MultiChannelImage3DT<double> imgData;
  1413. this->make3DImage( filelist, imgData );
  1414. int xsize = imgData.width();
  1415. int ysize = imgData.height();
  1416. int zsize = imgData.depth();
  1417. ////////////////////////// initialize variables /////////////////////////////
  1418. /////////////////////////////////////////////////////////////////////////////
  1419. firstiteration = true;
  1420. depth = 0;
  1421. Timer timer;
  1422. timer.start();
  1423. // classes occurred during training step
  1424. int classes = labelmapback.size();
  1425. // classes defined in config file
  1426. int numClasses = classNames->numClasses();
  1427. // class probabilities by pixel
  1428. probabilities.reInit ( xsize, ysize, zsize, numClasses );
  1429. probabilities.setAll ( 0 );
  1430. // class probabilities by region
  1431. vector<vector<double> > regionProbs;
  1432. // affiliation: pixel <-> (tree,node)
  1433. MultiChannelImage3DT<unsigned short int> nodeIndices ( xsize, ysize, zsize, nbTrees );
  1434. nodeIndices.setAll ( 0 );
  1435. // for categorization
  1436. SparseVector *globalCategorFeat;
  1437. globalCategorFeat = new SparseVector();
  1438. /////////////////////////// get feature values //////////////////////////////
  1439. /////////////////////////////////////////////////////////////////////////////
  1440. // Basic Features
  1441. int amountRegions;
  1442. addFeatureMaps ( imgData, filelist, amountRegions );
  1443. vector<int> rSize;
  1444. int shift = 0;
  1445. if ( useFeat1 )
  1446. {
  1447. shift = 1;
  1448. regionProbs = vector<vector<double> > ( amountRegions, vector<double> ( classes, 0.0 ) );
  1449. rSize = vector<int> ( amountRegions, 0 );
  1450. for ( int z = 0; z < zsize; z++ )
  1451. {
  1452. for ( int y = 0; y < ysize; y++ )
  1453. {
  1454. for ( int x = 0; x < xsize; x++ )
  1455. {
  1456. rSize[imgData ( x, y, z, rawChannels ) ]++;
  1457. }
  1458. }
  1459. }
  1460. }
  1461. ////////////////// traverse image example through trees /////////////////////
  1462. /////////////////////////////////////////////////////////////////////////////
  1463. bool noNewSplit = false;
  1464. for ( int d = 0; d < maxDepth && !noNewSplit; d++ )
  1465. {
  1466. depth++;
  1467. vector<vector<double> > lastRegionProbs = regionProbs;
  1468. if ( useFeat1 )
  1469. {
  1470. int numRegions = ( int ) regionProbs.size();
  1471. for ( int r = 0; r < numRegions; r++ )
  1472. for ( int c = 0; c < classes; c++ )
  1473. regionProbs[r][c] = 0.0;
  1474. }
  1475. if ( depth < maxDepth )
  1476. {
  1477. int firstChannel = rawChannels + shift;
  1478. if ( useFeat3 || useFeat4 )
  1479. this->updateProbabilityMaps ( nodeIndices, imgData, firstChannel );
  1480. }
  1481. double weight = computeWeight ( depth, maxDepth )
  1482. - computeWeight ( depth - 1, maxDepth );
  1483. noNewSplit = true;
  1484. int tree;
  1485. #pragma omp parallel for private(tree)
  1486. for ( tree = 0; tree < nbTrees; tree++ )
  1487. {
  1488. for ( int x = 0; x < xsize; x++ )
  1489. {
  1490. for ( int y = 0; y < ysize; y++ )
  1491. {
  1492. for ( int z = 0; z < zsize; z++ )
  1493. {
  1494. int node = nodeIndices.get ( x, y, z, tree );
  1495. if ( forest[tree][node].left > 0 )
  1496. {
  1497. noNewSplit = false;
  1498. Features feat;
  1499. feat.feats = &imgData;
  1500. feat.rProbs = &lastRegionProbs;
  1501. double val = forest[tree][node].feat->getVal ( feat, x, y, z );
  1502. if ( !isfinite ( val ) ) val = 0.0;
  1503. // traverse to left child
  1504. if ( val < forest[tree][node].decision )
  1505. {
  1506. int left = forest[tree][node].left;
  1507. nodeIndices.set ( x, y, z, left, tree );
  1508. #pragma omp critical
  1509. {
  1510. if ( fasthik != NULL
  1511. && useCategorization
  1512. && forest[tree][left].nodeNumber < uniquenumber )
  1513. ( *globalCategorFeat ) [forest[tree][left].nodeNumber] += weight;
  1514. }
  1515. }
  1516. // traverse to right child
  1517. else
  1518. {
  1519. int right = forest[tree][node].right;
  1520. nodeIndices.set ( x, y, z, right, tree );
  1521. #pragma omp critical
  1522. {
  1523. if ( fasthik != NULL
  1524. && useCategorization
  1525. && forest[tree][right].nodeNumber < uniquenumber )
  1526. ( *globalCategorFeat ) [forest[tree][right].nodeNumber] += weight;
  1527. }
  1528. }
  1529. }
  1530. }
  1531. }
  1532. }
  1533. }
  1534. if ( useFeat1 )
  1535. {
  1536. int xsize = nodeIndices.width();
  1537. int ysize = nodeIndices.height();
  1538. int zsize = nodeIndices.depth();
  1539. #pragma omp parallel for
  1540. for ( int x = 0; x < xsize; x++ )
  1541. {
  1542. for ( int y = 0; y < ysize; y++ )
  1543. {
  1544. for ( int z = 0; z < zsize; z++ )
  1545. {
  1546. for ( int tree = 0; tree < nbTrees; tree++ )
  1547. {
  1548. int node = nodeIndices.get ( x, y, z, tree );
  1549. for ( uint c = 0; c < forest[tree][node].dist.size(); c++ )
  1550. {
  1551. int r = (int) imgData ( x, y, z, rawChannels );
  1552. regionProbs[r][c] += forest[tree][node].dist[c];
  1553. }
  1554. }
  1555. }
  1556. }
  1557. }
  1558. int numRegions = (int) regionProbs.size();
  1559. for ( int r = 0; r < numRegions; r++ )
  1560. {
  1561. for ( int c = 0; c < (int) classes; c++ )
  1562. {
  1563. regionProbs[r][c] /= ( double ) ( rSize[r] );
  1564. }
  1565. }
  1566. }
  1567. if ( (depth < maxDepth) && firstiteration ) firstiteration = false;
  1568. }
  1569. vector<int> classesInImg;
  1570. if ( useCategorization )
  1571. {
  1572. if ( cndir != "" )
  1573. {
  1574. for ( int z = 0; z < zsize; z++ )
  1575. {
  1576. vector< string > list;
  1577. StringTools::split ( filelist[z], '/', list );
  1578. string orgname = list.back();
  1579. ifstream infile ( ( cndir + "/" + orgname + ".dat" ).c_str() );
  1580. while ( !infile.eof() && infile.good() )
  1581. {
  1582. int tmp;
  1583. infile >> tmp;
  1584. assert ( tmp >= 0 && tmp < numClasses );
  1585. classesInImg.push_back ( tmp );
  1586. }
  1587. }
  1588. }
  1589. else
  1590. {
  1591. globalCategorFeat->setDim ( uniquenumber );
  1592. globalCategorFeat->normalize();
  1593. ClassificationResult cr = fasthik->classify( globalCategorFeat);
  1594. for ( uint i = 0; i < ( uint ) classes; i++ )
  1595. {
  1596. cerr << cr.scores[i] << " ";
  1597. if ( cr.scores[i] > 0.0/*-0.3*/ )
  1598. {
  1599. classesInImg.push_back ( i );
  1600. }
  1601. }
  1602. }
  1603. cerr << "amount of classes: " << classes << " used classes: " << classesInImg.size() << endl;
  1604. }
  1605. if ( classesInImg.size() == 0 )
  1606. {
  1607. for ( uint i = 0; i < ( uint ) classes; i++ )
  1608. {
  1609. classesInImg.push_back ( i );
  1610. }
  1611. }
  1612. // final labeling step
  1613. if ( pixelWiseLabeling )
  1614. {
  1615. for ( int x = 0; x < xsize; x++ )
  1616. {
  1617. for ( int y = 0; y < ysize; y++ )
  1618. {
  1619. for ( int z = 0; z < zsize; z++ )
  1620. {
  1621. //TODO by nodes instead of pixel?
  1622. double maxProb = - numeric_limits<double>::max();
  1623. int maxClass = 0;
  1624. for ( uint c = 0; c < classesInImg.size(); c++ )
  1625. {
  1626. int i = classesInImg[c];
  1627. double curProb = getMeanProb ( x, y, z, i, nodeIndices );
  1628. probabilities ( x, y, z, labelmapback[i] ) = curProb;
  1629. if ( curProb > maxProb )
  1630. {
  1631. maxProb = curProb;
  1632. maxClass = labelmapback[i];
  1633. }
  1634. }
  1635. assert(maxProb <= 1);
  1636. // copy pixel labeling into segresults (output)
  1637. segresult.set ( x, y, maxClass, ( uint ) z );
  1638. }
  1639. }
  1640. }
  1641. #ifdef VISUALIZE
  1642. getProbabilityMap( probabilities );
  1643. #endif
  1644. }
  1645. else
  1646. {
  1647. // labeling by region
  1648. NICE::MultiChannelImageT<int> regions;
  1649. int xsize = imgData.width();
  1650. int ysize = imgData.height();
  1651. int zsize = imgData.depth();
  1652. regions.reInit ( xsize, ysize, zsize );
  1653. if ( useFeat1 )
  1654. {
  1655. int rchannel = -1;
  1656. for ( uint i = 0; i < channelType.size(); i++ )
  1657. {
  1658. if ( channelType[i] == 1 )
  1659. {
  1660. rchannel = i;
  1661. break;
  1662. }
  1663. }
  1664. assert ( rchannel > -1 );
  1665. for ( int z = 0; z < zsize; z++ )
  1666. {
  1667. for ( int y = 0; y < ysize; y++ )
  1668. {
  1669. for ( int x = 0; x < xsize; x++ )
  1670. {
  1671. regions.set ( x, y, imgData ( x, y, z, rchannel ), ( uint ) z );
  1672. }
  1673. }
  1674. }
  1675. }
  1676. else
  1677. {
  1678. amountRegions = segmentation->segRegions ( imgData, regions, imagetype );
  1679. #ifdef DEBUG
  1680. for ( unsigned int z = 0; z < ( uint ) zsize; z++ )
  1681. {
  1682. NICE::Matrix regmask;
  1683. NICE::ColorImage colorimg ( xsize, ysize );
  1684. NICE::ColorImage marked ( xsize, ysize );
  1685. regmask.resize ( xsize, ysize );
  1686. for ( int y = 0; y < ysize; y++ )
  1687. {
  1688. for ( int x = 0; x < xsize; x++ )
  1689. {
  1690. regmask ( x,y ) = regions ( x,y,z );
  1691. colorimg.setPixelQuick ( x, y, 0, imgData.get ( x,y,z,0 ) );
  1692. colorimg.setPixelQuick ( x, y, 1, imgData.get ( x,y,z,0 ) );
  1693. colorimg.setPixelQuick ( x, y, 2, imgData.get ( x,y,z,0 ) );
  1694. }
  1695. }
  1696. vector<int> colorvals;
  1697. colorvals.push_back ( 255 );
  1698. colorvals.push_back ( 0 );
  1699. colorvals.push_back ( 0 );
  1700. segmentation->markContours ( colorimg, regmask, colorvals, marked );
  1701. std::vector<string> list;
  1702. StringTools::split ( filelist[z], '/', list );
  1703. string savePath = StringTools::trim ( filelist[z], list.back() ) + "marked/" + list.back();
  1704. marked.write ( savePath );
  1705. }
  1706. #endif
  1707. }
  1708. regionProbs.clear();
  1709. regionProbs = vector<vector<double> > ( amountRegions, vector<double> ( classes, 0.0 ) );
  1710. vector<int> bestlabels ( amountRegions, labelmapback[classesInImg[0]] );
  1711. for ( int z = 0; z < zsize; z++ )
  1712. {
  1713. for ( int y = 0; y < ysize; y++ )
  1714. {
  1715. for ( int x = 0; x < xsize; x++ )
  1716. {
  1717. int r = regions ( x, y, ( uint ) z );
  1718. for ( uint i = 0; i < classesInImg.size(); i++ )
  1719. {
  1720. int c = classesInImg[i];
  1721. // get mean voting of all trees
  1722. regionProbs[r][c] += getMeanProb ( x, y, z, c, nodeIndices );
  1723. }
  1724. }
  1725. }
  1726. }
  1727. for ( int r = 0; r < amountRegions; r++ )
  1728. {
  1729. double maxProb = regionProbs[r][classesInImg[0]];
  1730. bestlabels[r] = classesInImg[0];
  1731. for ( int c = 1; c < classes; c++ )
  1732. {
  1733. if ( maxProb < regionProbs[r][c] )
  1734. {
  1735. maxProb = regionProbs[r][c];
  1736. bestlabels[r] = c;
  1737. }
  1738. }
  1739. bestlabels[r] = labelmapback[bestlabels[r]];
  1740. }
  1741. // copy region labeling into segresults (output)
  1742. for ( int z = 0; z < zsize; z++ )
  1743. {
  1744. for ( int y = 0; y < ysize; y++ )
  1745. {
  1746. for ( int x = 0; x < xsize; x++ )
  1747. {
  1748. segresult.set ( x, y, bestlabels[regions ( x,y, ( uint ) z ) ], ( uint ) z );
  1749. }
  1750. }
  1751. }
  1752. #ifdef WRITEREGIONS
  1753. for ( int z = 0; z < zsize; z++ )
  1754. {
  1755. RegionGraph rg;
  1756. NICE::ColorImage img ( xsize,ysize );
  1757. if ( imagetype == IMAGETYPE_RGB )
  1758. {
  1759. img = imgData.getColor ( z );
  1760. }
  1761. else
  1762. {
  1763. NICE::Image gray = imgData.getChannel ( z );
  1764. for ( int y = 0; y < ysize; y++ )
  1765. {
  1766. for ( int x = 0; x < xsize; x++ )
  1767. {
  1768. int val = gray.getPixelQuick ( x,y );
  1769. img.setPixelQuick ( x, y, val, val, val );
  1770. }
  1771. }
  1772. }
  1773. Matrix regions_tmp ( xsize,ysize );
  1774. for ( int y = 0; y < ysize; y++ )
  1775. {
  1776. for ( int x = 0; x < xsize; x++ )
  1777. {
  1778. regions_tmp ( x,y ) = regions ( x,y, ( uint ) z );
  1779. }
  1780. }
  1781. segmentation->getGraphRepresentation ( img, regions_tmp, rg );
  1782. for ( uint pos = 0; pos < regionProbs.size(); pos++ )
  1783. {
  1784. rg[pos]->setProbs ( regionProbs[pos] );
  1785. }
  1786. std::string s;
  1787. std::stringstream out;
  1788. std::vector< std::string > list;
  1789. StringTools::split ( filelist[z], '/', list );
  1790. out << "rgout/" << list.back() << ".graph";
  1791. string writefile = out.str();
  1792. rg.write ( writefile );
  1793. }
  1794. #endif
  1795. }
  1796. timer.stop();
  1797. cout << "\nTime for Classification: " << timer.getLastAbsolute() << endl;
  1798. // CLEANING UP
  1799. // TODO: operations in "forest"
  1800. while( !ops.empty() )
  1801. {
  1802. vector<Operation3D*> &tops = ops.back();
  1803. while ( !tops.empty() )
  1804. tops.pop_back();
  1805. ops.pop_back();
  1806. }
  1807. delete globalCategorFeat;
  1808. }
  1809. void SemSegContextTree3D::store ( std::ostream & os, int format ) const
  1810. {
  1811. os.precision ( numeric_limits<double>::digits10 + 1 );
  1812. os << nbTrees << endl;
  1813. classnames.store ( os );
  1814. map<int, int>::const_iterator it;
  1815. os << labelmap.size() << endl;
  1816. for ( it = labelmap.begin() ; it != labelmap.end(); it++ )
  1817. os << ( *it ).first << " " << ( *it ).second << endl;
  1818. os << labelmapback.size() << endl;
  1819. for ( it = labelmapback.begin() ; it != labelmapback.end(); it++ )
  1820. os << ( *it ).first << " " << ( *it ).second << endl;
  1821. int trees = forest.size();
  1822. os << trees << endl;
  1823. for ( int t = 0; t < trees; t++ )
  1824. {
  1825. int nodes = forest[t].size();
  1826. os << nodes << endl;
  1827. for ( int n = 0; n < nodes; n++ )
  1828. {
  1829. os << forest[t][n].left << " " << forest[t][n].right << " " << forest[t][n].decision << " " << forest[t][n].isleaf << " " << forest[t][n].depth << " " << forest[t][n].featcounter << " " << forest[t][n].nodeNumber << endl;
  1830. os << forest[t][n].dist << endl;
  1831. if ( forest[t][n].feat == NULL )
  1832. os << -1 << endl;
  1833. else
  1834. {
  1835. os << forest[t][n].feat->getOps() << endl;
  1836. forest[t][n].feat->store ( os );
  1837. }
  1838. }
  1839. }
  1840. os << channelType.size() << endl;
  1841. for ( int i = 0; i < ( int ) channelType.size(); i++ )
  1842. {
  1843. os << channelType[i] << " ";
  1844. }
  1845. os << endl;
  1846. os << rawChannels << endl;
  1847. os << uniquenumber << endl;
  1848. }
  1849. void SemSegContextTree3D::restore ( std::istream & is, int format )
  1850. {
  1851. is >> nbTrees;
  1852. classnames.restore ( is );
  1853. int lsize;
  1854. is >> lsize;
  1855. labelmap.clear();
  1856. for ( int l = 0; l < lsize; l++ )
  1857. {
  1858. int first, second;
  1859. is >> first;
  1860. is >> second;
  1861. labelmap[first] = second;
  1862. }
  1863. is >> lsize;
  1864. labelmapback.clear();
  1865. for ( int l = 0; l < lsize; l++ )
  1866. {
  1867. int first, second;
  1868. is >> first;
  1869. is >> second;
  1870. labelmapback[first] = second;
  1871. }
  1872. int trees;
  1873. is >> trees;
  1874. forest.clear();
  1875. for ( int t = 0; t < trees; t++ )
  1876. {
  1877. vector<TreeNode> tmptree;
  1878. forest.push_back ( tmptree );
  1879. int nodes;
  1880. is >> nodes;
  1881. for ( int n = 0; n < nodes; n++ )
  1882. {
  1883. TreeNode tmpnode;
  1884. forest[t].push_back ( tmpnode );
  1885. is >> forest[t][n].left;
  1886. is >> forest[t][n].right;
  1887. is >> forest[t][n].decision;
  1888. is >> forest[t][n].isleaf;
  1889. is >> forest[t][n].depth;
  1890. is >> forest[t][n].featcounter;
  1891. is >> forest[t][n].nodeNumber;
  1892. is >> forest[t][n].dist;
  1893. int feattype;
  1894. is >> feattype;
  1895. assert ( feattype < NBOPERATIONS );
  1896. forest[t][n].feat = NULL;
  1897. if ( feattype >= 0 )
  1898. {
  1899. for ( uint o = 0; o < ops.size(); o++ )
  1900. {
  1901. for ( uint o2 = 0; o2 < ops[o].size(); o2++ )
  1902. {
  1903. if ( forest[t][n].feat == NULL )
  1904. {
  1905. if ( ops[o][o2]->getOps() == feattype )
  1906. {
  1907. forest[t][n].feat = ops[o][o2]->clone();
  1908. break;
  1909. }
  1910. }
  1911. }
  1912. }
  1913. assert ( forest[t][n].feat != NULL );
  1914. forest[t][n].feat->restore ( is );
  1915. }
  1916. }
  1917. }
  1918. channelType.clear();
  1919. int ctsize;
  1920. is >> ctsize;
  1921. for ( int i = 0; i < ctsize; i++ )
  1922. {
  1923. int tmp;
  1924. is >> tmp;
  1925. switch (tmp)
  1926. {
  1927. case 0: useFeat0 = true; break;
  1928. case 1: useFeat1 = true; break;
  1929. case 2: useFeat2 = true; break;
  1930. case 3: useFeat3 = true; break;
  1931. case 4: useFeat4 = true; break;
  1932. case 5: useFeat5 = true; break;
  1933. }
  1934. channelType.push_back ( tmp );
  1935. }
  1936. // integralMap is deprecated but kept in RESTORE
  1937. // for downwards compatibility!
  1938. // std::vector<std::pair<int, int> > integralMap;
  1939. // integralMap.clear();
  1940. // int iMapSize;
  1941. // is >> iMapSize;
  1942. // for ( int i = 0; i < iMapSize; i++ )
  1943. // {
  1944. // int first;
  1945. // int second;
  1946. // is >> first;
  1947. // is >> second;
  1948. // integralMap.push_back ( pair<int, int> ( first, second ) );
  1949. // }
  1950. is >> rawChannels;
  1951. is >> uniquenumber;
  1952. }