MultiChannelImage3DT.tcc 23 KB

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