MultiChannelImage3DT.tcc 21 KB

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