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