MultiChannelImage3DT.tcc 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864
  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. void MultiChannelImage3DT<P>::calcIntegral( uint channel )
  459. {
  460. assert( channel < numChannels );
  461. assert( data[channel] != NULL );
  462. P *integralImage = data[channel];
  463. /** first column **/
  464. int k = xsize;
  465. for ( int y = 1 ; y < ysize; y++, k += xsize )
  466. integralImage[k] += integralImage[k-xsize];
  467. /** first row **/
  468. k = 1;
  469. for ( int x = 1 ; x < xsize; x++, k++ )
  470. integralImage[k] += integralImage[k-1];
  471. /** first stack (depth) **/
  472. k = xsize * ysize;
  473. for ( int z = 1 ; z < zsize; z++, k += (xsize*ysize) )
  474. integralImage[k] += integralImage[k-(xsize*ysize)];
  475. /** x-y plane **/
  476. k = xsize + 1;
  477. for ( int y = 1 ; y < ysize ; y++, k++ )
  478. for ( int x = 1 ; x < xsize ; x++, k++ )
  479. {
  480. integralImage[k] += integralImage[k-1];
  481. integralImage[k] += integralImage[k - xsize];
  482. integralImage[k] -= integralImage[k - xsize - 1];
  483. }
  484. /** y-z plane **/
  485. k = xsize*ysize + xsize;
  486. for ( int z = 1 ; z < zsize ; z++, k+=xsize )
  487. for ( int y = 1 ; y < zsize ; y++, k+=xsize )
  488. {
  489. integralImage[k] += integralImage[k-(xsize*ysize)];
  490. integralImage[k] += integralImage[k - xsize];
  491. integralImage[k] -= integralImage[k - xsize - (xsize*ysize)];
  492. }
  493. /** x-z plane **/
  494. k = xsize*ysize + 1;
  495. for ( int z = 1 ; z < zsize ; z++, k+=((xsize*ysize)-(xsize-1)) )
  496. for ( int x = 1 ; x < xsize ; x++, k++ )
  497. {
  498. integralImage[k] += integralImage[k-1];
  499. integralImage[k] += integralImage[k - (xsize*ysize)];
  500. integralImage[k] -= integralImage[k - (xsize*ysize) - 1];
  501. }
  502. /** all other pixels **/
  503. k = xsize*ysize + xsize + 1;
  504. for ( int z = 1 ; z < zsize ; z++, k+= xsize )
  505. {
  506. for ( int y = 1 ; y < ysize ; y++, k++ )
  507. {
  508. for ( int x = 1 ; x < xsize ; x++, k++ )
  509. {
  510. integralImage[k] += integralImage[k - (xsize*ysize)];
  511. integralImage[k] += integralImage[k - xsize];
  512. integralImage[k] += integralImage[k - 1];
  513. integralImage[k] += integralImage[k - (xsize*ysize) - xsize - 1];
  514. integralImage[k] -= integralImage[k - (xsize*ysize) - xsize];
  515. integralImage[k] -= integralImage[k - (xsize*ysize) - 1];
  516. integralImage[k] -= integralImage[k - xsize - 1];
  517. }
  518. }
  519. }
  520. }
  521. template<class P>
  522. void MultiChannelImage3DT<P>::calcVariance( uint srcchan, uint tarchan )
  523. {
  524. assert( srcchan < tarchan );
  525. assert( tarchan < numChannels );
  526. assert( data[srcchan] != NULL );
  527. assert( data[tarchan] != NULL );
  528. uint windowsize = 3;
  529. int win = (windowsize-1)/2;
  530. for ( int z = 0; z < zsize; z++ )
  531. {
  532. for ( int y = 0; y < ysize; y++ )
  533. {
  534. for ( int x = 0; x < xsize; x++ )
  535. {
  536. int meansum = 0;
  537. for ( int u = -win; u <= win; u++ )
  538. {
  539. for ( int v = -win; v <= win; v++ )
  540. {
  541. for ( int w = -win; w <= win; w++)
  542. {
  543. int u_tmp = u;
  544. int v_tmp = v;
  545. int w_tmp = w;
  546. if ( (x+u<0) || (x+u>=xsize) )
  547. u_tmp = -u_tmp;
  548. if ( (y+v<0) || (y+v>=ysize) )
  549. v_tmp = -v_tmp;
  550. if ( (z+w<0) || (z+w>=zsize) )
  551. w_tmp = -w_tmp;
  552. meansum += get( x+u_tmp, y+v_tmp, z+w_tmp, srcchan );
  553. }
  554. }
  555. }
  556. meansum /= (windowsize*windowsize*windowsize);
  557. unsigned long varsum = 0;
  558. for ( int u = -win; u <= win; u++ )
  559. {
  560. for ( int v = -win; v <= win; v++ )
  561. {
  562. for ( int w = -win; w <= win; w++)
  563. {
  564. int u_tmp = u;
  565. int v_tmp = v;
  566. int w_tmp = w;
  567. if ( (x+u<0) || (x+u>=xsize) )
  568. u_tmp = -u_tmp;
  569. if ( (y+v<0) || (y+v>=ysize) )
  570. v_tmp = -v_tmp;
  571. if ( (z+w<0) || (z+w>=zsize) )
  572. w_tmp = -w_tmp;
  573. long sdev = (get( x+u_tmp, y+v_tmp, z+w_tmp, srcchan ) - meansum );
  574. varsum += (sdev*sdev);
  575. }
  576. }
  577. }
  578. varsum /= (windowsize*windowsize+windowsize)-1;
  579. set( x, y, z, varsum, tarchan );
  580. }
  581. }
  582. }
  583. }
  584. template<class P>
  585. P MultiChannelImage3DT<P>::getIntegralValue(int ulfx, int ulfy, int ulfz, int lrbx, int lrby, int lrbz, int channel)
  586. {
  587. ulfx = std::max(ulfx-1, -1);
  588. ulfx = std::min(ulfx, xsize-1);
  589. ulfy = std::max(ulfy-1, -1);
  590. ulfy = std::min(ulfy, ysize-1);
  591. ulfz = std::max(ulfz-1, -1);
  592. ulfz = std::min(ulfz, zsize-1);
  593. lrbx = std::max(lrbx, 0);
  594. lrbx = std::min(lrbx, xsize-1);
  595. lrby = std::max(lrby, 0);
  596. lrby = std::min(lrby, ysize-1);
  597. lrbz = std::max(lrbz, 0);
  598. lrbz = std::min(lrbz, zsize-1);
  599. double val1, val2, val3, val4, val5, val6, val7, val8;
  600. val1 = get(lrbx, lrby, lrbz, channel);
  601. if( ulfz > -1 )
  602. val2 = get(lrbx, lrby, ulfz, channel);
  603. else
  604. val2 = 0;
  605. if( ulfx > -1 )
  606. val3 = get(ulfx, lrby, lrbz, channel);
  607. else
  608. val3 = 0;
  609. if( ulfx > -1 && ulfz > -1 )
  610. val4 = get(ulfx, lrby, ulfz, channel);
  611. else
  612. val4 = 0;
  613. if( ulfy > -1 )
  614. val5 = get(lrbx, ulfy, lrbz, channel);
  615. else
  616. val5 = 0;
  617. if( ulfy > -1 && ulfz > -1 )
  618. val6 = get(lrbx, ulfy, ulfz, channel);
  619. else
  620. val6 = 0;
  621. if( ulfx > -1 && ulfy > -1 )
  622. val7 = get(ulfx, ulfy, lrbz, channel);
  623. else
  624. val7 = 0;
  625. if(ulfx > -1 && ulfy > -1 && ulfz > -1)
  626. val8 = get(ulfx, ulfy, ulfz, channel);
  627. else
  628. val8 = 0;
  629. P volume = abs((lrbx-ulfx)*(lrby-ulfy)*(lrbz-ulfz));
  630. P val = val1 - val2 - val3 + val4 - ( val5 - val6 - val7 + val8 );
  631. if (volume != 0)
  632. return val/volume;
  633. else
  634. return 0.0;
  635. }
  636. template<class P>
  637. void MultiChannelImage3DT<P>::store( std::string filename ) const
  638. {
  639. // simple raw format
  640. FILE *f = fopen( filename.c_str(), "w" );
  641. if ( f == NULL ) {
  642. fprintf( stderr, "MultiChannelImage3DT::store: error writing to %s\n", filename.c_str() );
  643. exit( -1 );
  644. }
  645. fwrite( &xsize, sizeof( int ), 1, f );
  646. fwrite( &ysize, sizeof( int ), 1, f );
  647. fwrite( &zsize, sizeof( int ), 1, f );
  648. fwrite( &numChannels, sizeof( uint ), 1, f );
  649. for ( uint channel = 0 ; channel < numChannels ; channel++ )
  650. {
  651. assert( data[channel] != NULL );
  652. fwrite( data[channel], sizeof( P ), xsize*ysize*zsize, f );
  653. }
  654. fclose( f );
  655. }
  656. template<class P>
  657. void MultiChannelImage3DT<P>::restore( std::string filename )
  658. {
  659. // simple raw format
  660. FILE *f = fopen( filename.c_str(), "r" );
  661. if ( f == NULL ) {
  662. fprintf( stderr, "MultiChannelImage3DT::store: error reading from %s\n", filename.c_str() );
  663. exit( -1 );
  664. }
  665. fread( &xsize, sizeof( int ), 1, f );
  666. fread( &ysize, sizeof( int ), 1, f );
  667. fread( &zsize, sizeof( int ), 1, f );
  668. fread( &numChannels, sizeof( uint ), 1, f );
  669. if ( numChannels > 0 ) {
  670. reInit( xsize, ysize, zsize, numChannels );
  671. for ( uint channel = 0 ; channel < numChannels ; channel++ )
  672. {
  673. assert( data[channel] != NULL );
  674. fread( data[channel], sizeof( P ), xsize*ysize*zsize, f );
  675. }
  676. } else {
  677. freeData();
  678. data = NULL;
  679. }
  680. fclose( f );
  681. }
  682. template<class P>
  683. int MultiChannelImage3DT<P>::width() const
  684. {
  685. return xsize;
  686. }
  687. template<class P>
  688. int MultiChannelImage3DT<P>::height() const
  689. {
  690. return ysize;
  691. }
  692. template<class P>
  693. int MultiChannelImage3DT<P>::depth() const
  694. {
  695. return zsize;
  696. }
  697. template<class P>
  698. int MultiChannelImage3DT<P>::channels() const
  699. {
  700. return ( int )numChannels;
  701. }
  702. template<class P>
  703. int MultiChannelImage3DT<P>::getPixelInt( int x, int y, int z, int channel ) const
  704. {
  705. throw( "this type is not implemented\n" );
  706. return -1;
  707. }
  708. template<class P>
  709. double MultiChannelImage3DT<P>::getPixelFloat( int x, int y, int z, int channel ) const
  710. {
  711. throw( "this type is not implemented\n" );
  712. return -1.0;
  713. }
  714. template<class P>
  715. void MultiChannelImage3DT<P>::setPixelInt( int x, int y, int z, int channel, int pixel )
  716. {
  717. throw( "this type is not implemented\n" );
  718. }
  719. template<class P>
  720. void MultiChannelImage3DT<P>::setPixelFloat( int x, int y, int z, int channel, double pixel )
  721. {
  722. throw( "this type is not implemented\n" );
  723. }
  724. #define SET_FUNCS_PROTO_MACRO3D(MYTYPE) \
  725. template<>\
  726. int MultiChannelImage3DT<MYTYPE>::getPixelInt(int x, int y, int z, int channel) const;\
  727. template<>\
  728. double MultiChannelImage3DT<MYTYPE>::getPixelFloat(int x, int y, int z, int channel) const;\
  729. template<>\
  730. void MultiChannelImage3DT<MYTYPE>::setPixelInt(int x, int y, int z, int channel, int pixel);\
  731. template<>\
  732. void MultiChannelImage3DT<MYTYPE>::setPixelFloat(int x, int y, int z, int channel, double pixel);
  733. SET_FUNCS_PROTO_MACRO3D( double )
  734. SET_FUNCS_PROTO_MACRO3D( int )
  735. SET_FUNCS_PROTO_MACRO3D( long int )
  736. SET_FUNCS_PROTO_MACRO3D( float )
  737. SET_FUNCS_PROTO_MACRO3D( unsigned int )
  738. } // namespace