MultiChannelImage3DT.tcc 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943
  1. #include <iostream>
  2. #include <assert.h>
  3. #include <stdio.h>
  4. #include <vector>
  5. namespace NICE {
  6. template<class P>
  7. MultiChannelImage3DT<P>::MultiChannelImage3DT( int _xsize, int _ysize, int _zsize, uint _numChannels)
  8. {
  9. // data = NULL;
  10. numChannels = 0;
  11. xsize = 0;
  12. ysize = 0;
  13. zsize = 0;
  14. reInit( _xsize, _ysize, _zsize, _numChannels);
  15. }
  16. template<class P>
  17. MultiChannelImage3DT<P>::MultiChannelImage3DT()
  18. {
  19. xsize = 0;
  20. ysize = 0;
  21. zsize = 0;
  22. numChannels = 0;
  23. // data = NULL;
  24. }
  25. template<class P>
  26. P & MultiChannelImage3DT<P>::operator() (int x, int y, int z, uint channel)
  27. {
  28. assert( channel < numChannels );
  29. assert(( x < xsize ) && ( x >= 0 ) );
  30. assert(( y < ysize ) && ( y >= 0 ) );
  31. assert(( z < zsize ) && ( z >= 0 ) );
  32. assert( data[channel] != NULL );
  33. return data[channel][x + y*xsize + z*xsize*ysize];
  34. }
  35. template<class P>
  36. MultiChannelImageT<P> MultiChannelImage3DT<P>::operator[] (uint z)
  37. {
  38. MultiChannelImageT<P> img;
  39. for( int c = 0; c < numChannels; c++ )
  40. {
  41. P * datatmp = data[c];
  42. ImageT<P> tmp ( &datatmp[z*(xsize*ysize)], xsize, ysize, xsize*sizeof(P), GrayColorImageCommonImplementation::shallowCopy );
  43. img.addChannel(tmp);
  44. }
  45. return img;
  46. }
  47. template<class P>
  48. MultiChannelImage3DT<P>& MultiChannelImage3DT<P>::operator=(const MultiChannelImage3DT<P>& orig)
  49. {
  50. if( xsize == orig.xsize
  51. && ysize == orig.ysize
  52. && zsize == orig.zsize
  53. && numChannels == orig.numChannels)
  54. {
  55. int iMemSize = xsize*ysize*zsize;
  56. for(int c=0; c < numChannels; ++c)
  57. {
  58. std::copy(orig.data[c], orig.data[c]+iMemSize, data[c]);
  59. }
  60. }
  61. else
  62. {
  63. freeData();
  64. xsize = orig.xsize;
  65. ysize = orig.ysize;
  66. zsize = orig.zsize;
  67. numChannels = orig.numChannels;
  68. int iMemSize = xsize*ysize*zsize;
  69. for(int c = 0; c < numChannels; ++c)
  70. {
  71. P *t_newData = new P [iMemSize];
  72. std::copy(orig.data[c], orig.data[c]+iMemSize, t_newData);
  73. data.push_back( t_newData );
  74. }
  75. }
  76. return *this;
  77. }
  78. template<class P>
  79. MultiChannelImage3DT<P>::MultiChannelImage3DT( const MultiChannelImage3DT<P>& p )
  80. {
  81. xsize = p.xsize;
  82. ysize = p.ysize;
  83. zsize = p.zsize;
  84. numChannels = p.numChannels;
  85. int iMemSize = xsize*ysize*zsize;
  86. for(int c=0; c < numChannels; ++c)
  87. {
  88. P *t_newData = new P [iMemSize];
  89. std::copy(p.data[c], p.data[c]+iMemSize, t_newData);
  90. data.push_back( t_newData );
  91. }
  92. /* if(p.data != NULL)
  93. data = new P *[numChannels];
  94. else
  95. data = NULL;
  96. for ( int c = 0; c < ( int )numChannels; c++ )
  97. {
  98. if ( p.data[c] == NULL )
  99. {
  100. data[c] = NULL;
  101. }
  102. else
  103. {
  104. data[c] = new P [xsize*ysize*zsize];
  105. for ( int x = 0; x < xsize*ysize*zsize; x++ )
  106. {
  107. data[c][x] = p.data[c][x];
  108. }
  109. }
  110. }
  111. */
  112. }
  113. template<class P>
  114. void MultiChannelImage3DT<P>::addChannel( int newChans )
  115. {
  116. assert(xsize > 0);
  117. assert(ysize > 0);
  118. assert(zsize > 0);
  119. for (int i = 0; i < newChans; i++ )
  120. {
  121. this->data.push_back( new P [xsize*ysize*zsize] );
  122. }
  123. numChannels = this->data.size();
  124. }
  125. template<class P>
  126. template<class SrcP>
  127. void MultiChannelImage3DT<P>::addChannel( const NICE::ImageT<SrcP> &newImg )
  128. {
  129. int oldchan = numChannels;
  130. if(this->xsize > 0)
  131. {
  132. assert(newImg.width() == this->width() && newImg.height() == this->height());
  133. addChannel(1);
  134. }
  135. else
  136. {
  137. reInit( newImg.width(), newImg.height(), 1, 1 );
  138. }
  139. for ( int y = 0; y < ysize; y++ )
  140. for ( int x = 0; x < ysize; x++ )
  141. data[oldchan][x + y*xsize] = (P) newImg.getPixelQuick(x,y);
  142. }
  143. template<class P>
  144. template<class SrcP>
  145. void MultiChannelImage3DT<P>::addChannel(const NICE::MultiChannelImageT<SrcP> &newMCImg)
  146. {
  147. int oldchan = numChannels;
  148. if(this->xsize > 0)
  149. {
  150. assert(newMCImg.width() == this->width() && newMCImg.height() == this->height());
  151. assert(newMCImg.channels() == this->zsize);
  152. addChannel(1);
  153. }
  154. else
  155. {
  156. reInit( newMCImg.width(), newMCImg.height(), newMCImg.channels(), 1 );
  157. }
  158. for(int z = 0; z < this->zsize; z++)
  159. {
  160. //NICE::ImageT<SrcP> newImg = newMCImg[z];
  161. for(int y = 0; y < this->ysize; y++)
  162. {
  163. for(int x = 0; x < this->xsize; x++)
  164. {
  165. //data[oldchan][x + y*xsize + z*xsize*ysize] = (P)newImg(x,y);
  166. data[oldchan][x + y*xsize + z*xsize*ysize] = (P)newMCImg.get(x,y,(unsigned int)z);
  167. }
  168. }
  169. }
  170. }
  171. template<class P>
  172. template<class SrcP>
  173. void MultiChannelImage3DT<P>::addChannel(const NICE::MultiChannelImage3DT<SrcP> &newImg)
  174. {
  175. int oldchan = numChannels;
  176. if(numChannels > 0)
  177. {
  178. assert(newImg.width() == this->width() && newImg.height() == this->height() && newImg.depth() == this->depth());
  179. addChannel(newImg.channels());
  180. }
  181. else
  182. {
  183. reInit( newImg.width(), newImg.height(), newImg.depth(), newImg.channels() );
  184. }
  185. int chanNI = 0;
  186. for(int c = oldchan; c < (int)numChannels; c++, chanNI++)
  187. {
  188. int val = 0;
  189. for(int z = 0; z < this->zsize; z++)
  190. {
  191. for(int y = 0; y < this->ysize; y++)
  192. {
  193. for(int x = 0; x < this->xsize; x++, val++)
  194. {
  195. data[c][val] = newImg.get(x,y,z,chanNI);
  196. }
  197. }
  198. }
  199. }
  200. }
  201. template<class P>
  202. template<class SrcP>
  203. void MultiChannelImage3DT<P>::addChannelReferences(const NICE::MultiChannelImage3DT<SrcP> &newImg)
  204. {
  205. if(numChannels == 0)
  206. {
  207. xsize = newImg.width();
  208. ysize = newImg.height();
  209. zsize = newImg.depth();
  210. }
  211. else
  212. {
  213. if( !( newImg.width() == this->width()
  214. && newImg.height() == this->height()
  215. && newImg.depth() == this->depth() ) )
  216. {
  217. throw( " channelwise dimensions don't fit! Abort this crazy mixing of channels" );
  218. return;
  219. }
  220. }
  221. // add channel deep data references
  222. // -> wrap this MultiChannelImage3D around the other data
  223. const std::vector< P* > vecDataChannelPtrs = newImg.getDataPointer();
  224. for( int c=0; c < vecDataChannelPtrs.size() ; ++c)
  225. data.push_back( vecDataChannelPtrs[c] );
  226. numChannels = data.size();
  227. }
  228. template<class P>
  229. MultiChannelImage3DT<P>::~MultiChannelImage3DT()
  230. {
  231. freeData();
  232. }
  233. template<class P>
  234. void MultiChannelImage3DT<P>::freeData()
  235. {
  236. if ( !data.empty())
  237. {
  238. for ( int i = 0 ; i < (int)data.size() ; i++ )
  239. if ( data[i] != NULL )
  240. delete [] data[i];
  241. data.clear();
  242. }
  243. }
  244. template<class P>
  245. void MultiChannelImage3DT<P>::freeShallowData()
  246. {
  247. this->data.clear();
  248. }
  249. template<class P>
  250. void MultiChannelImage3DT<P>::reInit( int _xsize, int _ysize, int _zsize, int _numChannels )
  251. {
  252. freeData();
  253. xsize = _xsize;
  254. ysize = _ysize;
  255. zsize = _zsize;
  256. numChannels = _numChannels;
  257. this->addChannel(numChannels);
  258. }
  259. template<class P>
  260. template<class SrcP>
  261. void MultiChannelImage3DT<P>::reInitFrom( const MultiChannelImage3DT<SrcP> & src )
  262. {
  263. this->reInit(src.width(),src.height(), src.depth(), src.channels() );
  264. }
  265. template<class P>
  266. P MultiChannelImage3DT<P>::get( int x, int y, int z, uint channel ) const
  267. {
  268. assert( channel < numChannels );
  269. assert(( x < xsize ) && ( x >= 0 ) );
  270. assert(( y < ysize ) && ( y >= 0 ) );
  271. assert(( z < zsize ) && ( z >= 0 ) );
  272. assert( data[channel] != NULL );
  273. return data[channel][x + y*xsize + z*xsize*ysize];
  274. }
  275. template<class P>
  276. std::vector<P*> MultiChannelImage3DT<P>::getDataPointer() const
  277. {
  278. return data;
  279. }
  280. template<class P>
  281. void MultiChannelImage3DT<P>::set( int x, int y, int z, P val, uint channel )
  282. {
  283. assert( channel < numChannels );
  284. assert(( x < xsize ) && ( x >= 0 ) );
  285. assert(( y < ysize ) && ( y >= 0 ) );
  286. assert(( z < zsize ) && ( z >= 0 ) );
  287. assert( data[channel] != NULL );
  288. data[channel][x + y*xsize + z*xsize*ysize] = val;
  289. }
  290. template<class P>
  291. void MultiChannelImage3DT<P>::set( P val, uint channel )
  292. {
  293. assert( channel < numChannels );
  294. assert( data[channel] != NULL );
  295. for ( int k = 0 ; k < xsize*ysize*zsize ; k++ )
  296. data[channel][k] = val;
  297. }
  298. template<class P>
  299. void MultiChannelImage3DT<P>::setAll( P val )
  300. {
  301. for ( uint channel = 0 ; channel < numChannels ; channel++ )
  302. if ( data[channel] != NULL )
  303. set( val, channel );
  304. }
  305. template<class P>
  306. void MultiChannelImage3DT<P>::statistics( P & min, P & max, uint channel ) const
  307. {
  308. assert( channel < numChannels );
  309. P val = 0;
  310. for ( long k = 0 ; k < xsize*ysize*zsize ; k++ )
  311. {
  312. val = data [channel][k];
  313. if (( k == 0 ) || ( val > max ) ) max = val;
  314. if (( k == 0 ) || ( val < min ) ) min = val;
  315. }
  316. assert(finite(max));
  317. assert(finite(min));
  318. }
  319. template<class P>
  320. void MultiChannelImage3DT<P>::correctShading( uint channel ) const
  321. {
  322. assert( channel < numChannels );
  323. // some sort of correction trick hardly understandable :-)
  324. std::vector<double> meanVals;
  325. for( int z = 0; z < zsize; z++ )
  326. {
  327. double sumVal = 0;
  328. for( int y = 0; y < ysize; y++ )
  329. {
  330. for( int x = 0; x < xsize; x++ )
  331. {
  332. sumVal += data [channel][x + y*xsize + z*xsize*ysize];
  333. }
  334. }
  335. sumVal /= (xsize*ysize);
  336. meanVals.push_back( sumVal );
  337. }
  338. P newMax = std::numeric_limits<P>::min();
  339. const short int maxVal = 255;
  340. for ( int z = 0; z < zsize; z++ )
  341. {
  342. for ( int y = 0; y < ysize; y++ )
  343. {
  344. for ( int x = 0; x < xsize; x++ )
  345. {
  346. P tmp = data [channel][x + y*xsize + z*xsize*ysize];
  347. double newVal = maxVal * ( (double) tmp / meanVals[z] );
  348. if ( ( P ) newVal > newMax )
  349. newMax = ( P ) newVal;
  350. data [channel][x + y*xsize + z*xsize*ysize] = newVal;
  351. }
  352. }
  353. }
  354. for ( long k = 0 ; k < xsize*ysize*zsize ; k++ )
  355. {
  356. data [channel][k] = data [channel][k] / newMax * maxVal;
  357. }
  358. }
  359. template<class P>
  360. void MultiChannelImage3DT<P>::equalizeHistogram( uint channel ) const
  361. {
  362. assert(channel < numChannels );
  363. for( int z = 0; z < zsize; z++ )
  364. {
  365. NICE::Image img = getChannel(z, channel );
  366. NICE::Histogram hist(img,0,255,256);
  367. NICE::IntVector *histVec = NULL;
  368. histVec = hist.cumulative();
  369. for ( int i = 0; i < (int)histVec->size(); i++)
  370. {
  371. histVec->set(i, histVec->get(i) * 255 / (double)histVec->get(histVec->size()-1));
  372. }
  373. for ( int y = 0; y < ysize; y++ )
  374. {
  375. for ( int x = 0; x < xsize; x++ )
  376. {
  377. data [channel][x + y*xsize + z*xsize*ysize] = histVec->get( img.getPixel(x,y) );
  378. }
  379. }
  380. delete histVec;
  381. }
  382. }
  383. template<class P>
  384. Image MultiChannelImage3DT<P>::getChannel( int z, uint channel ) const
  385. {
  386. assert( channel < numChannels );
  387. NICE::Image img(xsize, ysize);
  388. convertToGrey( img, z, channel, true );
  389. return img;
  390. }
  391. template<class P>
  392. ImageT<P> MultiChannelImage3DT<P>::getChannelT( int z, uint channel ) const
  393. {
  394. assert( channel < numChannels );
  395. // P min, max;
  396. // statistics ( min, max, channel );
  397. // fprintf (stderr, "MultiChannelImage3DT<>::showChannel: max %f min %f\n", (double)max, (double)min );
  398. NICE::ImageT<P> img(xsize,ysize);
  399. long k = 0;
  400. for ( int y = 0; y < ysize; y++ )
  401. for( int x = 0; x < xsize; x++, k++ )
  402. {
  403. img.setPixel( x, y, data[channel][z*xsize*ysize + k] );
  404. }
  405. return img;
  406. }
  407. template<class P>
  408. ImageT<P> MultiChannelImage3DT<P>::getXSlice ( int x, uint channel ) const
  409. {
  410. assert( channel < numChannels );
  411. NICE::ImageT<P> img(zsize, ysize);
  412. for ( int y = 0; y < ysize; y++ )
  413. for ( int z = 0; z < zsize; z++ )
  414. img.setPixel( z, y, data[channel][z*xsize*ysize + y*xsize + x]);
  415. return img;
  416. }
  417. /** convert to ice image */
  418. template<class P>
  419. void MultiChannelImage3DT<P>::convertToGrey( NICE::Image & img, int z, uint channel, bool normalize ) const
  420. {
  421. assert( channel < numChannels );
  422. P min, max;
  423. if ( normalize )
  424. statistics( min, max, channel );
  425. bool skip_assignment = false;
  426. img.resize( xsize, ysize );
  427. if ( normalize )
  428. if ( max - min < std::numeric_limits<double>::min() )
  429. {
  430. fprintf( stderr, "MultiChannelImage3DT<>::showChannel: max %f min %f\n", ( double )max, ( double )min );
  431. img.set( max );
  432. skip_assignment = true;
  433. fprintf( stderr, "MultiChannelImage3DT<>::showChannel: image is uniform! (%f)\n", ( double )max );
  434. }
  435. if ( ! skip_assignment )
  436. {
  437. long k = 0;
  438. for ( int y = 0 ; y < ysize; y++ )
  439. {
  440. for ( int x = 0 ; x < xsize ; x++, k++ )
  441. {
  442. if ( normalize )
  443. {
  444. img.setPixel( x, y, ( int )(( data[channel][z*xsize*ysize + k] - min ) * 255 / ( max - min ) ) );
  445. }
  446. else
  447. {
  448. img.setPixel( x, y, ( int )( data[channel][z*xsize*ysize + k] ) );
  449. }
  450. }
  451. }
  452. }
  453. }
  454. template<class P>
  455. void MultiChannelImage3DT<P>::convertToColor( NICE::ColorImage & img, int z, const int chan1, const int chan2, const int chan3) const
  456. {
  457. assert( chan1 < numChannels && chan2 < numChannels && chan3 < numChannels);
  458. img.resize( xsize, ysize );
  459. long k = 0;
  460. for ( int y = 0 ; y < ysize; y++ )
  461. {
  462. for ( int x = 0 ; x < xsize ; x++, k++ )
  463. {
  464. img.setPixel( x, y, 0, ( int )( data[chan1][z*xsize*ysize + k] ) );
  465. img.setPixel( x, y, 1, ( int )( data[chan2][z*xsize*ysize + k] ) );
  466. img.setPixel( x, y, 2, ( int )( data[chan3][z*xsize*ysize + k] ) );
  467. }
  468. }
  469. }
  470. template<class P>
  471. NICE::MultiChannelImageT<P> MultiChannelImage3DT<P>::getColorMCI(int z) const
  472. {
  473. assert( z < zsize );
  474. assert( numChannels >= 3 );
  475. NICE::MultiChannelImageT<P> img( xsize, ysize, 3 );
  476. long k = 0;
  477. for ( int y = 0 ; y < ysize; y++ )
  478. for ( int x = 0; x < xsize; x++ )
  479. {
  480. img.set( x, y, 0, ( int )( data[0][z*xsize*ysize + k] ) );
  481. img.set( x, y, 1, ( int )( data[1][z*xsize*ysize + k] ) );
  482. img.set( x, y, 2, ( int )( data[2][z*xsize*ysize + k] ) );
  483. }
  484. return img;
  485. }
  486. template<class P>
  487. NICE::ColorImage MultiChannelImage3DT<P>::getColor(int z) const
  488. {
  489. assert( z < zsize );
  490. assert( numChannels >= 3 );
  491. NICE::ColorImage img( xsize, ysize );
  492. long k = 0;
  493. for ( int y = 0 ; y < ysize; y++ )
  494. {
  495. for ( int x = 0 ; x < xsize ; x++, k++ )
  496. {
  497. img.setPixel( x, y, 0, ( int )( data[0][z*xsize*ysize + k] ) );
  498. img.setPixel( x, y, 1, ( int )( data[1][z*xsize*ysize + k] ) );
  499. img.setPixel( x, y, 2, ( int )( data[2][z*xsize*ysize + k] ) );
  500. }
  501. }
  502. //showImage(img);
  503. //getchar();
  504. return img;
  505. }
  506. template<class P>
  507. NICE::ColorImage MultiChannelImage3DT<P>::getColorImageFromChannels(int z, int channel0, int channel1, int channel2) const
  508. {
  509. assert( z < zsize );
  510. assert( numChannels >= std::max( std::max(channel0,channel1),channel2 ) );
  511. NICE::ColorImage img( xsize, ysize );
  512. long k = 0;
  513. for ( int y = 0 ; y < ysize; y++ )
  514. {
  515. for ( int x = 0 ; x < xsize ; x++, k++ )
  516. {
  517. img.setPixel( x, y, 0, ( int )( data[channel0][z*xsize*ysize + k] ) );
  518. img.setPixel( x, y, 1, ( int )( data[channel1][z*xsize*ysize + k] ) );
  519. img.setPixel( x, y, 2, ( int )( data[channel2][z*xsize*ysize + k] ) );
  520. }
  521. }
  522. //showImage(img);
  523. //getchar();
  524. return img;
  525. }
  526. template<class P>
  527. void MultiChannelImage3DT<P>::calcIntegral( uint channel )
  528. {
  529. assert( channel < numChannels );
  530. assert( data[channel] != NULL );
  531. P *integralImage = data[channel];
  532. /** first column **/
  533. int k = xsize;
  534. for ( int y = 1 ; y < ysize; y++, k += xsize )
  535. integralImage[k] += integralImage[k-xsize];
  536. /** first row **/
  537. k = 1;
  538. for ( int x = 1 ; x < xsize; x++, k++ )
  539. integralImage[k] += integralImage[k-1];
  540. /** first stack (depth) **/
  541. k = xsize * ysize;
  542. for ( int z = 1 ; z < zsize; z++, k += (xsize*ysize) )
  543. integralImage[k] += integralImage[k-(xsize*ysize)];
  544. /** x-y plane **/
  545. k = xsize + 1;
  546. for ( int y = 1 ; y < ysize ; y++, k++ )
  547. for ( int x = 1 ; x < xsize ; x++, k++ )
  548. {
  549. integralImage[k] += integralImage[k-1];
  550. integralImage[k] += integralImage[k - xsize];
  551. integralImage[k] -= integralImage[k - xsize - 1];
  552. }
  553. /** y-z plane **/
  554. k = xsize*ysize + xsize;
  555. for ( int z = 1 ; z < zsize ; z++, k+=xsize )
  556. for ( int y = 1 ; y < ysize ; y++, k+=xsize )
  557. {
  558. integralImage[k] += integralImage[k-(xsize*ysize)];
  559. integralImage[k] += integralImage[k - xsize];
  560. integralImage[k] -= integralImage[k - xsize - (xsize*ysize)];
  561. }
  562. /** x-z plane **/
  563. k = xsize*ysize + 1;
  564. for ( int z = 1 ; z < zsize ; z++, k+=((xsize*ysize)-(xsize-1)) )
  565. for ( int x = 1 ; x < xsize ; x++, k++ )
  566. {
  567. integralImage[k] += integralImage[k-1];
  568. integralImage[k] += integralImage[k - (xsize*ysize)];
  569. integralImage[k] -= integralImage[k - (xsize*ysize) - 1];
  570. }
  571. /** all other pixels **/
  572. k = xsize*ysize + xsize + 1;
  573. for ( int z = 1 ; z < zsize ; z++, k+= xsize )
  574. {
  575. for ( int y = 1 ; y < ysize ; y++, k++ )
  576. {
  577. for ( int x = 1 ; x < xsize ; x++, k++ )
  578. {
  579. integralImage[k] += integralImage[k - (xsize*ysize)];
  580. integralImage[k] += integralImage[k - xsize];
  581. integralImage[k] += integralImage[k - 1];
  582. integralImage[k] += integralImage[k - (xsize*ysize) - xsize - 1];
  583. integralImage[k] -= integralImage[k - (xsize*ysize) - xsize];
  584. integralImage[k] -= integralImage[k - (xsize*ysize) - 1];
  585. integralImage[k] -= integralImage[k - xsize - 1];
  586. }
  587. }
  588. }
  589. }
  590. template<class P>
  591. void MultiChannelImage3DT<P>::calcVariance( uint srcchan, uint tarchan )
  592. {
  593. assert( srcchan < tarchan );
  594. assert( tarchan < numChannels );
  595. assert( data[srcchan] != NULL );
  596. assert( data[tarchan] != NULL );
  597. uint windowsize = 3;
  598. int win = (windowsize-1)/2;
  599. for ( int z = 0; z < zsize; z++ )
  600. {
  601. for ( int y = 0; y < ysize; y++ )
  602. {
  603. for ( int x = 0; x < xsize; x++ )
  604. {
  605. int meansum = 0;
  606. for ( int u = -win; u <= win; u++ )
  607. {
  608. for ( int v = -win; v <= win; v++ )
  609. {
  610. for ( int w = -win; w <= win; w++)
  611. {
  612. int u_tmp = u;
  613. int v_tmp = v;
  614. int w_tmp = w;
  615. if ( (x+u<0) || (x+u>=xsize) )
  616. u_tmp = -u_tmp;
  617. if ( (y+v<0) || (y+v>=ysize) )
  618. v_tmp = -v_tmp;
  619. if ( (z+w<0) || (z+w>=zsize) )
  620. w_tmp = -w_tmp;
  621. meansum += get( x+u_tmp, y+v_tmp, z+w_tmp, srcchan );
  622. }
  623. }
  624. }
  625. meansum /= (windowsize*windowsize*windowsize);
  626. unsigned long varsum = 0;
  627. for ( int u = -win; u <= win; u++ )
  628. {
  629. for ( int v = -win; v <= win; v++ )
  630. {
  631. for ( int w = -win; w <= win; w++)
  632. {
  633. int u_tmp = u;
  634. int v_tmp = v;
  635. int w_tmp = w;
  636. if ( (x+u<0) || (x+u>=xsize) )
  637. u_tmp = -u_tmp;
  638. if ( (y+v<0) || (y+v>=ysize) )
  639. v_tmp = -v_tmp;
  640. if ( (z+w<0) || (z+w>=zsize) )
  641. w_tmp = -w_tmp;
  642. long sdev = (get( x+u_tmp, y+v_tmp, z+w_tmp, srcchan ) - meansum );
  643. varsum += (sdev*sdev);
  644. }
  645. }
  646. }
  647. varsum /= (windowsize*windowsize+windowsize)-1;
  648. set( x, y, z, varsum, tarchan );
  649. }
  650. }
  651. }
  652. }
  653. template<class P>
  654. P MultiChannelImage3DT<P>::getIntegralValue(int ulfx, int ulfy, int ulfz, int lrbx, int lrby, int lrbz, int channel) const
  655. {
  656. ulfx = std::max(ulfx-1, -1);
  657. ulfx = std::min(ulfx, xsize-1);
  658. ulfy = std::max(ulfy-1, -1);
  659. ulfy = std::min(ulfy, ysize-1);
  660. ulfz = std::max(ulfz-1, -1);
  661. ulfz = std::min(ulfz, zsize-1);
  662. lrbx = std::max(lrbx, 0);
  663. lrbx = std::min(lrbx, xsize-1);
  664. lrby = std::max(lrby, 0);
  665. lrby = std::min(lrby, ysize-1);
  666. lrbz = std::max(lrbz, 0);
  667. lrbz = std::min(lrbz, zsize-1);
  668. double val1, val2, val3, val4, val5, val6, val7, val8;
  669. val1 = get(lrbx, lrby, lrbz, channel);
  670. if( ulfz > -1 )
  671. val2 = get(lrbx, lrby, ulfz, channel);
  672. else
  673. val2 = 0;
  674. if( ulfx > -1 )
  675. val3 = get(ulfx, lrby, lrbz, channel);
  676. else
  677. val3 = 0;
  678. if( ulfx > -1 && ulfz > -1 )
  679. val4 = get(ulfx, lrby, ulfz, channel);
  680. else
  681. val4 = 0;
  682. if( ulfy > -1 )
  683. val5 = get(lrbx, ulfy, lrbz, channel);
  684. else
  685. val5 = 0;
  686. if( ulfy > -1 && ulfz > -1 )
  687. val6 = get(lrbx, ulfy, ulfz, channel);
  688. else
  689. val6 = 0;
  690. if( ulfx > -1 && ulfy > -1 )
  691. val7 = get(ulfx, ulfy, lrbz, channel);
  692. else
  693. val7 = 0;
  694. if(ulfx > -1 && ulfy > -1 && ulfz > -1)
  695. val8 = get(ulfx, ulfy, ulfz, channel);
  696. else
  697. val8 = 0;
  698. P volume = abs((lrbx-ulfx)*(lrby-ulfy)*(lrbz-ulfz));
  699. P val = val1 - val2 - val3 + val4 - ( val5 - val6 - val7 + val8 );
  700. if (volume != 0)
  701. return val/volume;
  702. else
  703. return 0.0;
  704. }
  705. template<class P>
  706. void MultiChannelImage3DT<P>::store( std::string filename ) const
  707. {
  708. // simple raw format
  709. FILE *f = fopen( filename.c_str(), "w" );
  710. if ( f == NULL ) {
  711. fprintf( stderr, "MultiChannelImage3DT::store: error writing to %s\n", filename.c_str() );
  712. exit( -1 );
  713. }
  714. fwrite( &xsize, sizeof( int ), 1, f );
  715. fwrite( &ysize, sizeof( int ), 1, f );
  716. fwrite( &zsize, sizeof( int ), 1, f );
  717. fwrite( &numChannels, sizeof( uint ), 1, f );
  718. for ( uint channel = 0 ; channel < numChannels ; channel++ )
  719. {
  720. assert( data[channel] != NULL );
  721. fwrite( data[channel], sizeof( P ), xsize*ysize*zsize, f );
  722. }
  723. fclose( f );
  724. }
  725. template<class P>
  726. void MultiChannelImage3DT<P>::restore( std::string filename )
  727. {
  728. // simple raw format
  729. FILE *f = fopen( filename.c_str(), "r" );
  730. if ( f == NULL ) {
  731. fprintf( stderr, "MultiChannelImage3DT::store: error reading from %s\n", filename.c_str() );
  732. exit( -1 );
  733. }
  734. fread( &xsize, sizeof( int ), 1, f );
  735. fread( &ysize, sizeof( int ), 1, f );
  736. fread( &zsize, sizeof( int ), 1, f );
  737. fread( &numChannels, sizeof( uint ), 1, f );
  738. if ( numChannels > 0 ) {
  739. reInit( xsize, ysize, zsize, numChannels );
  740. for ( uint channel = 0 ; channel < numChannels ; channel++ )
  741. {
  742. assert( data[channel] != NULL );
  743. fread( data[channel], sizeof( P ), xsize*ysize*zsize, f );
  744. }
  745. } else {
  746. freeData();
  747. //data = NULL;
  748. }
  749. fclose( f );
  750. }
  751. template<class P>
  752. int MultiChannelImage3DT<P>::width() const
  753. {
  754. return xsize;
  755. }
  756. template<class P>
  757. int MultiChannelImage3DT<P>::height() const
  758. {
  759. return ysize;
  760. }
  761. template<class P>
  762. int MultiChannelImage3DT<P>::depth() const
  763. {
  764. return zsize;
  765. }
  766. template<class P>
  767. int MultiChannelImage3DT<P>::channels() const
  768. {
  769. return ( int )numChannels;
  770. }
  771. template<class P>
  772. int MultiChannelImage3DT<P>::getPixelInt( int x, int y, int z, int channel ) const
  773. {
  774. throw( "this type is not implemented\n" );
  775. return -1;
  776. }
  777. template<class P>
  778. double MultiChannelImage3DT<P>::getPixelFloat( int x, int y, int z, int channel ) const
  779. {
  780. throw( "this type is not implemented\n" );
  781. return -1.0;
  782. }
  783. template<class P>
  784. void MultiChannelImage3DT<P>::setPixelInt( int x, int y, int z, int channel, int pixel )
  785. {
  786. throw( "this type is not implemented\n" );
  787. }
  788. template<class P>
  789. void MultiChannelImage3DT<P>::setPixelFloat( int x, int y, int z, int channel, double pixel )
  790. {
  791. throw( "this type is not implemented\n" );
  792. }
  793. #define SET_FUNCS_PROTO_MACRO3D(MYTYPE) \
  794. template<>\
  795. int MultiChannelImage3DT<MYTYPE>::getPixelInt(int x, int y, int z, int channel) const;\
  796. template<>\
  797. double MultiChannelImage3DT<MYTYPE>::getPixelFloat(int x, int y, int z, int channel) const;\
  798. template<>\
  799. void MultiChannelImage3DT<MYTYPE>::setPixelInt(int x, int y, int z, int channel, int pixel);\
  800. template<>\
  801. void MultiChannelImage3DT<MYTYPE>::setPixelFloat(int x, int y, int z, int channel, double pixel);
  802. SET_FUNCS_PROTO_MACRO3D( double )
  803. SET_FUNCS_PROTO_MACRO3D( int )
  804. SET_FUNCS_PROTO_MACRO3D( long int )
  805. SET_FUNCS_PROTO_MACRO3D( float )
  806. SET_FUNCS_PROTO_MACRO3D( unsigned int )
  807. } // namespace