1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23
24 #ifndef __VOLUME_H__
25 #define __VOLUME_H__
26
27 #include "voxel.h"
28 #include <vcg/space/index/grid_static_ptr.h>
29
30 namespace vcg {
31
32 // forward definition
33 template < class VOL >
34 class VolumeIterator;
35
36 //******************************************
37 //******************************************
38 //typedef Voxel<float> Voxelf;
39
SFormat(const char * f,...)40 const char *SFormat( const char * f, ... )
41 {
42 static char buf[4096];
43 va_list marker;
44 va_start( marker, f );
45 vsprintf(buf,f,marker);
46 va_end( marker );
47 return buf;
48 }
49
50
51 template<class VOX_TYPE, class SCALAR_TYPE=float>
52 class Volume {
53 public:
54 typedef SCALAR_TYPE scalar;
55 typedef Point3<scalar> Point3x;
56 typedef Box3<scalar> Box3x;
57
58 typedef VOX_TYPE voxel_type;
59
BLOCKSIDE()60 static int BLOCKSIDE() { return 8;}
Volume()61 Volume(){
62 SetDefaultParam();
63 }
64 // I dati veri e propri
65 // Sono contenuti in un vettore di blocchi.
66 std::vector< std::vector<VOX_TYPE> > rv;
67 Box3x bbox;
68
69 _int64 AskedCells;
70 Point3x dim; /// Dimensione spaziale (lunghezza lati) del bbox
71
72 Point3i sz; /// Dimensioni griglia come numero di celle per lato
73
74 Point3i ssz; /// Dimensioni sottoblocco in esame come numero di celle per lato
75
76 Point3i rsz; /// Dimensioni macro griglia dei blocchi in cui e' suddiviso il volume (ogni blocco e' BLOCKSIDE()^3 celle)
77
78 Point3i asz; /// Dimensioni macro griglia dei blocchi relativa al sottoblocco in questione (quello allocato davvero!)
79
80
81 Point3x voxel; /// Dimensioni di una cella
82
83
84 int WN,WP; // di quanti vox mi devo allargare per settare la manhattan distance in neg e pos
85 int DeltaVoxelSafe; // di quanti vox mi devo allargare per stare sicuro nel fare solo una sottoparte.
86
ISize()87 const Point3i ISize() { return sz; }
88 private :
89 // offset e distanze varie precalcolate una volta per tutte
90 Point3f nnf[26];
91 Point3i nni[26];
92 float len[26];
93 float slen[26];
94
95 /// Gestione sottoparte
96 Point3i div;
97 Point3i pos;
98 public:
99 Box3i SubPart; // Sottoparte del volume da considerare ufficialmente
100 Box3x SubBox; // BBox della sottoparte del volume da considerare in coord assolute
101 Box3i SubPartSafe; // come sopra ma aumentati per sicurezza.
102 Box3x SubBoxSafe;
103
104 FILE *LogFP;
105 bool Verbose; // se true stampa un sacco di info in piu su logfp;
106
SetDefaultParam()107 void SetDefaultParam(){
108 WN=0;
109 WP=1;
110 //WN=-2;//
111 //WP=3;
112 DeltaVoxelSafe=BLOCKSIDE();
113 Verbose=true;
114 LogFP=stderr;
115 }
116
Init(const Volume & VV)117 void Init(const Volume &VV)
118 {
119 SetDefaultParam();
120 WN=VV.WN;
121 WP=VV.WP;
122 DeltaVoxelSafe=VV.DeltaVoxelSafe;
123 Init(VV.AskedCells,VV.bbox,VV.div,VV.pos);
124 }
125
126 void Init(_int64 cells, Box3x bb, Point3i _div=Point3i(1,1,1), Point3i _pos=Point3i(0,0,0))
127 {
128 Point3d voxdim;voxdim.Import(bb.max-bb.min);
129 AskedCells=cells;
130 vcg::BestDim<double>( cells, voxdim, sz );
131 bbox=bb;
132 /*
133 printf("grid of ~%i kcells: %d x %d x %d \n",int(cells/1000),sz[0],sz[1],sz[2]);
134 printf("grid voxel size of %f %f %f\n",voxdim[0]/sz[0],voxdim[1]/sz[1],voxdim[2]/sz[2]);
135 */
136 // il box deve essere multipli di BLOCKSIDE()
137 sz=((sz/BLOCKSIDE())+Point3i(1,1,1))*BLOCKSIDE();
138
139
140 rsz=sz/BLOCKSIDE();
141 if(sz!=rsz*BLOCKSIDE()) {
142 assert(0); // il box deve essere multipli di BLOCKSIDE()
143 exit(-1);
144 }
145
146 dim=bbox.max-bbox.min;
147 voxel[0]=dim[0]/sz[0];
148 voxel[1]=dim[1]/sz[1];
149 voxel[2]=dim[2]/sz[2];
150
151 SetSubPart(_div,_pos);
152 ssz=SubPartSafe.max-SubPartSafe.min;
153 asz=ssz/BLOCKSIDE() + Point3i(1,1,1);
154 rv.clear();
155 rv.resize(asz[0]*asz[1]*asz[2]);
156 for(size_t i=0;i<rv.size();++i)
157 rv[i].resize(0,VOX_TYPE::Zero());
158 SetDim(bb);
159 }
160
161 private:
162
163 // Da chiamare sempre DOPO la resize...
SetDim(const Box3x &)164 void SetDim(const Box3x & /*bb*/)
165 {
166
167 // Setup the precomputed offsets and offset normals
168 int cnt=0,x,y,z;
169 for(z=-1;z<=1;++z){
170 for(y=-1;y<=1;++y){
171 for(x=-1;x<=1;++x)
172 if(x!=0 || y!=0 || z!=0)
173 {
174 nnf[cnt]=Point3f(x,y,z);
175 len[cnt]=nnf[cnt].Norm();
176 slen[cnt]=nnf[cnt].SquaredNorm();
177 nnf[cnt].Normalize();
178 nni[cnt]=Point3i(x,y,z);
179 cnt++;
180 }
181 }
182 }
183 }
184
185 /*
186 Parametri
187 div indica il numero di blocchi da fare lungo i vari assi (sempre >=1)
188 pos indicano la coord del subbox da prendere in considerazione (sempre >=0 && < xdiv,ydiv,zdiv)
189 */
190
SetSubPart(Point3i _div,Point3i _pos)191 void SetSubPart(Point3i _div, Point3i _pos)
192 {
193 int k;
194 // Controllo correttezza parametri.
195 for(k=0;k<3;++k)
196 {
197 assert(_div[k]>0);
198 if(_div[k]==0){
199 printf("Error in subbox definition:\n the subdiv settings must be greater than 0; it was %i %i %i\n",_div[0],_div[1],_div[2]);
200 exit(-1);
201 }
202 if(_pos[k]<0 || _pos[k]>=_div[k]){
203 printf("Error in subbox definition:\n the Position of the subbox must be between (0,0,0) and (%i,%i,%i); it was %i %i %i\n",_div[0],_div[1],_div[2],_pos[0],_pos[1],_pos[2]);
204 exit(-1);
205 }
206 }
207
208 div=_div;
209 pos=_pos;
210
211 // Setting the subpart under analisys
212 for(k=0;k<3;++k)
213 {
214 SubPart.min[k]= pos[k]*sz[k]/div[k];
215 SubPart.max[k]=(pos[k]+1)*sz[k]/div[k];
216 SubBox.min[k]= bbox.min[k]+SubPart.min[k]*voxel[k];
217 SubBox.max[k]= bbox.min[k]+SubPart.max[k]*voxel[k];
218 }
219
220 // Setting the Safe Subpart under analisys
221 SubPartSafe=SubPart;
222 for(k=0;k<3;++k)
223 {
224 SubPartSafe.min[k] -= DeltaVoxelSafe;
225 SubPartSafe.max[k] += DeltaVoxelSafe;
226
227 if( SubPartSafe.min[k]< 0 ) SubPartSafe.min[k] = 0;
228 if( SubPartSafe.max[k]> sz[k] ) SubPartSafe.max[k] = sz[k];
229 SubBoxSafe.min[k]= bbox.min[k]+SubPartSafe.min[k]*voxel[k];
230 SubBoxSafe.max[k]= bbox.min[k]+SubPartSafe.max[k]*voxel[k];
231 }
232 /*
233 printf(" Computing only subvolume: (%d x %d x %d)= %dk cells \n"
234 " %d,%d,%d -> %d,%d,%d\n"
235 ,SubPart.DimX(),SubPart.DimY(),SubPart.DimZ(),(int)(((__int64)SubPart.DimX()*(__int64)SubPart.DimY()*(__int64)SubPart.DimZ())/1000)
236 ,SubPart.min[0] ,SubPart.min[1] ,SubPart.min[2]
237 ,SubPart.max[0] ,SubPart.max[1] ,SubPart.max[2] );
238 */
239 }
240
241 public:
242
243 // Sa
244 /*bool Write(string filename, const float &minv, const float &maxv )
245 {
246 FILE *fp;
247 if(div!=Point3i(1,1,1)) {
248 string subvoltag;
249 GetSubVolumeTag(subvoltag);
250 filename+=subvoltag;
251 }
252 string datname=filename;
253 string rawname=filename;
254 datname+=".dat";
255 rawname+=".raw";
256
257 fp=fopen(datname,"w");
258
259 fprintf(fp,"ObjectFileName: %s\n",rawname.c_str());
260 fprintf(fp,"TaggedFileName: ---\n");
261 fprintf(fp,"Resolution: %i %i %i\n",SubPart.max[0]-SubPart.min[0],SubPart.max[1]-SubPart.min[1],SubPart.max[2]-SubPart.min[2]);
262 fprintf(fp,"SliceThickness: %f %f %f\n",voxel[2]/voxel[0],voxel[2]/voxel[1],voxel[2]/voxel[2]);
263 fprintf(fp,"Format: UCHAR\n");
264 fprintf(fp,"NbrTags: 0\n");
265 fprintf(fp,"ObjectType: TEXTURE_VOLUME_OBJECT\n");
266 fprintf(fp,"ObjectModel: RGBA\n");
267 fprintf(fp,"GridType: EQUIDISTANT\n");
268
269 fclose(fp);
270 fp=fopen(rawname,"wb");
271 if(!fp)
272 {
273 printf("Error: unable ro open output volume file '%s'\n",filename);
274 return false;
275 }
276
277 int i,j,k;
278 for(k=SubPart.min[2];k<SubPart.max[2];++k)
279 for(j=SubPart.min[1];j<SubPart.max[1];++j)
280 for(i=SubPart.min[0];i<SubPart.max[0];++i)
281 {
282 float fv=V(i,j,k).V();
283 fv=(fv-minv)/(maxv-minv);
284 if(fv<0) fv=0;
285 else if(fv>1) fv=1;
286 fv=((fv*2.0f)-1.0f)*127;
287 char fs= (char) fv;
288 fwrite(&fs,sizeof(char),1,fp);
289 }
290 fclose(fp);
291 return true;
292 }*/
AbsPos(Point3i pi,Point3x & p)293 void AbsPos(Point3i pi, Point3x &p)
294 {
295 p[0]=bbox.min[0]+pi[0]*voxel[0];
296 p[1]=bbox.min[1]+pi[1]*voxel[1];
297 p[2]=bbox.min[2]+pi[2]*voxel[2];
298 }
299
300
GetSubVolumeTag(std::string & subtag)301 void GetSubVolumeTag(std::string &subtag)
302 {
303 char buf[32];
304 if (div[0]<= 10 && div[1]<= 10 && div[2]<= 10 ) sprintf(buf,"_%01d%01d%01d",pos[0],pos[1],pos[2]);
305 else if(div[0]<= 100 && div[1]<= 100 && div[2]<= 100 ) sprintf(buf,"_%02d%02d%02d",pos[0],pos[1],pos[2]);
306 else sprintf(buf,"_%03d%03d%03d",pos[0],pos[1],pos[2]);
307 subtag=buf;
308 }
309
310 /*
311 * Compute the offset <lpos> inside the subblock <rpos> of voxel (x,y,z).
312 * return true if the subblock is allocated.
313 */
Pos(const int & _x,const int & _y,const int & _z,int & rpos,int & lpos)314 bool Pos(const int &_x,const int &_y,const int &_z, int & rpos,int &lpos) const
315 {
316 int x=_x-SubPartSafe.min[0]; int y=_y-SubPartSafe.min[1]; int z=_z-SubPartSafe.min[2];
317
318 assert(_x>=SubPartSafe.min[0] && _x<SubPartSafe.max[0] &&
319 _y>=SubPartSafe.min[1] && _y<SubPartSafe.max[1] &&
320 _z>=SubPartSafe.min[2] && _z<SubPartSafe.max[2]);
321
322 // assert(x>=0 && x<sz[0] && y>=0 && y<sz[1] && z>=0 && z<sz[2]);
323
324 int rx=x/BLOCKSIDE(); int ry=y/BLOCKSIDE(); int rz=z/BLOCKSIDE();
325 assert(rx>=0 && rx<asz[0] && ry>=0 && ry<asz[1] && rz>=0 && rz<asz[2]);
326 rpos = rz*asz[0]*asz[1]+ry*asz[0]+rx;
327 assert(rpos < int(rv.size()));
328 int lx = x%BLOCKSIDE(); int ly = y%BLOCKSIDE(); int lz = z % BLOCKSIDE();
329 lpos = lz*BLOCKSIDE()*BLOCKSIDE()+ly*BLOCKSIDE()+lx;
330 if(rv[rpos].empty()) return false;
331 return true;
332 }
333
334 /*
335 Funzione inversa della precedente
336 Date due posizioni rpos e lpos restituisce x,y,z assoluti
337 */
IPos(int & x,int & y,int & z,const int & rpos,const int & lpos)338 bool IPos(int &x,int &y,int &z, const int & rpos, const int &lpos) const
339 {
340 assert (rpos>=0 && lpos >=0);
341
342 int rz = rpos / (asz[0]*asz[1]); int remainder = rpos % (asz[0]*asz[1]);
343 int ry = ( remainder ) / asz[0] ;
344 int rx = remainder % asz[0];
345
346 assert(rx>=0 && rx<asz[0] && ry>=0 && ry<asz[1] && rz>=0 && rz<asz[2]);
347
348 int lz = lpos / (BLOCKSIDE()*BLOCKSIDE()); int lemaindel = lpos % (BLOCKSIDE()*BLOCKSIDE());
349 int ly = ( lemaindel ) / BLOCKSIDE();
350 int lx = lemaindel % BLOCKSIDE();
351
352 x = rx*BLOCKSIDE()+lx;
353 y = ry*BLOCKSIDE()+ly;
354 z = rz*BLOCKSIDE()+lz;
355
356 x+=SubPartSafe.min[0];
357 y+=SubPartSafe.min[1];
358 z+=SubPartSafe.min[2];
359
360 assert(x>=0 && x<sz[0] && y>=0 && y<sz[1] && z>=0 && z<sz[2]);
361 //int trpos,tlpos;
362 //assert(rv[rpos].size()>0);
363 //Pos(x,y,z,trpos,tlpos);
364 //assert(trpos==rpos && tlpos == lpos);
365 return true;
366 }
367
Alloc(int rpos,const VOX_TYPE & zeroval)368 void Alloc(int rpos, const VOX_TYPE &zeroval)
369 {
370 rv[rpos].resize(BLOCKSIDE()*BLOCKSIDE()*BLOCKSIDE(),zeroval);
371 }
372 /************************************/
373 // Funzioni di accesso ai dati
ValidCell(const Point3i & p1,const Point3i & p2)374 bool ValidCell(const Point3i &p1, const Point3i &p2) const
375 {
376 if(!cV(p1.X(),p1.Y(),p1.Z()).B() ) return false;
377 if(!cV(p2.X(),p1.Y(),p1.Z()).B() ) return false;
378 if(!cV(p1.X(),p2.Y(),p1.Z()).B() ) return false;
379 if(!cV(p2.X(),p2.Y(),p1.Z()).B() ) return false;
380 if(!cV(p1.X(),p1.Y(),p2.Z()).B() ) return false;
381 if(!cV(p2.X(),p1.Y(),p2.Z()).B() ) return false;
382 if(!cV(p1.X(),p2.Y(),p2.Z()).B() ) return false;
383 if(!cV(p2.X(),p2.Y(),p2.Z()).B() ) return false;
384 return true;
385 }
386
Val(const int & x,const int & y,const int & z)387 float Val(const int &x,const int &y,const int &z) const {
388 if(!cV(x,y,z).B()) return 1000;
389 return cV(x,y,z).V();
390 //else return numeric_limits<float>::quiet_NaN( );
391 }
392
V(const int & x,const int & y,const int & z)393 VOX_TYPE &V(const int &x,const int &y,const int &z) {
394 int rpos,lpos;
395 if(!Pos(x,y,z,rpos,lpos)) Alloc(rpos,VOX_TYPE::Zero());
396 return rv[rpos][lpos];
397 }
398
cV(const int & x,const int & y,const int & z)399 const VOX_TYPE &cV(const int &x,const int &y,const int &z) const
400 {
401 int rpos,lpos;
402 if(!Pos(x,y,z,rpos,lpos)) return VOX_TYPE::Zero();
403 else return rv[rpos][lpos];
404 }
V(const int & x,const int & y,const int & z)405 const VOX_TYPE &V(const int &x,const int &y,const int &z) const
406 {
407 int rpos,lpos;
408 if(!Pos(x,y,z,rpos,lpos)) return VOX_TYPE::Zero();
409 else return rv[rpos][lpos];
410 }
411 /************************************/
Fill(VOX_TYPE (__cdecl * func)(const Point3i & p))412 void Fill(VOX_TYPE (__cdecl *func)(const Point3i &p) )
413 {
414 int x,y,z;
415 for(z=0;z<sz[2];++z)
416 for(y=0;y<sz[1];++y)
417 for(x=0;x<sz[0];++x)
418 {
419 Point3i p(x,y,z);
420 V(x,y,z)=func(p);
421 }
422 }
423
Fill(VOX_TYPE const p)424 void Fill(VOX_TYPE const p)
425 {
426 int x,y,z;
427 for(z=0;z<sz[2];++z)
428 for(y=0;y<sz[1];++y)
429 for(x=0;x<sz[0];++x)
430 {
431 V(x,y,z)=p;
432 }
433 }
434
435
436 // Copia sul volume corrente una versione smoothed del volume in ingresso S.
437 // il parametro serve a specificare il range di valori di campo vicini allo zero che non vanno mediati!
438 // questo perche se si smootha anche sullo zero si smoota anche dove e' bene allineato
439
440 void CopySmooth( Volume<VOX_TYPE> &S, scalar SafeZone=1, scalar SafeQuality=0)
441 {
442 if(sz!=S.sz)
443 {
444 printf("Error");
445 exit(-1);
446 }
447 int lcnt=0;
448 VolumeIterator< Volume > vi(S);
449 vi.Restart();
450 vi.FirstNotEmpty();
451 // const Voxelf *VC;
452 while(vi.IsValid())
453 // scandisci il volume in ingresso, per ogni voxel non vuoto del volume
454 // in ingresso calcola la media con gli adiacenti
455 {
456 if((*vi).B())
457 {
458 int x,y,z;
459 IPos(x,y,z,vi.rpos,vi.lpos);
460 if(Bound1(x,y,z))
461 {
462 VOX_TYPE &VC = V(x,y,z);
463 for(int i=0;i<26;++i)
464 {
465 VOX_TYPE &VV= S.V(x+nni[i][0],y+nni[i][1],z+nni[i][2]);
466 if(VV.B()) VC+=VV;
467 }
468 lcnt++;
469
470 /*
471 Voxelf &VV=V(x,y,z);
472 //Voxelf &VV=rv[vi.rpos][vi.lpos];
473 for(int i=0;i<26;++i)
474 {
475 VC = &(S.V(x+nni[i][0],y+nni[i][1],z+nni[i][2]));
476 if(VC->b)
477 {
478 VV+=*VC;
479 lcnt++;
480 }
481 }*/
482 }
483 }
484 vi.Next();
485 if(vi.IsValid()) vi.FirstNotEmpty();
486 //if((lcnt%100)==0) vi.Dump();
487
488 }
489 // Step 2,
490 // dopo aver calcolato la media,
491
492 VolumeIterator< Volume > svi(*this);
493 svi.Restart();
494 svi.FirstNotEmpty();
495 int smoothcnt=0;
496 int preservedcnt=0;
497 int blendedcnt=0;
498 const float FieldBorder = 1; // dove finisce la transizione tra la zona safe e quella smoothed
499 const float EndFBorderZone = SafeZone+FieldBorder;
500 const float EndQBorderZone = SafeQuality*1.5;
501 const float QBorder = EndQBorderZone-SafeQuality; // dove finisce la transizione tra la zona safe e quella smoothed
502 while(svi.IsValid())
503 {
504 if((*svi).Cnt()>0)
505 {
506 VOX_TYPE &sv=S.rv[svi.rpos][svi.lpos];
507 (*svi).Normalize(1); // contiene il valore mediato
508 float SafeThr = fabs(sv.V());
509
510 // Se la qualita' e' bassa o se siamo distanti si smootha sempre
511 // se siamo vicini allo zero e con buona qualita' si deve fare attenzione
512 if(SafeThr<EndFBorderZone && sv.Q() > EndQBorderZone)
513 { // se il voxel corrente aveva un valore < safezone E qualita' > SafeQuality
514 // allora si copia il valore originale di S
515 if((SafeThr <= SafeZone) && sv.Q() > SafeQuality )
516 {
517 (*svi)=sv;
518 (*svi).SetB(true);
519 ++preservedcnt;
520 }
521 else
522 { // Siamo nella zona di transizione o per field o per quality
523 float blendq= std::max(0.0f,std::min(1.0f,(EndQBorderZone-sv.Q())/QBorder));
524 float blendf= std::max(0.0f,std::min(1.0f,(EndFBorderZone-SafeThr)/FieldBorder));
525 float BlendFactor = 1.0-std::max(blendf,blendq); // quanto del voxel originale <sv> si prende;
526 (*svi).Blend(sv,BlendFactor);
527 ++blendedcnt;
528 }
529 }
530 ++smoothcnt;
531 }
532 svi.Next();
533 if(svi.IsValid()) svi.FirstNotEmpty();
534 }
535
536 if(Verbose) fprintf(LogFP,"CopySmooth %i voxels, %i preserved, %i blended\n",smoothcnt,preservedcnt,blendedcnt);
537 }
538
Merge(Volume<VOX_TYPE> & S)539 void Merge(Volume<VOX_TYPE> &S)
540 {
541 VolumeIterator< Volume > svi(S);
542 svi.Restart();
543 svi.FirstNotEmpty();
544 int loccnt=0;
545
546 while(svi.IsValid())
547 {
548 if((*svi).B())
549 {
550 int x,y,z;
551 IPos(x,y,z,svi.rpos,svi.lpos);
552 if(cV(x,y,z).B()) V(x,y,z).Merge( (*svi));
553 else {
554 V(x,y,z).Set((*svi));
555 V(x,y,z).SetB(true);
556 }
557 ++loccnt;
558 }
559 svi.Next();
560 if(svi.IsValid()) svi.FirstNotEmpty();
561 }
562
563 printf("Merge2 %i voxels\n",loccnt);
564
565 }
566
567
Interize(Point3x & vert)568 void Interize( Point3x & vert ) const // OK
569 {
570 for(int j=0;j<3;++j)
571 {
572 assert(vert[j]>=bbox.min[j]);
573 assert(vert[j]<=bbox.max[j]);
574 vert[j] = (vert[j] - bbox.min[j]) * sz[j] / (bbox.max[j] - bbox.min[j]);
575 }
576 }
577
578
DeInterize(Point3x & vert)579 void DeInterize( Point3x & vert ) const // OK
580 {
581 for(int j=0;j<3;++j)
582 vert[j] = vert[j] * (bbox.max[j] - bbox.min[j]) / sz[j] + bbox.min[j];
583 }
584
SplatVert(const Point3x & v0,double quality,const Point3x & nn,Color4b c)585 bool SplatVert( const Point3x & v0, double quality, const Point3x & nn, Color4b c)
586 {
587 Box3i ibox;
588
589 assert(math::Abs(SquaredNorm(nn) - 1.0) < 0.0001); // Just a safety check that the vertex normals are NORMALIZED!
590 ibox.min=Point3i(floor(v0[0]),floor(v0[1]),floor(v0[2]));
591 ibox.max=Point3i( ceil(v0[0]), ceil(v0[1]), ceil(v0[2]));
592 ibox.Intersect(SubPartSafe);
593
594 ibox.max[0] = std::min(SubPartSafe.max[0]-1,ibox.max[0]);
595 ibox.max[1] = std::min(SubPartSafe.max[1]-1,ibox.max[1]);
596 ibox.max[2] = std::min(SubPartSafe.max[2]-1,ibox.max[2]);
597
598
599 // Skip faces not colliding current subvolume.
600 if(ibox.IsEmpty())
601 {
602 // point outside the box do nothing
603 return false;
604 }
605
606 Point3x iV, deltaIV;
607
608 // Now scan the eight voxel surrounding the splat
609 // and fill them with the distance from the plane
610 for(iV[0]=ibox.min[0]; iV[0]<=ibox.max[0]; ++iV[0])
611 for(iV[1]=ibox.min[1]; iV[1]<=ibox.max[1]; ++iV[1])
612 for(iV[2]=ibox.min[2]; iV[2]<=ibox.max[2]; ++iV[2])
613 {
614 deltaIV = v0-iV;
615 scalar offset = deltaIV * nn;
616 VOX_TYPE &VV=V(iV[0],iV[1],iV[2]);
617 VV=VOX_TYPE(offset,nn,quality,c);
618 }
619 return true;
620 }
621
622 template <const int CoordZ>
RasterFace(const int sx,const int ex,const int sy,const int ey,scalar dist,const Point3x & norm,scalar quality,const Point3x & v0,const Point3x & v1,const Point3x & v2,const Point3x & d10,const Point3x & d21,const Point3x & d02)623 void RasterFace(const int sx, const int ex, const int sy, const int ey,
624 scalar dist, const Point3x &norm, scalar quality,
625 const Point3x &v0, const Point3x &v1, const Point3x &v2,
626 const Point3x &d10, const Point3x &d21, const Point3x &d02)
627 {
628 const scalar EPS = scalar(1e-12);
629 const int crd0 = CoordZ;
630 const int crd1 = (CoordZ+1)%3;
631 const int crd2 = (CoordZ+2)%3;
632 assert(fabs(norm[crd0])+0.001 > fabs(norm[crd1]));
633 assert(fabs(norm[crd0])+0.001 > fabs(norm[crd2]));
634 scalar x,y;
635 for(x=sx;x<=ex;++x)
636 for(y=sy;y<=ey;++y)
637 {
638 scalar n0 = (x-v0[crd1])*d10[crd2] - (y-v0[crd2])*d10[crd1];
639 scalar n1 = (x-v1[crd1])*d21[crd2] - (y-v1[crd2])*d21[crd1];
640 scalar n2 = (x-v2[crd1])*d02[crd2] - (y-v2[crd2])*d02[crd1];
641
642 if( (n0>-EPS && n1>-EPS && n2>-EPS) ||
643 (n0< EPS && n1< EPS && n2< EPS ) )
644 {
645 scalar iz = ( dist - x*norm[crd1] - y*norm[crd2] ) / norm[crd0];
646 //assert(iz>=fbox.min[2] && iz<=fbox.max[2]);
647 AddIntercept<CoordZ>(x,y,iz, quality, norm );
648 }
649 }
650 }
651
652 // Si sa che la faccia ha una intercetta sull'asse z-dir di coord xy alla posizione z;
653 // quindi si setta nei 2 vertici prima e 2 dopo la distanza corrispondente.
654 template<int CoordZ>
AddIntercept(const int x,const int y,const scalar z,const scalar q,const Point3f & n)655 void AddIntercept( const int x, const int y, const scalar z, const scalar q, const Point3f &n )
656 {
657 scalar esgn = (n[CoordZ] > 0 ? -1 : 1);
658 int zint = floor(z);
659 scalar dist=z-zint; // sempre positivo e compreso tra zero e uno
660
661 for(int k=WN;k<=WP;k++)
662 {
663 if(zint+k >= SubPartSafe.min[CoordZ] && zint+k < SubPartSafe.max[CoordZ])
664 {
665 VOX_TYPE *VV;
666 if(CoordZ==2) VV=&V(x,y,zint+k);
667 if(CoordZ==1) VV=&V(y,zint+k,x);
668 if(CoordZ==0) VV=&V(zint+k,x,y);
669 scalar nvv= esgn*( k-dist);
670 if(!VV->B() || fabs(VV->V()) > fabs(nvv)) {
671 *VV=VOX_TYPE(nvv,n,q);
672 }
673 }
674 }
675 }
676
677 // assume che i punti della faccia in ingresso siano stati interized
ScanFace2(const Point3x & v0,const Point3x & v1,const Point3x & v2,scalar quality,const Point3x & norm)678 bool ScanFace2( const Point3x & v0, const Point3x & v1, const Point3x & v2,
679 scalar quality, const Point3x & norm)//, const int name ) // OK
680 {
681
682 Box3x fbox; // Bounding Box della faccia (double)
683 Box3i ibox; // Bounding Box della faccia (int)
684 int sx,sy,sz; // Bounding Box intero
685 int ex,ey,ez;
686
687 // Calcolo bbox della faccia
688 fbox.Set(v0);
689 fbox.Add(v1);
690 fbox.Add(v2);
691
692 // BBox intero (nota che il cast a int fa truncation (funge solo perche' v0,v1,v2 sono positivi)
693
694 ibox.min[0] =sx = floor(fbox.min[0]); if( ((scalar)sx)!=fbox.min[0] ) ++sx; // necessario se il punto e'approx a .9999
695 ibox.min[1] =sy = floor(fbox.min[1]); if( ((scalar)sy)!=fbox.min[1] ) ++sy;
696 ibox.min[2] =sz = floor(fbox.min[2]); if( ((scalar)sz)!=fbox.min[2] ) ++sz;
697 ibox.max[0] =ex = floor(fbox.max[0]);
698 ibox.max[1] =ey = floor(fbox.max[1]);
699 ibox.max[2] =ez = floor(fbox.max[2]);
700 // Skip faces not colliding current subvolume.
701 if(!ibox.Collide(SubPartSafe)) return false;
702
703 Point3x d10 = v1 - v0;
704 Point3x d21 = v2 - v1;
705 Point3x d02 = v0 - v2;
706
707 assert(norm.Norm() >0.999f && norm.Norm()<1.001f);
708 // assert(0);
709 scalar dist = norm * v0;
710
711
712 /**** Rasterizzazione bbox ****/
713
714 // Clamping dei valori di rasterizzazione al subbox corrente
715 sx = std::max(SubPartSafe.min[0],sx); ex = std::min(SubPartSafe.max[0]-1,ex);
716 sy = std::max(SubPartSafe.min[1],sy); ey = std::min(SubPartSafe.max[1]-1,ey);
717 sz = std::max(SubPartSafe.min[2],sz); ez = std::min(SubPartSafe.max[2]-1,ez);
718
719 if(fabs(norm[0]) > fabs(norm[1]) && fabs(norm[0])>fabs(norm[2])) RasterFace<0>(sy,ey,sz,ez,dist,norm,quality,v0,v1,v2,d10,d21,d02);
720 else if(fabs(norm[1]) > fabs(norm[0]) && fabs(norm[1])>fabs(norm[2])) RasterFace<1>(sz,ez,sx,ex,dist,norm,quality,v0,v1,v2,d10,d21,d02);
721 else RasterFace<2>(sx,ex,sy,ey,dist,norm,quality,v0,v1,v2,d10,d21,d02);
722
723
724 return true;
725 }
726
727 // assume che i punti della faccia in ingresso siano stati interized
ScanFace(const Point3x & v0,const Point3x & v1,const Point3x & v2,double quality,const Point3x & nn)728 bool ScanFace( const Point3x & v0, const Point3x & v1, const Point3x & v2,
729 double quality, const Point3x & nn)//, const int name ) // OK
730 {
731 const scalar EPS = scalar(1e-12);
732 // const scalar EPS_INT = scalar(1e-20);
733 const scalar EPS_INT = .3f;//scalar(1e-20);
734
735 #ifndef _NDEBUG
736 if(quality==0.0)
737 {
738 printf("Zero quality face are not allowed!\n");
739 exit(-1);
740 }
741 #endif
742
743 // ++nfaces;
744
745 Box3x fbox; // Bounding Box della faccia (double)
746 Box3i ibox; // Bounding Box della faccia (int)
747 int sx,sy,sz; // Bounding Box intero
748 int ex,ey,ez;
749
750 // Calcolo bbox della faccia
751 fbox.Set(v0);
752 fbox.Add(v1);
753 fbox.Add(v2);
754
755 // BBox intero (nota che il cast a int fa truncation (funge solo perche' v0,v1,v2 sono positivi)
756
757 ibox.min[0] =sx = floor(fbox.min[0]); if( ((scalar)sx)!=fbox.min[0] ) ++sx; // necessario se il punto e'approx a .9999
758 ibox.min[1] =sy = floor(fbox.min[1]); if( ((scalar)sy)!=fbox.min[1] ) ++sy;
759 ibox.min[2] =sz = floor(fbox.min[2]); if( ((scalar)sz)!=fbox.min[2] ) ++sz;
760 ibox.max[0] =ex = floor(fbox.max[0]);
761 ibox.max[1] =ey = floor(fbox.max[1]);
762 ibox.max[2] =ez = floor(fbox.max[2]);
763 // Skip faces not colliding current subvolume.
764 if(!ibox.Collide(SubPartSafe)) return false;
765
766 /**** Dati per controllo intersezione ****/
767
768 // Versori degli spigoli della faccia
769
770 Point3x d10 = v1 - v0;
771 Point3x d21 = v2 - v1;
772 Point3x d02 = v0 - v2;
773
774 // Normale al piano della faccia e distanza origine
775
776 Point3x norm = d10 ^ d21;
777 norm.Normalize();
778 double dist = norm * v0;
779
780 /**** Rasterizzazione bbox ****/
781
782 int x,y,z;
783
784
785 // Clamping dei valori di rasterizzazione al subbox corrente
786 sx = std::max(SubPartSafe.min[0],sx); ex = std::min(SubPartSafe.max[0]-1,ex);
787 sy = std::max(SubPartSafe.min[1],sy); ey = std::min(SubPartSafe.max[1]-1,ey);
788 sz = std::max(SubPartSafe.min[2],sz); ez = std::min(SubPartSafe.max[2]-1,ez);
789
790 // Rasterizzazione xy
791
792 if(fabs(norm[2])>EPS_INT)
793 for(x=sx;x<=ex;++x)
794 for(y=sy;y<=ey;++y)
795 {
796 double n0 = ((double)x-v0[0])*d10[1] - ((double)y-v0[1])*d10[0];
797 double n1 = ((double)x-v1[0])*d21[1] - ((double)y-v1[1])*d21[0];
798 double n2 = ((double)x-v2[0])*d02[1] - ((double)y-v2[1])*d02[0];
799
800 if( (n0>-EPS && n1>-EPS && n2>-EPS) ||
801 (n0< EPS && n1< EPS && n2< EPS ))
802 {
803 double iz = ( dist - double(x)*norm[0] - double(y)*norm[1] ) / norm[2];
804 //assert(iz>=fbox.min[2] && iz<=fbox.max[2]);
805 AddXYInt(x,y,iz,-norm[2], quality, nn );
806 }
807 }
808
809 // Rasterizzazione xz
810
811 if(fabs(norm[1])>EPS_INT)
812 for(x=sx;x<=ex;++x)
813 for(z=sz;z<=ez;++z)
814 {
815 double n0 = ((double)x-v0[0])*d10[2] - ((double)z-v0[2])*d10[0];
816 double n1 = ((double)x-v1[0])*d21[2] - ((double)z-v1[2])*d21[0];
817 double n2 = ((double)x-v2[0])*d02[2] - ((double)z-v2[2])*d02[0];
818
819 if( (n0>-EPS && n1>-EPS && n2>-EPS) ||
820 (n0< EPS && n1< EPS && n2< EPS ))
821 {
822 double iy = ( dist - double(x)*norm[0] - double(z)*norm[2] ) / norm[1];
823 //assert(iy>=fbox.min[1] && iy<=fbox.max[1]);
824 AddXZInt(x,z,iy,-norm[1], quality,nn );
825 }
826 }
827
828 // Rasterizzazione yz
829
830 if(fabs(norm[0])>EPS_INT)
831 for(y=sy;y<=ey;++y)
832 for(z=sz;z<=ez;++z)
833 {
834 double n0 = ((double)y-v0[1])*d10[2] - ((double)z-v0[2])*d10[1];
835 double n1 = ((double)y-v1[1])*d21[2] - ((double)z-v1[2])*d21[1];
836 double n2 = ((double)y-v2[1])*d02[2] - ((double)z-v2[2])*d02[1];
837
838 if( (n0>-EPS && n1>-EPS && n2>-EPS) ||
839 (n0< EPS && n1< EPS && n2< EPS ) )
840 {
841 double ix = ( dist - double(y)*norm[1] - double(z)*norm[2] ) / norm[0];
842 //assert(ix>=fbox.min[0] && ix<=fbox.max[0]);
843 AddYZInt(y,z,ix,-norm[0], quality, nn );
844 }
845 }
846 return true;
847 }
848 // Si sa che la faccia ha una intercetta sull'asse z-dir di coord xy alla posizione z;
849 // quindi si setta nei 2 vertici prima e 2 dopo la distanza corrispondente.
850
AddXYInt(const int x,const int y,const double z,const double sgn,const double q,const Point3f & n)851 void AddXYInt( const int x, const int y, const double z, const double sgn, const double q, const Point3f &n )
852 { double esgn = (sgn<0 ? -1 : 1);//*max(fabs(sgn),0.001);
853 double dist=z-floor(z); // sempre positivo e compreso tra zero e uno
854 int zint = floor(z);
855 for(int k=WN;k<=WP;k++)
856 if(zint+k >= SubPartSafe.min[2] && zint+k < SubPartSafe.max[2])
857 {
858 VOX_TYPE &VV=V(x,y,zint+k);
859 double nvv= esgn*( k-dist);
860 if(!VV.B() || fabs(VV.V()) > fabs(nvv)) {
861 VV=VOX_TYPE(nvv,n,q);
862 }
863 }
864 }
AddYZInt(const int y,const int z,const double x,const double sgn,const double q,const Point3f & n)865 void AddYZInt( const int y, const int z, const double x, const double sgn, const double q, const Point3f &n )
866 { double esgn = (sgn<0 ? -1 : 1);//*max(fabs(sgn),0.001);
867 double dist=x-floor(x); // sempre positivo e compreso tra zero e uno
868 int xint = int(floor(x));
869 for(int k=WN;k<=WP;k++)
870 if(xint+k >= SubPartSafe.min[0] && xint+k < SubPartSafe.max[0])
871 {
872 VOX_TYPE &VV=V(xint+k,y,z);
873 double nvv= esgn*( k-dist);
874 if(!VV.B() || fabs(VV.V()) > fabs(nvv)) {
875 VV=VOX_TYPE(nvv,n,q);
876 }
877 }
878 }
AddXZInt(const int x,const int z,const double y,const double sgn,const double q,const Point3f & n)879 void AddXZInt( const int x, const int z, const double y, const double sgn, const double q, const Point3f &n )
880 { double esgn = (sgn<0 ? -1 : 1);//*max(fabs(sgn),0.001);
881 double dist=y-scalar(floor(y)); // sempre positivo e compreso tra zero e uno
882 int yint = floor(y);
883 for(int k=WN;k<=WP;k++)
884 if(yint+k >= SubPartSafe.min[1] && yint+k < SubPartSafe.max[1])
885 {
886 VOX_TYPE &VV=V(x,yint+k,z);
887 double nvv= esgn*( k-dist);
888 if(!VV.B() || fabs(VV.V()) > fabs(nvv)) {
889 VV=VOX_TYPE(nvv,n,q);
890 }
891 }
892 }
893
894
895
Dump(FILE * fp)896 void Dump(FILE *fp)
897 {
898 fprintf(fp,"Volume Info:\n");
899 fprintf(fp," BBbox %7.4f %7.4f %7.4f - %7.4f %7.4f %7.4f:\n",bbox.min[0],bbox.min[1],bbox.min[2],bbox.max[0],bbox.max[1],bbox.max[2]);
900 fprintf(fp," Size in voxels %7i %7i %7i (tot: %7.3f M):\n",sz[0],sz[1],sz[2],(double(sz[0]*sz[1])/1000000.0)*sz[2]);
901 fprintf(fp," Voxel dimension %7.4f %7.4f %7.4f \n",voxel[0],voxel[1],voxel[2]);
902
903 fprintf(fp," Size in MacroCell %7i %7i %7i (tot: %7.3f M):\n",rsz[0],rsz[1],rsz[2],double(rsz[0]*rsz[1]*rsz[2])/1000000.0);
904 fprintf(fp," Memory Info: \n Voxel Size %8li b Virtually needed mem %8i Mb\n",
905 sizeof(VOX_TYPE),int(sizeof(VOX_TYPE)*(_int64)(sz[0])*(_int64)(sz[1])*(_int64)(sz[2])/(1024*1024)));
906 if(div!=Point3i(1,1,1))
907 {
908 fprintf(fp," Subdivided in %6i %6i %6i (tot: %12i block):\n",div[0],div[1],div[2],div[0]*div[1]*div[2]);
909 fprintf(fp," Computing subblock %6i %6i %6i :\n",pos[0],pos[1],pos[2]);
910 fprintf(fp," %6i %6i %6i - %6i %6i %6i :\n",SubPart.min[0],SubPart.min[1],SubPart.min[2],SubPart.max[0],SubPart.max[1],SubPart.max[2]);
911 fprintf(fp," Safe %6i %6i %6i - %6i %6i %6i :\n",SubPartSafe.min[0],SubPartSafe.min[1],SubPartSafe.min[2],SubPartSafe.max[0],SubPartSafe.max[1],SubPartSafe.max[2]);
912
913 }
914 fprintf(fp,"\n");
915 }
916
Allocated()917 int Allocated()
918 {int cnt=0;
919 for(size_t i=0;i<rv.size();++i)
920 if(!rv[i].empty()) cnt++;
921 return cnt;
922 }
923
Bound1(const int x,const int y,const int z)924 bool Bound1(const int x, const int y, const int z)
925 {
926 return (x>SubPartSafe.min[0] && x < SubPartSafe.max[0]-1 ) &&
927 (y>SubPartSafe.min[1] && y < SubPartSafe.max[1]-1 ) &&
928 (z>SubPartSafe.min[2] && z < SubPartSafe.max[2]-1 ) ;
929 }
930
931 /*
932 Note sull'algoritmo di espansione:
933
934 Si riempie i voxel vuoti
935
936 Se il volume e' inizialmente riempito con ii valori delle intercette alla superficie
937 nei 2 vertici immediatamente adiacenti all'edge intersecato dalla superficie
938 si deve espadnere tale coampo in maniera sensata.
939
940 Notare che e' importante che non tutto il campo sia riempito con un approx ella distanza di hausdorf:
941 il campo deve essere "tagliato" sui bordi della superficie per evitare pasticci. Levoy riempie il campo
942 solo lungo la direzione dello scanner, io invece riempio lungo la normale alla superficie. In questo modo
943 si evita il problema che l'espansione e' legata all'acquisizione iniziale
944
945 */
946
947
Expand(scalar AngleThrRad)948 void Expand(scalar AngleThrRad)
949 {
950 int i;
951 VolumeIterator< Volume > vi(*this);
952
953 float CosThr=math::Cos(AngleThrRad);
954 // printf("Expand2 angle %f, %f\n",AngleThrRad,CosThr);
955 int loccnt=0;
956
957 vi.Restart();
958 vi.FirstNotEmpty();
959 while(vi.IsValid())
960 {
961 if((*vi).B()) // si espande solo i voxel con valori "validi"
962 {
963 int x,y,z;
964 IPos(x,y,z,vi.rpos,vi.lpos);
965 Point3f n=(*vi).N();
966 VOX_TYPE vtmp = (*vi);
967 if(Bound1(x,y,z))
968 for(i=0;i<26;++i)
969 {
970 float angle = -(nnf[i]*n); // cos angolo tra la normale alla superficie e la direzione di espansione
971 if( fabs(angle)> CosThr )
972 {
973 //bbfloat tt=(*vi).V();
974 vtmp.SetV((*vi).V()+len[i]*angle); // la nuova distanza e' la distanza rispetto al piano passante per il punto a cui si riferisce VV;
975 VOX_TYPE &VV= V(x+nni[i][0],y+nni[i][1],z+nni[i][2]);
976 if(!VV.B()){
977 VV+=vtmp;
978 loccnt++;
979 }
980 }
981 }
982 }
983 vi.Next();
984 if(vi.IsValid()) vi.FirstNotEmpty();
985 }
986 printf("Expand %8i ",loccnt);
987 Normalize(1);
988 }
989
990 // filla i buchi vuoti di un volume;
991 // scorre tutti i voxel pieni e aggiunge agli adiacenti vuoti tale valore;
992 // si riscorre tutti i voxel vuoti e se si sono sommati almeno thr valori sitiene.
993 // in pratica maggiore il valore di thr meno buchi si riempiono.
994
995 void Refill(const int thr,float maxdistance = std::numeric_limits<float>::max() )
996 {
997 int lcnt=0;
998 VolumeIterator< Volume > vi(*this);
999 vi.Restart();
1000 vi.FirstNotEmpty();
1001
1002 while(vi.IsValid())
1003 {
1004 if((*vi).B())
1005 {
1006 int x,y,z;
1007 IPos(x,y,z,vi.rpos,vi.lpos);
1008 if(Bound1(x,y,z))
1009 {
1010 for(int i=0;i<26;++i)
1011 {
1012 VOX_TYPE &VC= V(x+nni[i][0],y+nni[i][1],z+nni[i][2]);
1013 if(!VC.B()){
1014 if(VC.Cnt()==0) lcnt++;
1015 VC+=(*vi);
1016 }
1017 }
1018 }
1019 }
1020
1021 vi.Next();
1022
1023 if(vi.IsValid()) vi.FirstNotEmpty();
1024
1025
1026 }
1027 printf("ReFill %8i ",lcnt);
1028 Normalize(thr,maxdistance);
1029 }
1030
1031 /*
1032 Attraversa il volume e modifica il valore di campo in modo da creare una offset surface che si connetta
1033 bene con la superficie originale
1034 il parametro specifica dove dovrebbe passare l'offset surface
1035
1036 L'idea e' quella di fare un altro zero del distance field al threshold specificato
1037
1038 La cosa deve essere smooth
1039 quindi scelgo una funzione che abbia 2 zeri (in zero e in thr) e
1040 */
Offset(const float thr)1041 void Offset(const float thr)
1042 {
1043 int lcnt=0;
1044 VolumeIterator< Volume > vi(*this);
1045 vi.Restart();
1046 vi.FirstNotEmpty();
1047 float thr2=thr/2.0;
1048 while(vi.IsValid())
1049 {
1050 if((*vi).B())
1051 {
1052 float vv=(*vi).V();
1053 if(thr<0) if(vv<thr2) vv=thr-vv;
1054 if(thr>0) if(vv>thr2) vv=thr-vv;
1055
1056 (*vi).SetV(vv);
1057 }
1058
1059 vi.Next();
1060
1061 if(vi.IsValid()) vi.FirstNotEmpty();
1062
1063
1064 }
1065 printf("ReFill %8i ",lcnt);
1066 Normalize(thr);
1067 }
1068
1069 // prende un volume e mette a true il campo b di tutti i voxel che hanno un valore significativo.
1070 // thr indica il numero minimo di valori che si devono essere sommati sul voxel;
1071 // di solito e' uguale a 1;
1072 int Normalize(int thr, float maxdistance=std::numeric_limits<float>::max() )
1073 {
1074 VolumeIterator< Volume > vi(*this);
1075 vi.Restart();
1076 vi.FirstNotEmpty();
1077 int loccnt=0;
1078 while(vi.IsValid())
1079 {
1080 if(!(*vi).B())
1081 {
1082 if((*vi).Normalize(thr))
1083 ++loccnt;
1084 if(math::Abs((*vi).V())>maxdistance) *vi=VOX_TYPE::Zero();
1085 }
1086 vi.Next();
1087 if(vi.IsValid()) vi.FirstNotEmpty();
1088 }
1089 printf("Normalize(%i) %8i voxels\n",thr,loccnt);
1090 return loccnt;
1091 }
1092
1093
1094
1095 // Salva
SlicedPPMQ(const char * filename,const char * tag,int SliceNum)1096 void SlicedPPMQ( const char * filename,const char *tag,int SliceNum)
1097 {
1098 std::string basename=filename;
1099 std::string name;
1100 int ix,iy,iz;
1101 Color4b Tab[100];
1102 for(int ii=1;ii<100;++ii)
1103 Tab[ii].SetColorRamp(0,100,ii);
1104 Tab[0]=Color4b::Gray;
1105
1106 int ZStep=sz[2]/(SliceNum+1);
1107 for(iz=ZStep;iz<sz[2];iz+=ZStep)
1108 if(iz>=SubPartSafe.min[2] && iz<SubPartSafe.max[2])
1109 {
1110 name=SFormat("%s%03i%s_q.ppm",filename,iz,tag);
1111 FILE * fp = fopen(name.c_str(),"wb");
1112 fprintf(fp,
1113 "P6\n"
1114 "%d %d\n"
1115 "255\n"
1116 ,sz[1]
1117 ,sz[0]
1118 );
1119 unsigned char rgb[3];
1120 for(ix=0;ix<sz[0];++ix)
1121 {
1122 for(iy=0;iy<sz[1];++iy)
1123 {
1124 if( ix>=SubPartSafe.min[0] && ix<SubPartSafe.max[0] &&
1125 iy>=SubPartSafe.min[1] && iy<SubPartSafe.max[1])
1126 {
1127 if(!V(ix,iy,iz).B()) {
1128 rgb[0]=rgb[1]=rgb[2]=64;
1129 }
1130 else
1131 {
1132 float vv=V(ix,iy,iz).Q();
1133 int qi=std::min(V(ix,iy,iz).Q()*100.0f,100.0f);
1134
1135 if( vv>0) {
1136 rgb[0]=Tab[qi][0];
1137 rgb[1]=Tab[qi][1];
1138 rgb[2]=Tab[qi][2];
1139 }
1140 else if(vv<0)
1141 {
1142 rgb[0]=128;
1143 rgb[1]=255+32*vv;
1144 rgb[2]=0;//V(ix,iy,iz).Q()*256;
1145 }
1146 else {
1147 rgb[0]=255; rgb[1]=255; rgb[2]=255;
1148 }
1149 }
1150 }
1151 else{
1152 rgb[0]=rgb[1]=rgb[2]=64;
1153 }
1154 fwrite(rgb,3,1,fp);
1155 }
1156 }
1157 fclose(fp);
1158 }
1159 }
1160
1161 void SlicedPPM( const char * filename,const char *tag,int SliceNum=1)
1162 {
1163 std::string basename=filename;
1164 std::string name;
1165 int ix,iy,iz;
1166 int ZStep=sz[2]/(SliceNum+1);
1167 for(iz=ZStep;iz<sz[2];iz+=ZStep)
1168 if(iz>=SubPartSafe.min[2] && iz<SubPartSafe.max[2])
1169 {
1170 name=SFormat("%s_%03i_%s.ppm",filename,iz,tag);
1171 printf("Saving slice '%s'",name.c_str());
1172 FILE * fp = fopen(name.c_str(),"wb");
1173 if(!fp) return;
1174 fprintf(fp,
1175 "P6\n"
1176 "%d %d\n"
1177 "255\n"
1178 ,sz[1]
1179 ,sz[0]
1180 );
1181 unsigned char rgb[3];
1182 for(ix=0;ix<sz[0];++ix)
1183 {
1184 for(iy=0;iy<sz[1];++iy)
1185 {
1186 if( ix>=SubPartSafe.min[0] && ix<SubPartSafe.max[0] &&
1187 iy>=SubPartSafe.min[1] && iy<SubPartSafe.max[1])
1188 {
1189 if(!V(ix,iy,iz).B()) {
1190 rgb[0]=rgb[1]=rgb[2]=64;
1191 }
1192 else
1193 {
1194 float vv=V(ix,iy,iz).V();
1195 if( vv>0) {
1196 rgb[0]=255-32*vv;
1197 rgb[1]=128;
1198 rgb[2]=0;//V(ix,iy,iz).Q()*256;
1199 }
1200 else if(vv<0)
1201 {
1202 rgb[0]=128;
1203 rgb[1]=255+32*vv;
1204 rgb[2]=0;//V(ix,iy,iz).Q()*256;
1205 }
1206 else {
1207 rgb[0]=255; rgb[1]=255; rgb[2]=255;
1208 }
1209 }
1210 }
1211 else{
1212 rgb[0]=rgb[1]=rgb[2]=64;
1213 }
1214 fwrite(rgb,3,1,fp);
1215 }
1216 }
1217 fclose(fp);
1218 }
1219 }
1220 template < class VertexPointerType >
GetXIntercept(const vcg::Point3i & p1,const vcg::Point3i & p2,VertexPointerType & v,float)1221 void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, float /*thr*/)
1222 {
1223 float f1 = Val(p1.X(), p1.Y(), p1.Z());
1224 float f2 = Val(p2.X(), p2.Y(), p2.Z());
1225 float u = (float) f1/(f1-f2);
1226 v->P().X() = (float) p1.X()*(1-u) + u*p2.X();
1227 v->P().Y() = (float) p1.Y();
1228 v->P().Z() = (float) p1.Z();
1229 v->Q()=cV(p1.X(), p1.Y(), p1.Z()).Q();
1230 v->C()=cV(p1.X(), p1.Y(), p1.Z()).C4b();
1231 }
1232
1233 template < class VertexPointerType >
GetYIntercept(const vcg::Point3i & p1,const vcg::Point3i & p2,VertexPointerType & v,float)1234 void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, float /*thr*/)
1235 {
1236 float f1 = Val(p1.X(), p1.Y(), p1.Z());
1237 float f2 = Val(p2.X(), p2.Y(), p2.Z());
1238 float u = (float) f1/(f1-f2);
1239 v->P().X() = (float) p1.X();
1240 v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y();
1241 v->P().Z() = (float) p1.Z();
1242 v->Q()=cV(p1.X(), p1.Y(), p1.Z()).Q();
1243 v->C()=cV(p1.X(), p1.Y(), p1.Z()).C4b();
1244 }
1245
1246 template < class VertexPointerType>
GetZIntercept(const vcg::Point3i & p1,const vcg::Point3i & p2,VertexPointerType & v,float)1247 void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, float /*thr*/)
1248 {
1249 float f1 = Val(p1.X(), p1.Y(), p1.Z());
1250 float f2 = Val(p2.X(), p2.Y(), p2.Z());
1251 float u = (float) f1/(f1-f2);
1252 v->P().X() = (float) p1.X();
1253 v->P().Y() = (float) p1.Y();
1254 v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z();
1255 v->Q()=cV(p1.X(), p1.Y(), p1.Z()).Q();
1256 v->C()=cV(p1.X(), p1.Y(), p1.Z()).C4b();
1257 }
1258
1259 };
1260
1261
1262
1263 template < class VOL >
1264 class VolumeIterator
1265 {
1266 public:
1267 VOL &V;
1268 //vector<VOL::voxel_type> vi;
VolumeIterator(VOL & _VV)1269 VolumeIterator(VOL &_VV):V(_VV) {}
1270
1271 //Point3i curPos;
1272 int rpos;
1273 int lpos;
1274
Restart()1275 void Restart(){rpos=0;lpos=0;}
1276 private :
1277 public:
1278
1279
Set(const Point3i & p)1280 void Set(const Point3i &p)
1281 {
1282 //curPos=p;
1283 V.Pos(p[0],p[1],p[2],rpos,lpos);
1284 }
FirstNotEmpty()1285 bool FirstNotEmpty()
1286 {
1287
1288 //Dump();
1289 typename std::vector<std::vector<typename VOL::voxel_type> >::iterator rvi=V.rv.begin()+rpos;
1290 do
1291 {
1292 if((*rvi).empty())
1293 {
1294 while(rvi!=V.rv.end() && (*rvi).empty()) ++rvi;
1295 if(rvi==V.rv.end())
1296 {
1297 rpos=-1;
1298 return false;
1299 }
1300 rpos= rvi-V.rv.begin();
1301 lpos=0;
1302 }
1303 typename std::vector<typename VOL::voxel_type>::iterator lvi= (*rvi).begin()+lpos;
1304 // un voxel e' non vuoto se ha b!=0;
1305 while(lvi!=(*rvi).end() && !((*lvi).B() || (*lvi).Cnt()>0)) {
1306 ++lvi;
1307 }
1308 if(lvi!=(*rvi).end())
1309 {
1310 lpos= lvi-(*rvi).begin();
1311 //V.IPos(p[0],p[1],p[2],rpos,lpos);
1312 //Dump();
1313 return true;
1314 }
1315 else lpos=0;
1316 ++rvi;
1317 rpos= rvi-V.rv.begin();
1318
1319 } while (rvi!=V.rv.end());
1320 rpos=-1;
1321 return false;
1322 }
1323
1324 typename VOL::voxel_type &operator *()
1325 {
1326 assert(rpos>=0 && lpos >=0);
1327 return V.rv[rpos][lpos];
1328 }
Next()1329 bool Next()
1330 {
1331 assert(IsValid());
1332 if(lpos< VOL::BLOCKSIDE() * VOL::BLOCKSIDE() * VOL::BLOCKSIDE() -1)
1333 {
1334 ++lpos;
1335 //V.IPos(p[0],p[1],p[2],rpos,lpos);
1336 return true;
1337 }
1338 if(rpos < int(V.rv.size()-1))
1339 {
1340 lpos=0;
1341 ++rpos;
1342 //V.IPos(p[0],p[1],p[2],rpos,lpos);
1343 return true;
1344 }
1345 rpos=-1;
1346 lpos=-1;
1347 return false;
1348 }
1349
IsValid()1350 bool IsValid() {
1351 return rpos>=0;
1352 }
Dump()1353 void Dump()
1354 {
1355 int x,y,z;
1356 V.IPos(x,y,z,rpos,lpos);
1357 printf("Iterator r %4i l %4i (%3i %3i %3i)\n",rpos,lpos,x,y,z);
1358 }
1359
1360
1361 };
1362 }
1363 #endif
1364