diff --git a/apps/plymc/plymc.pro b/apps/plymc/plymc.pro new file mode 100755 index 00000000..9e5a3f0f --- /dev/null +++ b/apps/plymc/plymc.pro @@ -0,0 +1,10 @@ +TARGET = plymc +DEPENDPATH += . +INCLUDEPATH += ../.. +CONFIG += console c++11 +CONFIG -= app_bundle +TEMPLATE = app + +SOURCES += ../../wrap/ply/plylib.cpp \ + plymc_main.cpp + diff --git a/apps/plymc/plymc_main.cpp b/apps/plymc/plymc_main.cpp new file mode 100644 index 00000000..417ea42d --- /dev/null +++ b/apps/plymc/plymc_main.cpp @@ -0,0 +1,223 @@ +/* + * plymc_main.cpp + * filter_plymc + * + * Created by Paolo Cignoni on 10/23/09. + * Copyright 2009 ISTI - CNR. All rights reserved. + * + */ + +#include +#include "simplemeshprovider.h" +#define _PLYMC_VER "4.0" + +using namespace std; +using namespace vcg; + + + +string MYbasename = "plymcout"; +string VolumeBaseName; +string subtag=""; // la stringa da appendere al nome di file per generare i nomi di file relativi a vari sottopezzi +string alnfile; + +/************ Command Line Parameters *******/ + +int saveMask=vcg::tri::io::Mask::IOM_ALL; +void usage() +{ + printf( + "\nUsage:\n" + " plymc [options] filein.ply [filein.ply...]\n" + " plymc [options] filein.aln \n" + "Options: (no leading space before numeric values!) \n" + " -oname Set the base output name (default plymcout, without 'ply' extension)\n" + " -vname Enable the saving of the final volume with the specified filename\n" + " -C## Set numbers mesh that can be cached (default: 6)\n" + " -c## Set numbers of k cells (default: 10000)\n" + " -V# Set the required voxel size (override -c)\n" + " -s... Compute only a subvolume (specify 6 integers) \n" + " -S... Compute all the subvolumes of a partition (specify 3 int) \n" + " -X... Compute a range of the the subvolumes of a partition (specify 9 int)\n" + " -M Apply a 'safe' simplification step that removes only the unecessary triangles\n" + " -w# Set distance field Expansion factor in voxel (default 3)\n" + " -W# Set distance field Exp. as an absolute dist (override -w)\n" + " -a# Set angle threshold for distance field expansion (default 30)\n" + " -f# Set the fill threshold (default 12: all voxel having less \n" + " than 12 adjacent are not automaticall filled)\n" + " -L# Set Number of smoothing passes to be done after merging of all meshes\n" + " -R# Set Number of Refill passes to be done after merging of all meshes\n" + " -l# Make a single smoothing step after each expansion\n" + " -G# Disable Geodesic Quality\n" + " -F# Use per vertex quality defined in plyfile (geodesic is disabled)\n" + " -O# Set an Offset (<0!) threshold and build a double surface\n" + " -q# Set Quality threshold for smoothing. Only whose distance (in voxel)\n" + " is lower than the specified one are smoothed (default 3 voxel)\n" + " -Q# Same of above but expressed in absolute units.\n" + " -p use vertex splatting instead face rasterizing\n" + " -d# set as verbose level (default 0)\n" + " -D# save debug slices during processing\n" + + "\nNotes:\n\n" + "The Quality threshold can be expressed in voxel unit or in absolute units.\n" + "It represents the geodetic distance from the mesh border.\n" + "I.e. -q3 means that all voxels that are within 3voxel from the mesh border\n" + "are smoothed.\n\n" + + "A partition in the working volume is defined using 3 integers, that \n" + "specify the subdivision along the three axis.\n" + "To automatically compute ALL subvolumes of a given subdivision use '-S' \n" + "-S 1 1 1 default no subdivision at all\n" + "-S 2 2 2 compute all the octant of a octree-like subdivision\n" + "-S 1 1 4 compute four Z-slices\n\n" + + "To work only on a SINGLE subvolume of the partition you have to specify \n" + "six integers: the first three ints specify the subdivision along the\n" + "three axis and the last three ones which subvolume you desire.\n" + "the parameter to be used is '-s' \n" + "-s 1 1 1 0 0 0 default no subdivision at all\n" + "-s 1 1 3 0 0 1 make three Z-slices and take the middle one \n" + "-s 2 2 2 1 1 1 the last octant in a octree-like subdivision\n\n" + + "To START FROM a specific subvolume of the partition you have to specify \n" + "six integers: the first three ints specify the subdivision along the\n" + "three axis and the last three ones which subvolume you want to start\n" + "the process will go on using lexicographic order. Subvolumes are ordered\n" + "by Z, then by Y, then by X\n" + "The parameter to be used is '-K' \n" + "-K 4 4 4 0 0 0 a partition of 64 blocks, starting from the first one\n" + "-K 4 4 4 1 0 3 a partition of 64 blocks, starting from block 19 (index 1 0 3)\n\n" + + "To work only on a specific subvolume range of the partition you have \n" + "to specify nine integers: the first three ints specify the subdivision\n" + "along the three axis and, the next three which is the starting subvolume\n" + "and the last three which is the last subvolume to be computed.\n" + "the process will compute all blocks with x,y,z index included in the interval\n" + "specified: Xstart<=X<=Xend Ystart<=Y<=Yend Zstart<=Z<=Zend\n" + "-X 3 3 3 0 0 0 2 2 2 three subdivision on each axis, all subvolumes\n" + "-X 2 2 2 1 0 0 1 1 1 three subdivision on each axis, only the 'right' part\n\n" + ); + exit(-1); +} + + + +int main(int argc, char *argv[]) +{ + + Histogram h; + tri::PlyMC > pmc; + tri::PlyMC >::Parameter &p = pmc.p; + + + // This line is required to be sure that the decimal separatore is ALWAYS the . and not the , + // see the comment at the beginning of the file + setlocale(LC_ALL, "En_US"); + + printf( "\n PlyMC " _PLYMC_VER " (" __DATE__ ")\n" + " Copyright 2002-2016 Visual Computing Group I.S.T.I. C.N.R.\n" + " Paolo Cignoni (p.cignoni@isti.cnr.it)\n\n"); + //// Parameters + int i=1; + if(argc<2) usage(); + while(argv[i][0]=='-') + { + switch(argv[i][1]) + { + case 'o' : p.basename=argv[i]+2;printf("Setting Basename to %s\n",MYbasename.c_str());break; + case 'C' : pmc.MP.setCacheSize(atoi(argv[i]+2));printf("Setting MaxSize of MeshCache to %i\n",atoi(argv[i]+2)); break; + case 'c' : p.NCell =atoi(argv[i]+2);printf("Setting NCell to %i\n",p.NCell); break; + case 'v' : p.SaveVolumeFlag=true; VolumeBaseName=argv[i]+2; printf("Saving Volume enabled: volume Basename to %s\n",VolumeBaseName.c_str());break; + case 'V' : p.VoxSize =atof(argv[i]+2);printf("Setting VoxSize to %f; overridden NCell\n",p.VoxSize);p.NCell=0;break; + case 'w' : p.WideNum =atoi(argv[i]+2);printf("Setting WideNum to %i\n",p.WideNum);break; + case 'W' : p.WideSize=atof(argv[i]+2);printf("Setting WideSize to %f;overridden WideNum\n",p.WideSize);break; + case 'L' : p.SmoothNum =atoi(argv[i]+2);printf("Setting Laplacian SmoothNum to %i\n",p.SmoothNum);break; + case 'R' : p.RefillNum =atoi(argv[i]+2);printf("Setting Refilling Num to %i\n",p.RefillNum);break; + case 'q' : p.QualitySmoothVox=atof(argv[i]+2);printf("Setting QualitySmoothThr to %f; \n",p.QualitySmoothVox);break; + case 'Q' : p.QualitySmoothAbs=atof(argv[i]+2);printf("Setting QualitySmoothAbsolute to %f; it will override the default %f voxel value\n",p.QualitySmoothAbs,p.QualitySmoothVox);break; + case 'l' : p.IntraSmoothFlag=true; printf("Setting Laplacian Smooth after expansion \n");break; + case 'G' : p.GeodesicQualityFlag=false; printf("Disabling Geodesic Quality\n");break; + case 'F' : p.PLYFileQualityFlag=true; p.GeodesicQualityFlag=false; printf("Enabling PlyFile (and disabling Geodesic) Quality\n");break; + case 'f' : p.FillThr=atoi(argv[i]+2);printf("Setting Fill Threshold to %i\n",p.FillThr);break; + case 'a' : p.ExpAngleDeg=atoi(argv[i]+2);printf("Setting Expanding Angle Threshold to %f Degree\n",p.ExpAngleDeg);break; + case 'O' : p.OffsetThr=atof(argv[i]+2);printf("Setting Offset Threshold to %f \n",p.OffsetThr);p.OffsetFlag=true;break; + case 's' : + p.IDiv[0]=atoi(argv[++i]); p.IDiv[1]=atoi(argv[++i]); p.IDiv[2]=atoi(argv[++i]); + p.IPosS[0]=atoi(argv[++i]);p.IPosS[1]=atoi(argv[++i]);p.IPosS[2]=atoi(argv[++i]); + p.IPosE[0]=p.IPosS[0]; p.IPosE[1]=p.IPosS[1]; p.IPosE[2]=p.IPosS[2]; + if((p.IPosS[0]>=p.IDiv[0]) || (p.IPosS[1]>=p.IDiv[1]) || (p.IPosS[2]>=p.IDiv[2])) + { + printf("the subvolume you have requested is invalid (out of bounds)"); + exit(-1); + } + printf("Computing ONLY subvolume [%i,%i,%i] on a %ix%ix%i\n",p.IPosS[0],p.IPosS[1],p.IPosS[2],p.IDiv[0],p.IDiv[1],p.IDiv[2]); + break; + case 'S' : + p.IDiv[0]=atoi(argv[++i]);p.IDiv[1]=atoi(argv[++i]);p.IDiv[2]=atoi(argv[++i]); + p.IPosS=Point3i(0,0,0); + p.IPosE[0]=p.IDiv[0]-1; p.IPosE[1]=p.IDiv[1]-1; p.IPosE[2]=p.IDiv[2]-1; + printf("Autocomputing ALL subvolumes on a %ix%ix%i\n",p.IDiv[0],p.IDiv[1],p.IDiv[2]); + break; + case 'K' : + p.IDiv[0]=atoi(argv[++i]); p.IDiv[1]=atoi(argv[++i]);p.IDiv[2]=atoi(argv[++i]); + p.IPosB[0]=atoi(argv[++i]);p.IPosB[1]=atoi(argv[++i]);p.IPosB[2]=atoi(argv[++i]); + p.IPosS=Point3i(0,0,0); + p.IPosE[0]=p.IDiv[0]-1; p.IPosE[1]=p.IDiv[1]-1; p.IPosE[2]=p.IDiv[2]-1; + if((p.IPosB[0]>=p.IDiv[0]) || (p.IPosB[1]>=p.IDiv[1]) || (p.IPosB[2]>=p.IDiv[2])) + { + printf("the start subvolume you have requested is invalid (out of bounds)"); + exit(-1); + } + printf("Autocomputing ALL subvolumes FROM [%i,%i,%i] on a %ix%ix%i\n",p.IPosB[0],p.IPosB[1],p.IPosB[2],p.IDiv[0],p.IDiv[1],p.IDiv[2]); + break; + case 'X' : + p.IDiv[0]=atoi(argv[++i]); p.IDiv[1]=atoi(argv[++i]);p.IDiv[2]=atoi(argv[++i]); + p.IPosS[0]=atoi(argv[++i]);p.IPosS[1]=atoi(argv[++i]);p.IPosS[2]=atoi(argv[++i]); + p.IPosE[0]=atoi(argv[++i]);p.IPosE[1]=atoi(argv[++i]);p.IPosE[2]=atoi(argv[++i]); + // test if the interval is ok + int Istart,Iend; + Istart = p.IPosS[2] + (p.IPosS[1]*p.IDiv[2]) + (p.IPosS[0]*p.IDiv[2]*p.IDiv[1]); + Iend = p.IPosE[2] + (p.IPosE[1]*p.IDiv[2]) + (p.IPosE[0]*p.IDiv[2]*p.IDiv[1]); + if((Iend-Istart)<=0) + { + printf("the range you have requested is invalid (reversed or empty)"); + exit(-1); + } + if((p.IPosS[0]>=p.IDiv[0]) || (p.IPosS[1]>=p.IDiv[1]) || (p.IPosS[2]>=p.IDiv[2]) || + (p.IPosE[0]>=p.IDiv[0]) || (p.IPosE[1]>=p.IDiv[1]) || (p.IPosE[2]>=p.IDiv[2])) + { + printf("the subvolume you have requested is invalid (out of bounds)"); + exit(-1); + } + printf("Autocomputing subvolumes FROM [%i,%i,%i] TO [%i,%i,%i] on a %ix%ix%i\n",p.IPosS[0],p.IPosS[1],p.IPosS[2],p.IPosE[0],p.IPosE[1],p.IPosE[2],p.IDiv[0],p.IDiv[1],p.IDiv[2]); + break; + // case 'B' : p.SafeBorder =atoi(argv[i]+2);printf("Setting SafeBorder among blocks to %i*%i (default 1)\n",p.SafeBorder,Volume::BLOCKSIDE());break; + case 'p' : p.VertSplatFlag =true; printf("Enabling VertexSplatting instead of face rasterization\n");break; + case 'd' : p.VerboseLevel=atoi(argv[i]+2);printf("Enabling VerboseLevel= %i )\n",p.VerboseLevel);break; + case 'D' : p.VerboseLevel=1; p.SliceNum=atoi(argv[i]+2);printf("Enabling Debug Volume saving of %i slices (VerboseLevel=1)\n",p.SliceNum);break; + case 'M' : p.SimplificationFlag =true; printf("Enabling PostReconstruction simplification\n"); break; + default : {printf("Error unable to parse option '%s'\n",argv[i]); exit(0);} + } + ++i; + } + + + + Matrix44f Identity; + Identity.SetIdentity(); + string alnfile; + while(i +#include +#include +#include +#include + +using namespace std; +using namespace vcg; + +template + class MeshCache +{ + class Pair + { + public: + Pair(){used=0;} + TriMeshType *M; + std::string Name; + int used; // 'data' dell'ultimo accesso. si butta fuori quello lru + }; + + std::list MV; + +public: + void clear(); + + MeshCache() {MaxSize=6;} + ~MeshCache() { + typename std::list::iterator mi; + for(mi=MV.begin();mi!=MV.end();++mi) + delete (*mi).M; + } +// Restituisce true se la mesh e' in cache; +// restituisce in ogni caso il puntatore dove sta (o dovrebbe stare) la mesh +// Gestione LRU + + bool Find(std::string &name, TriMeshType * &sm) + { + typename std::list::iterator mi; + typename std::list::iterator oldest; // quello che e' piu' tempo che non viene acceduto. + int last; + + last = std::numeric_limits::max(); + oldest = MV.begin(); + + for(mi=MV.begin();mi!=MV.end();++mi) + { + if((*mi).usedMaxSize) { + sm=(*oldest).M; + (*oldest).used=0; + (*oldest).Name=name; + } else { + MV.push_back(Pair()); + MV.back().Name=name; + MV.back().M=new TriMeshType(); + sm=MV.back().M; + } + return false; +} + + + size_t MaxSize; + size_t size() const {return MV.size();} +}; + +template + class SimpleMeshProvider +{ + std::vector< std::string > meshnames; + std::vector TrV; + std::vector WV; // vettore con i pesi da applicare alla mesh. + std::vector BBV; // vettore con i bbox trasformati delle mesh da scannare. + vcg::Box3f fullBBox; + MeshCache MC; + + public: + + int size() {return meshnames.size();} + + int getCacheSize() {return MC.MaxSize;} + int setCacheSize(size_t newsize) + { + if(newsize == MC.MaxSize) + return MC.MaxSize; + if(newsize <= 0) + return MC.MaxSize; + + MC.MaxSize = newsize; + return newsize; + } + + bool openALN (const char* alnName) + { +// vector rmaps; +// ALNParser::ParseALN(rmaps, alnName); + +// for(size_t i=0; i::FileExtension(meshnames[i],"PLY") || tri::io::Importer::FileExtension(meshnames[i],"ply")) + { + if(!(TrV[i]==Id)) + ret=ply::ScanBBox(meshnames[i].c_str(),BBV[i],mt,true,0); + else + ret=vcg::ply::ScanBBox(meshnames[i].c_str(),BBV[i]); + + } + else + { printf("Trying to import a non-ply file %s\n",meshnames[i].c_str());fflush(stdout); + TriMeshType m; + int retVal=tri::io::Importer::Open(m,meshnames[i].c_str()); + ret = (retVal==0); + tri::UpdateBounding::Box(m); + BBV[i].Import(m.bbox); + } + if( ! ret) + { + printf("\n\nwarning:\n file '%s' not found\n",meshnames[i].c_str());fflush(stdout); + continue; + } + fullBBox.Add(BBV[i]); + } + return true; + } + +}; + + class SVertex; + class SFace; + class SUsedTypes: public vcg::UsedTypes < vcg::Use::AsVertexType, + vcg::Use::AsFaceType >{}; + +class SVertex : public Vertex< SUsedTypes, vertex::Coord3f, vertex::Normal3f,vertex::VFAdj, vertex::BitFlags, vertex::Color4b, vertex::Qualityf>{}; +class SFace : public Face< SUsedTypes, face::VertexRef, face::Normal3f,face::Qualityf, face::VFAdj, face::BitFlags> {}; +class SMesh : public vcg::tri::TriMesh< std::vector< SVertex>, std::vector< SFace > > {}; + +#endif // SIMPLEMESHPROVIDER_H diff --git a/vcg/complex/algorithms/create/plymc/plymc.h b/vcg/complex/algorithms/create/plymc/plymc.h index 14fb0cc5..ed43bfe7 100644 --- a/vcg/complex/algorithms/create/plymc/plymc.h +++ b/vcg/complex/algorithms/create/plymc/plymc.h @@ -54,6 +54,8 @@ #include #include +#include + // local optimization #include #include diff --git a/vcg/complex/algorithms/create/plymc/svoxel.h b/vcg/complex/algorithms/create/plymc/svoxel.h new file mode 100644 index 00000000..13af65b5 --- /dev/null +++ b/vcg/complex/algorithms/create/plymc/svoxel.h @@ -0,0 +1,174 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004-2016 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +#ifndef __SVOXEL_H__ +#define __SVOXEL_H__ +template +class SVoxel +{ + +private: + short v; // distanza dalla superficie espressa in centesimi di voxel; + short q; // come distanza dal bordo espressa in centesimi di voxel; // nota che questo implica che non si potrebbe rasterizzare una singola mesh in un volume di lato > 3000^3 + char n[3]; + unsigned char cnt; + // se != 0 b deve essere true; + // b e'tenuto implicitamente usando cnt; + // se cnt == 0 || cnt>0 -> b=false + // se cnt == 255 -> b=true + +// b==false cnt==0 totalmente non inzializzato (Zero) +// b==false cnt >0 da normalizzare +// b==true cnt==0 gia' normalizzato +// b==true cnt >0 Errore!!! + + + +public: + typedef SCALAR_TYPE scalar; + SVoxel(SCALAR_TYPE vv, bool /*bb*/, Point3 &nn, scalar qq) {SetV(vv);SetN(nn);SetQ(qq);} + SVoxel(SCALAR_TYPE vv, const Point3 &nn, scalar qq) {SetV(vv);SetN(nn);SetQ(qq);cnt=255;} + bool B() const {return cnt==255;} // puo' essere a true solo se cnt == 0; (il che significa che e' stato gia' normalizzato + + void SetB(bool val) { + assert( val == (cnt==255 || cnt==0) ); + if(val) cnt=255; + else if(cnt==255) cnt=0; + } + + int Cnt() { + if(cnt==255) return 0; + else return int(cnt); + } + void SetCnt(int val) { cnt=(unsigned char)val;} + + + Point3 N() const { + return Point3(scalar(n[0])/127.0f,scalar(n[1])/127.0f,scalar(n[2])/127.0f); + } + scalar N(const int i) const + { + return scalar(n[i])/127.0f; + } + + void SetN(const Point3 &nn) + { + n[0]=char( nn[0]*127.0f ); + n[1]=char( nn[1]*127.0f ); + n[2]=char( nn[2]*127.0f ); + } + scalar V() const + { + return scalar(v)/100.0f; + } + + inline void Blend( SVoxel const & vx, scalar w) + { + float w1=1.0-w; + SetV(V()*w1+vx.V()*w); + SetQ(Q()*w1+vx.Q()*w); + SetN(N()*w1+vx.N()*w); + //return *this; + } + + void SetV(const float &vv) + { + v= short(vv*100.0f); + if(v==0) v=1; + } + + scalar Q() const + { + return scalar(q)/20.0f; + } + + void SetQ(const float &qq) + { + int qi = qq * 20.0f; + if (qi>32767) qi =32767; + q = qi; + } + + void Merge(const SVoxel &VOX) + { + v=(v*Q()+VOX.Q()*VOX.v)/(Q()+VOX.Q()); + SetQ(Q()+VOX.Q()); + } + void Set(const SVoxel &VOX) + { + v=VOX.v; + n[0]=VOX.n[0]; + n[1]=VOX.n[1]; + n[2]=VOX.n[2]; + q=VOX.q; + + } + + inline SVoxel & operator += ( SVoxel const & vx) + { + if(cnt==0) + { + v=vx.v; + q=vx.q; + n[0]=vx.n[0]; + n[1]=vx.n[1]; + n[2]=vx.n[2]; + cnt=1; + } + else + { + const int cnt1=cnt+1; + v+=vx.v; + q = (q*cnt+vx.q)/(cnt1); + n[0]=(vx.n[0]*int(cnt) +vx.n[0])/(cnt1) ; + n[1]=(vx.n[1]*int(cnt) +vx.n[1])/(cnt1) ; + n[2]=(vx.n[2]*int(cnt) +vx.n[2])/(cnt1) ; + if(cnt==255) cnt=1; + else ++cnt; + } + return *this; + } + + inline bool Normalize(int thr) + { + assert(cnt>0); + if(cnt + +class TriEdgeCollapseMCParameter : public BaseParameterClass +{ +public: + Box3f bb; + bool preserveBBox; + float areaThr; + void SetDefaultParams() + { + bb.SetNull(); + preserveBBox=true; + areaThr=0; + } + + TriEdgeCollapseMCParameter() {SetDefaultParams();} +}; + + + +template +// Specialized Simplification classes for removing all small pieces on meshes. +class MCTriEdgeCollapse: public tri::TriEdgeCollapse< MCTriMesh, VertexPair, MYTYPE> { + + public: + typedef typename MCTriMesh::VertexPointer VertexPointer; + typedef typename MCTriMesh::FaceType FaceType; + typedef typename MCTriMesh::VertexType::CoordType CoordType; + typedef typename MCTriMesh::VertexType::ScalarType ScalarType; + + inline MCTriEdgeCollapse( const VertexPair &p, int mark, BaseParameterClass *pp) + { + this->localMark = mark; + this->pos=p; + this->_priority = ComputePriority(pp); + } + + + // The priority is simply the edge lenght, + // but we consider collapsing edges that lies on one plane of the MC grid. + virtual inline ScalarType ComputePriority(BaseParameterClass *_pp) + { + TriEdgeCollapseMCParameter *pp=(TriEdgeCollapseMCParameter *)_pp; + const CoordType & p0 = this->pos.V(0)->cP(); + const CoordType & p1 = this->pos.V(1)->cP(); + + int diffCnt=0; + if( (p0[0] != p1[0]) ) diffCnt++; + if( (p0[1] != p1[1]) ) diffCnt++; + if( (p0[2] != p1[2]) ) diffCnt++; + + // non MC plane collapse return highest cost + // if(diffCnt>2) return this->_priority=std::numeric_limits::max() ; + if(pp->preserveBBox) + { + const Box3f &bb=pp->bb; + // collapse on the bbox border. Avoid it. + if(p0[0]==bb.min[0] || p0[0]==bb.max[0]) return this->_priority=std::numeric_limits::max() ; + if(p0[1]==bb.min[1] || p0[1]==bb.max[1]) return this->_priority=std::numeric_limits::max() ; + if(p0[2]==bb.min[2] || p0[2]==bb.max[2]) return this->_priority=std::numeric_limits::max() ; + if(p1[0]==bb.min[0] || p1[0]==bb.max[0]) return this->_priority=std::numeric_limits::max() ; + if(p1[1]==bb.min[1] || p1[1]==bb.max[1]) return this->_priority=std::numeric_limits::max() ; + if(p1[2]==bb.min[2] || p1[2]==bb.max[2]) return this->_priority=std::numeric_limits::max() ; + } + return this->_priority=Distance(p0,p1); + } + + inline void Execute(MCTriMesh &m,BaseParameterClass *) + { + const CoordType & p0 = this->pos.V(0)->cP(); + const CoordType & p1 = this->pos.V(1)->cP(); + std::vector starVec0; + std::vector starVec1; +// int count0=0,count1=0; + + vcg::face::VVStarVF(this->pos.V(0),starVec0); +// for(size_t i=0;icP()[0]) ) count0++; +// if( (p0[1]==starVec[i]->cP()[1]) ) count0++; +// if( (p0[2]==starVec[i]->cP()[2]) ) count0++; +// } + vcg::face::VVStarVF(this->pos.V(1),starVec1); +// for(size_t i=0;icP()[0]) ) count1++; +// if( (p1[1]==starVec[i]->cP()[1]) ) count1++; +// if( (p1[2]==starVec[i]->cP()[2]) ) count1++; +// } + CoordType MidPoint=(p0+p1)/2.0; + + if(starVec0.size()>starVec1.size()) MidPoint=p0; + if(starVec0.size()pos, MidPoint); + vcg::tri::EdgeCollapser::Do(m, this->pos, MidPoint); + } + + +}; + +#endif diff --git a/vcg/complex/algorithms/create/plymc/volume.h b/vcg/complex/algorithms/create/plymc/volume.h new file mode 100644 index 00000000..d543c32e --- /dev/null +++ b/vcg/complex/algorithms/create/plymc/volume.h @@ -0,0 +1,1379 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004-2016 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +#ifndef __VOLUME_H__ +#define __VOLUME_H__ + +#ifdef __MINGW32__ +#define _int64 __int64 +#endif + +#include "voxel.h" +#include "svoxel.h" +#include +#include + +//#define BLOCKSIDE() 8 + +// Stato di un voxel + +// B() dice se ci sono dati in uno stadio usabile. +// Cnt() dice quanti ce ne sono stati sommati (per la normalizzazione) + +// b==false cnt==0 totalmente non inzializzato (Zero) +// b==false cnt >0 da normalizzare +// b==true cnt==0 gia' normalizzato +// b==true cnt >0 Errore!!! + +// forward definition +template < class VOL > +class VolumeIterator; + +//****************************************** +//****************************************** +//typedef Voxel Voxelf; + +const char *SFormat( const char * f, ... ) + { + static char buf[4096]; + va_list marker; + va_start( marker, f ); + vsprintf(buf,f,marker); + va_end( marker ); + return buf; + } + + +template +class Volume { +public: + typedef SCALAR_TYPE scalar; + typedef Point3 Point3x; + typedef Box3 Box3x; + + typedef VOX_TYPE voxel_type; + + static int BLOCKSIDE() { return 8;} + Volume(){ + SetDefaultParam(); + } + // I dati veri e propri + // Sono contenuti in un vettore di blocchi. + std::vector< std::vector > rv; + Box3x bbox; + + _int64 AskedCells; + Point3x dim; /// Dimensione spaziale (lunghezza lati) del bbox + + Point3i sz; /// Dimensioni griglia come numero di celle per lato + + Point3i ssz; /// Dimensioni sottoblocco in esame come numero di celle per lato + + Point3i rsz; /// Dimensioni macro griglia dei blocchi in cui e' suddiviso il volume (ogni blocco e' BLOCKSIDE()^3 celle) + + Point3i asz; /// Dimensioni macro griglia dei blocchi relativa al sottoblocco in questione (quello allocato davvero!) + + + Point3x voxel; /// Dimensioni di una cella + + + int WN,WP; // di quanti vox mi devo allargare per settare la manhattan distance in neg e pos + int DeltaVoxelSafe; // di quanti vox mi devo allargare per stare sicuro nel fare solo una sottoparte. + + const Point3i ISize() { return sz; } +private : + // offset e distanze varie precalcolate una volta per tutte + Point3f nnf[26]; + Point3i nni[26]; + float len[26]; + float slen[26]; + + /// Gestione sottoparte + Point3i div; + Point3i pos; +public: + Box3i SubPart; // Sottoparte del volume da considerare ufficialmente + Box3x SubBox; // BBox della sottoparte del volume da considerare in coord assolute + Box3i SubPartSafe; // come sopra ma aumentati per sicurezza. + Box3x SubBoxSafe; + + FILE *LogFP; +bool Verbose; // se true stampa un sacco di info in piu su logfp; + + void SetDefaultParam(){ + WN=0; + WP=1; + //WN=-2;// + //WP=3; + DeltaVoxelSafe=BLOCKSIDE(); + Verbose=true; + LogFP=stderr; + } + + void Init(const Volume &VV) + { + SetDefaultParam(); + WN=VV.WN; + WP=VV.WP; + DeltaVoxelSafe=VV.DeltaVoxelSafe; + Init(VV.AskedCells,VV.bbox,VV.div,VV.pos); + } + + void Init(_int64 cells, Box3x bb, Point3i _div=Point3i(1,1,1), Point3i _pos=Point3i(0,0,0)) + { + Point3d voxdim;voxdim.Import(bb.max-bb.min); + AskedCells=cells; + vcg::BestDim( cells, voxdim, sz ); + bbox=bb; +/* + printf("grid of ~%i kcells: %d x %d x %d \n",int(cells/1000),sz[0],sz[1],sz[2]); + printf("grid voxel size of %f %f %f\n",voxdim[0]/sz[0],voxdim[1]/sz[1],voxdim[2]/sz[2]); +*/ + // il box deve essere multipli di BLOCKSIDE() + sz=((sz/BLOCKSIDE())+Point3i(1,1,1))*BLOCKSIDE(); + + + rsz=sz/BLOCKSIDE(); + if(sz!=rsz*BLOCKSIDE()) { + assert(0); // il box deve essere multipli di BLOCKSIDE() + exit(-1); + } + + dim=bbox.max-bbox.min; + voxel[0]=dim[0]/sz[0]; + voxel[1]=dim[1]/sz[1]; + voxel[2]=dim[2]/sz[2]; + + SetSubPart(_div,_pos); + ssz=SubPartSafe.max-SubPartSafe.min; + asz=ssz/BLOCKSIDE() + Point3i(1,1,1); + rv.clear(); + rv.resize(asz[0]*asz[1]*asz[2]); + for(size_t i=0;i=1) +pos indicano la coord del subbox da prendere in considerazione (sempre >=0 && < xdiv,ydiv,zdiv) +*/ + +void SetSubPart(Point3i _div, Point3i _pos) +{ + int k; + // Controllo correttezza parametri. + for(k=0;k<3;++k) + { + assert(_div[k]>0); + if(_div[k]==0){ + 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]); + exit(-1); + } + if(_pos[k]<0 || _pos[k]>=_div[k]){ + 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]); + exit(-1); + } + } + + div=_div; + pos=_pos; + + // Setting the subpart under analisys + for(k=0;k<3;++k) + { + SubPart.min[k]= pos[k]*sz[k]/div[k]; + SubPart.max[k]=(pos[k]+1)*sz[k]/div[k]; + SubBox.min[k]= bbox.min[k]+SubPart.min[k]*voxel[k]; + SubBox.max[k]= bbox.min[k]+SubPart.max[k]*voxel[k]; + } + + // Setting the Safe Subpart under analisys + SubPartSafe=SubPart; + for(k=0;k<3;++k) + { + SubPartSafe.min[k] -= DeltaVoxelSafe; + SubPartSafe.max[k] += DeltaVoxelSafe; + + if( SubPartSafe.min[k]< 0 ) SubPartSafe.min[k] = 0; + if( SubPartSafe.max[k]> sz[k] ) SubPartSafe.max[k] = sz[k]; + SubBoxSafe.min[k]= bbox.min[k]+SubPartSafe.min[k]*voxel[k]; + SubBoxSafe.max[k]= bbox.min[k]+SubPartSafe.max[k]*voxel[k]; + } +/* + printf(" Computing only subvolume: (%d x %d x %d)= %dk cells \n" + " %d,%d,%d -> %d,%d,%d\n" + ,SubPart.DimX(),SubPart.DimY(),SubPart.DimZ(),(int)(((__int64)SubPart.DimX()*(__int64)SubPart.DimY()*(__int64)SubPart.DimZ())/1000) + ,SubPart.min[0] ,SubPart.min[1] ,SubPart.min[2] + ,SubPart.max[0] ,SubPart.max[1] ,SubPart.max[2] ); +*/ +} + +public: + + // Sa + /*bool Write(string filename, const float &minv, const float &maxv ) + { + FILE *fp; + if(div!=Point3i(1,1,1)) { + string subvoltag; + GetSubVolumeTag(subvoltag); + filename+=subvoltag; + } + string datname=filename; + string rawname=filename; + datname+=".dat"; + rawname+=".raw"; + + fp=fopen(datname,"w"); + + fprintf(fp,"ObjectFileName: %s\n",rawname.c_str()); + fprintf(fp,"TaggedFileName: ---\n"); + 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]); + fprintf(fp,"SliceThickness: %f %f %f\n",voxel[2]/voxel[0],voxel[2]/voxel[1],voxel[2]/voxel[2]); + fprintf(fp,"Format: UCHAR\n"); + fprintf(fp,"NbrTags: 0\n"); + fprintf(fp,"ObjectType: TEXTURE_VOLUME_OBJECT\n"); + fprintf(fp,"ObjectModel: RGBA\n"); + fprintf(fp,"GridType: EQUIDISTANT\n"); + + fclose(fp); + fp=fopen(rawname,"wb"); + if(!fp) + { + printf("Error: unable ro open output volume file '%s'\n",filename); + return false; + } + + int i,j,k; + for(k=SubPart.min[2];k1) fv=1; + fv=((fv*2.0f)-1.0f)*127; + char fs= (char) fv; + fwrite(&fs,sizeof(char),1,fp); + } + fclose(fp); + return true; + }*/ + void AbsPos(Point3i pi, Point3x &p) + { + p[0]=bbox.min[0]+pi[0]*voxel[0]; + p[1]=bbox.min[1]+pi[1]*voxel[1]; + p[2]=bbox.min[2]+pi[2]*voxel[2]; + } + + + void GetSubVolumeTag(std::string &subtag) + { + char buf[32]; + if (div[0]<= 10 && div[1]<= 10 && div[2]<= 10 ) sprintf(buf,"_%01d%01d%01d",pos[0],pos[1],pos[2]); + else if(div[0]<= 100 && div[1]<= 100 && div[2]<= 100 ) sprintf(buf,"_%02d%02d%02d",pos[0],pos[1],pos[2]); + else sprintf(buf,"_%03d%03d%03d",pos[0],pos[1],pos[2]); + subtag=buf; + } + + /* + Data una posizione x,y,z restituisce true se tale posizione appartiene a un blocco gia' allocato + In ogni caso mette in rpos la posizione del subbloc e in lpos la posizione all'interno del sottoblocco + */ + bool Pos(const int &_x,const int &_y,const int &_z, int & rpos,int &lpos) const + { + int x=_x-SubPartSafe.min[0]; int y=_y-SubPartSafe.min[1]; int z=_z-SubPartSafe.min[2]; + + assert(_x>=SubPartSafe.min[0] && _x=SubPartSafe.min[1] && _y=SubPartSafe.min[2] && _z=0 && x=0 && y=0 && z=0 && rx=0 && ry=0 && rz=0 && lpos >=0); + + int rz = rpos / (asz[0]*asz[1]); int remainder = rpos % (asz[0]*asz[1]); + int ry = ( remainder ) / asz[0] ; + int rx = remainder % asz[0]; + + assert(rx>=0 && rx=0 && ry=0 && rz=0 && x=0 && y=0 && z0); + //Pos(x,y,z,trpos,tlpos); + //assert(trpos==rpos && tlpos == lpos); + return true; + } + + void Alloc(int rpos, const VOX_TYPE &zeroval) + { + rv[rpos].resize(BLOCKSIDE()*BLOCKSIDE()*BLOCKSIDE(),zeroval); + } + /************************************/ + // Funzioni di accesso ai dati + bool ValidCell(const Point3i &p1, const Point3i &p2) const + { + if(!cV(p1.X(),p1.Y(),p1.Z()).B() ) return false; + if(!cV(p2.X(),p1.Y(),p1.Z()).B() ) return false; + if(!cV(p1.X(),p2.Y(),p1.Z()).B() ) return false; + if(!cV(p2.X(),p2.Y(),p1.Z()).B() ) return false; + if(!cV(p1.X(),p1.Y(),p2.Z()).B() ) return false; + if(!cV(p2.X(),p1.Y(),p2.Z()).B() ) return false; + if(!cV(p1.X(),p2.Y(),p2.Z()).B() ) return false; + if(!cV(p2.X(),p2.Y(),p2.Z()).B() ) return false; + return true; + } + + float Val(const int &x,const int &y,const int &z) const { + if(!cV(x,y,z).B()) return 1000; + return cV(x,y,z).V(); + //else return numeric_limits::quiet_NaN( ); + } + + VOX_TYPE &V(const int &x,const int &y,const int &z) { + int rpos,lpos; + if(!Pos(x,y,z,rpos,lpos)) Alloc(rpos,VOX_TYPE::Zero()); + return rv[rpos][lpos]; + } + + const VOX_TYPE &cV(const int &x,const int &y,const int &z) const + { + int rpos,lpos; + if(!Pos(x,y,z,rpos,lpos)) return VOX_TYPE::Zero(); + else return rv[rpos][lpos]; + } + const VOX_TYPE &V(const int &x,const int &y,const int &z) const + { + int rpos,lpos; + if(!Pos(x,y,z,rpos,lpos)) return VOX_TYPE::Zero(); + else return rv[rpos][lpos]; + } + /************************************/ + void Fill(VOX_TYPE (__cdecl *func)(const Point3i &p) ) +{ + int x,y,z; + for(z=0;z &S, scalar SafeZone=1, scalar SafeQuality=0) +{ + if(sz!=S.sz) + { + printf("Error"); + exit(-1); + } + int lcnt=0; + VolumeIterator< Volume > vi(S); + vi.Restart(); + vi.FirstNotEmpty(); + // const Voxelf *VC; + while(vi.IsValid()) + // scandisci il volume in ingresso, per ogni voxel non vuoto del volume + // in ingresso calcola la media con gli adiacenti + { + if((*vi).B()) + { + int x,y,z; + IPos(x,y,z,vi.rpos,vi.lpos); + if(Bound1(x,y,z)) + { + VOX_TYPE &VC = V(x,y,z); + for(int i=0;i<26;++i) + { + VOX_TYPE &VV= S.V(x+nni[i][0],y+nni[i][1],z+nni[i][2]); + if(VV.B()) VC+=VV; + } + lcnt++; + + /* + Voxelf &VV=V(x,y,z); + //Voxelf &VV=rv[vi.rpos][vi.lpos]; + for(int i=0;i<26;++i) + { + VC = &(S.V(x+nni[i][0],y+nni[i][1],z+nni[i][2])); + if(VC->b) + { + VV+=*VC; + lcnt++; + } + }*/ + } + } + vi.Next(); + if(vi.IsValid()) vi.FirstNotEmpty(); + //if((lcnt%100)==0) vi.Dump(); + + } + // Step 2, + // dopo aver calcolato la media, + + VolumeIterator< Volume > svi(*this); + svi.Restart(); + svi.FirstNotEmpty(); + int smoothcnt=0; + int preservedcnt=0; + int blendedcnt=0; + const float FieldBorder = 1; // dove finisce la transizione tra la zona safe e quella smoothed + const float EndFBorderZone = SafeZone+FieldBorder; + const float EndQBorderZone = SafeQuality*1.5; + const float QBorder = EndQBorderZone-SafeQuality; // dove finisce la transizione tra la zona safe e quella smoothed + while(svi.IsValid()) + { + if((*svi).Cnt()>0) + { + VOX_TYPE &sv=S.rv[svi.rpos][svi.lpos]; + (*svi).Normalize(1); // contiene il valore mediato + float SafeThr = fabs(sv.V()); + + // Se la qualita' e' bassa o se siamo distanti si smootha sempre + // se siamo vicini allo zero e con buona qualita' si deve fare attenzione + if(SafeThr EndQBorderZone) + { // se il voxel corrente aveva un valore < safezone E qualita' > SafeQuality + // allora si copia il valore originale di S + if((SafeThr <= SafeZone) && sv.Q() > SafeQuality ) + { + (*svi)=sv; + (*svi).SetB(true); + ++preservedcnt; + } + else + { // Siamo nella zona di transizione o per field o per quality + float blendq= std::max(0.0f,std::min(1.0f,(EndQBorderZone-sv.Q())/QBorder)); + float blendf= std::max(0.0f,std::min(1.0f,(EndFBorderZone-SafeThr)/FieldBorder)); + float BlendFactor = 1.0-std::max(blendf,blendq); // quanto del voxel originale si prende; + (*svi).Blend(sv,BlendFactor); + ++blendedcnt; + } + } + ++smoothcnt; + } + svi.Next(); + if(svi.IsValid()) svi.FirstNotEmpty(); + } + + if(Verbose) fprintf(LogFP,"CopySmooth %i voxels, %i preserved, %i blended\n",smoothcnt,preservedcnt,blendedcnt); +} + +void Merge(Volume &S) +{ + VolumeIterator< Volume > svi(S); + svi.Restart(); + svi.FirstNotEmpty(); + int loccnt=0; + + while(svi.IsValid()) + { + if((*svi).B()) + { + int x,y,z; + IPos(x,y,z,svi.rpos,svi.lpos); + if(cV(x,y,z).B()) V(x,y,z).Merge( (*svi)); + else { + V(x,y,z).Set((*svi)); + V(x,y,z).SetB(true); + } + ++loccnt; + } + svi.Next(); + if(svi.IsValid()) svi.FirstNotEmpty(); + } + + printf("Merge2 %i voxels\n",loccnt); + +} + + +void Interize( Point3x & vert ) const // OK +{ + for(int j=0;j<3;++j) + { + assert(vert[j]>=bbox.min[j]); + assert(vert[j]<=bbox.max[j]); + vert[j] = (vert[j] - bbox.min[j]) * sz[j] / (bbox.max[j] - bbox.min[j]); + } +} + + +void DeInterize( Point3x & vert ) const // OK +{ + for(int j=0;j<3;++j) + vert[j] = vert[j] * (bbox.max[j] - bbox.min[j]) / sz[j] + bbox.min[j]; +} + +bool SplatVert( const Point3x & v0, double quality, const Point3x & nn, Color4b c) +{ + Box3i ibox; + + assert(math::Abs(SquaredNorm(nn) - 1.0) < 0.0001); // Just a safety check that the vertex normals are NORMALIZED! + ibox.min=Point3i(floor(v0[0]),floor(v0[1]),floor(v0[2])); + ibox.max=Point3i( ceil(v0[0]), ceil(v0[1]), ceil(v0[2])); + ibox.Intersect(SubPartSafe); + + ibox.max[0] = std::min(SubPartSafe.max[0]-1,ibox.max[0]); + ibox.max[1] = std::min(SubPartSafe.max[1]-1,ibox.max[1]); + ibox.max[2] = std::min(SubPartSafe.max[2]-1,ibox.max[2]); + + + // Skip faces not colliding current subvolume. + if(ibox.IsEmpty()) + { + // point outside the box do nothing + return false; + } + + Point3x iV, deltaIV; + + // Now scan the eight voxel surrounding the splat + // and fill them with the distance from the plane + for(iV[0]=ibox.min[0]; iV[0]<=ibox.max[0]; ++iV[0]) + for(iV[1]=ibox.min[1]; iV[1]<=ibox.max[1]; ++iV[1]) + for(iV[2]=ibox.min[2]; iV[2]<=ibox.max[2]; ++iV[2]) + { + deltaIV = v0-iV; + scalar offset = deltaIV * nn; + VOX_TYPE &VV=V(iV[0],iV[1],iV[2]); + VV=VOX_TYPE(offset,nn,quality,c); + } + return true; +} + +template + void 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) +{ + const scalar EPS = scalar(1e-12); + const int crd0 = CoordZ; + const int crd1 = (CoordZ+1)%3; + const int crd2 = (CoordZ+2)%3; + assert(fabs(norm[crd0])+0.001 > fabs(norm[crd1])); + assert(fabs(norm[crd0])+0.001 > fabs(norm[crd2])); + scalar x,y; + for(x=sx;x<=ex;++x) + for(y=sy;y<=ey;++y) + { + scalar n0 = (x-v0[crd1])*d10[crd2] - (y-v0[crd2])*d10[crd1]; + scalar n1 = (x-v1[crd1])*d21[crd2] - (y-v1[crd2])*d21[crd1]; + scalar n2 = (x-v2[crd1])*d02[crd2] - (y-v2[crd2])*d02[crd1]; + + if( (n0>-EPS && n1>-EPS && n2>-EPS) || + (n0< EPS && n1< EPS && n2< EPS ) ) + { + scalar iz = ( dist - x*norm[crd1] - y*norm[crd2] ) / norm[crd0]; + //assert(iz>=fbox.min[2] && iz<=fbox.max[2]); + AddIntercept(x,y,iz, quality, norm ); + } + } +} + +// Si sa che la faccia ha una intercetta sull'asse z-dir di coord xy alla posizione z; +// quindi si setta nei 2 vertici prima e 2 dopo la distanza corrispondente. +template + void AddIntercept( const int x, const int y, const scalar z, const scalar q, const Point3f &n ) +{ + scalar esgn = (n[CoordZ] > 0 ? -1 : 1); + int zint = floor(z); + scalar dist=z-zint; // sempre positivo e compreso tra zero e uno + + for(int k=WN;k<=WP;k++) + { + if(zint+k >= SubPartSafe.min[CoordZ] && zint+k < SubPartSafe.max[CoordZ]) + { + VOX_TYPE *VV; + if(CoordZ==2) VV=&V(x,y,zint+k); + if(CoordZ==1) VV=&V(y,zint+k,x); + if(CoordZ==0) VV=&V(zint+k,x,y); + scalar nvv= esgn*( k-dist); + if(!VV->B() || fabs(VV->V()) > fabs(nvv)) { + *VV=VOX_TYPE(nvv,n,q); + } + } + } +} + +// assume che i punti della faccia in ingresso siano stati interized +bool ScanFace2( const Point3x & v0, const Point3x & v1, const Point3x & v2, + scalar quality, const Point3x & norm)//, const int name ) // OK +{ + + Box3x fbox; // Bounding Box della faccia (double) + Box3i ibox; // Bounding Box della faccia (int) + int sx,sy,sz; // Bounding Box intero + int ex,ey,ez; + + // Calcolo bbox della faccia + fbox.Set(v0); + fbox.Add(v1); + fbox.Add(v2); + + // BBox intero (nota che il cast a int fa truncation (funge solo perche' v0,v1,v2 sono positivi) + + ibox.min[0] =sx = floor(fbox.min[0]); if( ((scalar)sx)!=fbox.min[0] ) ++sx; // necessario se il punto e'approx a .9999 + ibox.min[1] =sy = floor(fbox.min[1]); if( ((scalar)sy)!=fbox.min[1] ) ++sy; + ibox.min[2] =sz = floor(fbox.min[2]); if( ((scalar)sz)!=fbox.min[2] ) ++sz; + ibox.max[0] =ex = floor(fbox.max[0]); + ibox.max[1] =ey = floor(fbox.max[1]); + ibox.max[2] =ez = floor(fbox.max[2]); + // Skip faces not colliding current subvolume. + if(!ibox.Collide(SubPartSafe)) return false; + + Point3x d10 = v1 - v0; + Point3x d21 = v2 - v1; + Point3x d02 = v0 - v2; + + assert(norm.Norm() >0.999f && norm.Norm()<1.001f); +// assert(0); + scalar dist = norm * v0; + + + /**** Rasterizzazione bbox ****/ + + // Clamping dei valori di rasterizzazione al subbox corrente + sx = std::max(SubPartSafe.min[0],sx); ex = std::min(SubPartSafe.max[0]-1,ex); + sy = std::max(SubPartSafe.min[1],sy); ey = std::min(SubPartSafe.max[1]-1,ey); + sz = std::max(SubPartSafe.min[2],sz); ez = std::min(SubPartSafe.max[2]-1,ez); + + 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); + 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); + else RasterFace<2>(sx,ex,sy,ey,dist,norm,quality,v0,v1,v2,d10,d21,d02); + + +return true; +} + +// assume che i punti della faccia in ingresso siano stati interized +bool ScanFace( const Point3x & v0, const Point3x & v1, const Point3x & v2, + double quality, const Point3x & nn)//, const int name ) // OK +{ + const scalar EPS = scalar(1e-12); +// const scalar EPS_INT = scalar(1e-20); + const scalar EPS_INT = .3f;//scalar(1e-20); + +#ifndef _NDEBUG + if(quality==0.0) + { + printf("Zero quality face are not allowed!\n"); + exit(-1); + } +#endif + +// ++nfaces; + + Box3x fbox; // Bounding Box della faccia (double) + Box3i ibox; // Bounding Box della faccia (int) + int sx,sy,sz; // Bounding Box intero + int ex,ey,ez; + + // Calcolo bbox della faccia + fbox.Set(v0); + fbox.Add(v1); + fbox.Add(v2); + + // BBox intero (nota che il cast a int fa truncation (funge solo perche' v0,v1,v2 sono positivi) + + ibox.min[0] =sx = floor(fbox.min[0]); if( ((scalar)sx)!=fbox.min[0] ) ++sx; // necessario se il punto e'approx a .9999 + ibox.min[1] =sy = floor(fbox.min[1]); if( ((scalar)sy)!=fbox.min[1] ) ++sy; + ibox.min[2] =sz = floor(fbox.min[2]); if( ((scalar)sz)!=fbox.min[2] ) ++sz; + ibox.max[0] =ex = floor(fbox.max[0]); + ibox.max[1] =ey = floor(fbox.max[1]); + ibox.max[2] =ez = floor(fbox.max[2]); + // Skip faces not colliding current subvolume. + if(!ibox.Collide(SubPartSafe)) return false; + + /**** Dati per controllo intersezione ****/ + + // Versori degli spigoli della faccia + + Point3x d10 = v1 - v0; + Point3x d21 = v2 - v1; + Point3x d02 = v0 - v2; + + // Normale al piano della faccia e distanza origine + + Point3x norm = d10 ^ d21; + norm.Normalize(); + double dist = norm * v0; + + /**** Rasterizzazione bbox ****/ + + int x,y,z; + + + // Clamping dei valori di rasterizzazione al subbox corrente + sx = std::max(SubPartSafe.min[0],sx); ex = std::min(SubPartSafe.max[0]-1,ex); + sy = std::max(SubPartSafe.min[1],sy); ey = std::min(SubPartSafe.max[1]-1,ey); + sz = std::max(SubPartSafe.min[2],sz); ez = std::min(SubPartSafe.max[2]-1,ez); + + // Rasterizzazione xy + + if(fabs(norm[2])>EPS_INT) + for(x=sx;x<=ex;++x) + for(y=sy;y<=ey;++y) + { + double n0 = ((double)x-v0[0])*d10[1] - ((double)y-v0[1])*d10[0]; + double n1 = ((double)x-v1[0])*d21[1] - ((double)y-v1[1])*d21[0]; + double n2 = ((double)x-v2[0])*d02[1] - ((double)y-v2[1])*d02[0]; + + if( (n0>-EPS && n1>-EPS && n2>-EPS) || + (n0< EPS && n1< EPS && n2< EPS )) + { + double iz = ( dist - double(x)*norm[0] - double(y)*norm[1] ) / norm[2]; + //assert(iz>=fbox.min[2] && iz<=fbox.max[2]); + AddXYInt(x,y,iz,-norm[2], quality, nn ); + } + } + + // Rasterizzazione xz + + if(fabs(norm[1])>EPS_INT) + for(x=sx;x<=ex;++x) + for(z=sz;z<=ez;++z) + { + double n0 = ((double)x-v0[0])*d10[2] - ((double)z-v0[2])*d10[0]; + double n1 = ((double)x-v1[0])*d21[2] - ((double)z-v1[2])*d21[0]; + double n2 = ((double)x-v2[0])*d02[2] - ((double)z-v2[2])*d02[0]; + + if( (n0>-EPS && n1>-EPS && n2>-EPS) || + (n0< EPS && n1< EPS && n2< EPS )) + { + double iy = ( dist - double(x)*norm[0] - double(z)*norm[2] ) / norm[1]; + //assert(iy>=fbox.min[1] && iy<=fbox.max[1]); + AddXZInt(x,z,iy,-norm[1], quality,nn ); + } + } + + // Rasterizzazione yz + + if(fabs(norm[0])>EPS_INT) + for(y=sy;y<=ey;++y) + for(z=sz;z<=ez;++z) + { + double n0 = ((double)y-v0[1])*d10[2] - ((double)z-v0[2])*d10[1]; + double n1 = ((double)y-v1[1])*d21[2] - ((double)z-v1[2])*d21[1]; + double n2 = ((double)y-v2[1])*d02[2] - ((double)z-v2[2])*d02[1]; + + if( (n0>-EPS && n1>-EPS && n2>-EPS) || + (n0< EPS && n1< EPS && n2< EPS ) ) + { + double ix = ( dist - double(y)*norm[1] - double(z)*norm[2] ) / norm[0]; + //assert(ix>=fbox.min[0] && ix<=fbox.max[0]); + AddYZInt(y,z,ix,-norm[0], quality, nn ); + } + } + return true; +} +// Si sa che la faccia ha una intercetta sull'asse z-dir di coord xy alla posizione z; +// quindi si setta nei 2 vertici prima e 2 dopo la distanza corrispondente. + +void AddXYInt( const int x, const int y, const double z, const double sgn, const double q, const Point3f &n ) +{ double esgn = (sgn<0 ? -1 : 1);//*max(fabs(sgn),0.001); + double dist=z-floor(z); // sempre positivo e compreso tra zero e uno + int zint = floor(z); + for(int k=WN;k<=WP;k++) + if(zint+k >= SubPartSafe.min[2] && zint+k < SubPartSafe.max[2]) + { + VOX_TYPE &VV=V(x,y,zint+k); + double nvv= esgn*( k-dist); + if(!VV.B() || fabs(VV.V()) > fabs(nvv)) { + VV=VOX_TYPE(nvv,n,q); + } + } +} +void AddYZInt( const int y, const int z, const double x, const double sgn, const double q, const Point3f &n ) +{ double esgn = (sgn<0 ? -1 : 1);//*max(fabs(sgn),0.001); + double dist=x-floor(x); // sempre positivo e compreso tra zero e uno + int xint = int(floor(x)); + for(int k=WN;k<=WP;k++) + if(xint+k >= SubPartSafe.min[0] && xint+k < SubPartSafe.max[0]) + { + VOX_TYPE &VV=V(xint+k,y,z); + double nvv= esgn*( k-dist); + if(!VV.B() || fabs(VV.V()) > fabs(nvv)) { + VV=VOX_TYPE(nvv,n,q); + } + } +} +void AddXZInt( const int x, const int z, const double y, const double sgn, const double q, const Point3f &n ) +{ double esgn = (sgn<0 ? -1 : 1);//*max(fabs(sgn),0.001); + double dist=y-scalar(floor(y)); // sempre positivo e compreso tra zero e uno + int yint = floor(y); + for(int k=WN;k<=WP;k++) + if(yint+k >= SubPartSafe.min[1] && yint+k < SubPartSafe.max[1]) + { + VOX_TYPE &VV=V(x,yint+k,z); + double nvv= esgn*( k-dist); + if(!VV.B() || fabs(VV.V()) > fabs(nvv)) { + VV=VOX_TYPE(nvv,n,q); + } + } +} + + + + void Dump(FILE *fp) + { + fprintf(fp,"Volume Info:\n"); + 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]); + 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]); + fprintf(fp," Voxel dimension %7.4f %7.4f %7.4f \n",voxel[0],voxel[1],voxel[2]); + + 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); + fprintf(fp," Memory Info: \n Voxel Size %8li b Virtually needed mem %8i Mb\n", + sizeof(VOX_TYPE),int(sizeof(VOX_TYPE)*(_int64)(sz[0])*(_int64)(sz[1])*(_int64)(sz[2])/(1024*1024))); + if(div!=Point3i(1,1,1)) + { + fprintf(fp," Subdivided in %6i %6i %6i (tot: %12i block):\n",div[0],div[1],div[2],div[0]*div[1]*div[2]); + fprintf(fp," Computing subblock %6i %6i %6i :\n",pos[0],pos[1],pos[2]); + 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]); + 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]); + + } + fprintf(fp,"\n"); + } + + int Allocated() + {int cnt=0; + for(size_t i=0;iSubPartSafe.min[0] && x < SubPartSafe.max[0]-1 ) && + (y>SubPartSafe.min[1] && y < SubPartSafe.max[1]-1 ) && + (z>SubPartSafe.min[2] && z < SubPartSafe.max[2]-1 ) ; +} + +/* +Note sull'algoritmo di espansione: + +Si riempie i voxel vuoti + +Se il volume e' inizialmente riempito con ii valori delle intercette alla superficie +nei 2 vertici immediatamente adiacenti all'edge intersecato dalla superficie +si deve espadnere tale coampo in maniera sensata. + +Notare che e' importante che non tutto il campo sia riempito con un approx ella distanza di hausdorf: +il campo deve essere "tagliato" sui bordi della superficie per evitare pasticci. Levoy riempie il campo +solo lungo la direzione dello scanner, io invece riempio lungo la normale alla superficie. In questo modo +si evita il problema che l'espansione e' legata all'acquisizione iniziale + +*/ + + +void Expand(scalar AngleThrRad) +{ + int i; + VolumeIterator< Volume > vi(*this); + + float CosThr=math::Cos(AngleThrRad); +// printf("Expand2 angle %f, %f\n",AngleThrRad,CosThr); + int loccnt=0; + + vi.Restart(); + vi.FirstNotEmpty(); + while(vi.IsValid()) + { + if((*vi).B()) // si espande solo i voxel con valori "validi" + { + int x,y,z; + IPos(x,y,z,vi.rpos,vi.lpos); + Point3f n=(*vi).N(); + VOX_TYPE vtmp = (*vi); + if(Bound1(x,y,z)) + for(i=0;i<26;++i) + { + float angle = -(nnf[i]*n); // cos angolo tra la normale alla superficie e la direzione di espansione + if( fabs(angle)> CosThr ) + { + //bbfloat tt=(*vi).V(); + 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; + VOX_TYPE &VV= V(x+nni[i][0],y+nni[i][1],z+nni[i][2]); + if(!VV.B()){ + VV+=vtmp; + loccnt++; + } + } + } + } + vi.Next(); + if(vi.IsValid()) vi.FirstNotEmpty(); + } + printf("Expand %8i ",loccnt); + Normalize(1); +} + +// filla i buchi vuoti di un volume; +// scorre tutti i voxel pieni e aggiunge agli adiacenti vuoti tale valore; +// si riscorre tutti i voxel vuoti e se si sono sommati almeno thr valori sitiene. +// in pratica maggiore il valore di thr meno buchi si riempiono. + +void Refill(const int thr,float maxdistance = std::numeric_limits::max() ) +{ + int lcnt=0; + VolumeIterator< Volume > vi(*this); + vi.Restart(); + vi.FirstNotEmpty(); + + while(vi.IsValid()) + { + if((*vi).B()) + { + int x,y,z; + IPos(x,y,z,vi.rpos,vi.lpos); + if(Bound1(x,y,z)) + { + for(int i=0;i<26;++i) + { + VOX_TYPE &VC= V(x+nni[i][0],y+nni[i][1],z+nni[i][2]); + if(!VC.B()){ + if(VC.Cnt()==0) lcnt++; + VC+=(*vi); + } + } + } + } + + vi.Next(); + + if(vi.IsValid()) vi.FirstNotEmpty(); + + + } + printf("ReFill %8i ",lcnt); +Normalize(thr,maxdistance); +} + +/* +Attraversa il volume e modifica il valore di campo in modo da creare una offset surface che si connetta +bene con la superficie originale +il parametro specifica dove dovrebbe passare l'offset surface + +L'idea e' quella di fare un altro zero del distance field al threshold specificato + +La cosa deve essere smooth +quindi scelgo una funzione che abbia 2 zeri (in zero e in thr) e +*/ +void Offset(const float thr) +{ + int lcnt=0; + VolumeIterator< Volume > vi(*this); + vi.Restart(); + vi.FirstNotEmpty(); + float thr2=thr/2.0; + while(vi.IsValid()) + { + if((*vi).B()) + { + float vv=(*vi).V(); + if(thr<0) if(vv0) if(vv>thr2) vv=thr-vv; + + (*vi).SetV(vv); + } + + vi.Next(); + + if(vi.IsValid()) vi.FirstNotEmpty(); + + + } + printf("ReFill %8i ",lcnt); + Normalize(thr); +} + +// prende un volume e mette a true il campo b di tutti i voxel che hanno un valore significativo. +// thr indica il numero minimo di valori che si devono essere sommati sul voxel; +// di solito e' uguale a 1; +int Normalize(int thr, float maxdistance=std::numeric_limits::max() ) +{ + VolumeIterator< Volume > vi(*this); + vi.Restart(); + vi.FirstNotEmpty(); + int loccnt=0; + while(vi.IsValid()) + { + if(!(*vi).B()) + { + if((*vi).Normalize(thr)) + ++loccnt; + if(math::Abs((*vi).V())>maxdistance) *vi=VOX_TYPE::Zero(); + } + vi.Next(); + if(vi.IsValid()) vi.FirstNotEmpty(); + } + printf("Normalize(%i) %8i voxels\n",thr,loccnt); + return loccnt; +} + + + +// Salva +void SlicedPPMQ( const char * filename,const char *tag,int SliceNum) + { + std::string basename=filename; + std::string name; + int ix,iy,iz; + Color4b Tab[100]; + for(int ii=1;ii<100;++ii) + Tab[ii].SetColorRamp(0,100,ii); + Tab[0]=Color4b::Gray; + + int ZStep=sz[2]/(SliceNum+1); + for(iz=ZStep;iz=SubPartSafe.min[2] && iz=SubPartSafe.min[0] && ix=SubPartSafe.min[1] && iy0) { + rgb[0]=Tab[qi][0]; + rgb[1]=Tab[qi][1]; + rgb[2]=Tab[qi][2]; + } + else if(vv<0) + { + rgb[0]=128; + rgb[1]=255+32*vv; + rgb[2]=0;//V(ix,iy,iz).Q()*256; + } + else { + rgb[0]=255; rgb[1]=255; rgb[2]=255; + } + } + } + else{ + rgb[0]=rgb[1]=rgb[2]=64; + } + fwrite(rgb,3,1,fp); + } + } + fclose(fp); + } + } + +void SlicedPPM( const char * filename,const char *tag,int SliceNum=1) + { + std::string basename=filename; + std::string name; + int ix,iy,iz; + int ZStep=sz[2]/(SliceNum+1); + for(iz=ZStep;iz=SubPartSafe.min[2] && iz=SubPartSafe.min[0] && ix=SubPartSafe.min[1] && iy0) { + rgb[0]=255-32*vv; + rgb[1]=128; + rgb[2]=0;//V(ix,iy,iz).Q()*256; + } + else if(vv<0) + { + rgb[0]=128; + rgb[1]=255+32*vv; + rgb[2]=0;//V(ix,iy,iz).Q()*256; + } + else { + rgb[0]=255; rgb[1]=255; rgb[2]=255; + } + } + } + else{ + rgb[0]=rgb[1]=rgb[2]=64; + } + fwrite(rgb,3,1,fp); + } + } + fclose(fp); + } + } +template < class VertexPointerType > +void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, float /*thr*/) +{ + float f1 = Val(p1.X(), p1.Y(), p1.Z()); + float f2 = Val(p2.X(), p2.Y(), p2.Z()); + float u = (float) f1/(f1-f2); + v->P().X() = (float) p1.X()*(1-u) + u*p2.X(); + v->P().Y() = (float) p1.Y(); + v->P().Z() = (float) p1.Z(); + v->Q()=cV(p1.X(), p1.Y(), p1.Z()).Q(); + v->C()=cV(p1.X(), p1.Y(), p1.Z()).C4b(); +} + +template < class VertexPointerType > +void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, float /*thr*/) +{ + float f1 = Val(p1.X(), p1.Y(), p1.Z()); + float f2 = Val(p2.X(), p2.Y(), p2.Z()); + float u = (float) f1/(f1-f2); + v->P().X() = (float) p1.X(); + v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y(); + v->P().Z() = (float) p1.Z(); + v->Q()=cV(p1.X(), p1.Y(), p1.Z()).Q(); + v->C()=cV(p1.X(), p1.Y(), p1.Z()).C4b(); +} + +template < class VertexPointerType> +void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, float /*thr*/) +{ + float f1 = Val(p1.X(), p1.Y(), p1.Z()); + float f2 = Val(p2.X(), p2.Y(), p2.Z()); + float u = (float) f1/(f1-f2); + v->P().X() = (float) p1.X(); + v->P().Y() = (float) p1.Y(); + v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z(); + v->Q()=cV(p1.X(), p1.Y(), p1.Z()).Q(); + v->C()=cV(p1.X(), p1.Y(), p1.Z()).C4b(); +} + +}; + + + +template < class VOL > +class VolumeIterator +{ + public: + VOL &V; + //vector vi; + VolumeIterator(VOL &_VV):V(_VV) {} + + //Point3i curPos; + int rpos; + int lpos; + + void Restart(){rpos=0;lpos=0;} + private : + public: + + + void Set(const Point3i &p) + { + //curPos=p; + V.Pos(p[0],p[1],p[2],rpos,lpos); + } + bool FirstNotEmpty() + { + + //Dump(); + typename std::vector >::iterator rvi=V.rv.begin()+rpos; + do + { + if((*rvi).empty()) + { + while(rvi!=V.rv.end() && (*rvi).empty()) ++rvi; + if(rvi==V.rv.end()) + { + rpos=-1; + return false; + } + rpos= rvi-V.rv.begin(); + lpos=0; + } + typename std::vector::iterator lvi= (*rvi).begin()+lpos; + // un voxel e' non vuoto se ha b!=0; + while(lvi!=(*rvi).end() && !((*lvi).B() || (*lvi).Cnt()>0)) { + ++lvi; + } + if(lvi!=(*rvi).end()) + { + lpos= lvi-(*rvi).begin(); + //V.IPos(p[0],p[1],p[2],rpos,lpos); + //Dump(); + return true; + } + else lpos=0; + ++rvi; + rpos= rvi-V.rv.begin(); + + } while (rvi!=V.rv.end()); + rpos=-1; + return false; + } + + typename VOL::voxel_type &operator *() + { + assert(rpos>=0 && lpos >=0); + return V.rv[rpos][lpos]; + } + bool Next() + { + assert(IsValid()); + if(lpos< VOL::BLOCKSIDE() * VOL::BLOCKSIDE() * VOL::BLOCKSIDE() -1) + { + ++lpos; + //V.IPos(p[0],p[1],p[2],rpos,lpos); + return true; + } + if(rpos < int(V.rv.size()-1)) + { + lpos=0; + ++rpos; + //V.IPos(p[0],p[1],p[2],rpos,lpos); + return true; + } + rpos=-1; + lpos=-1; + return false; + } + + bool IsValid() { + return rpos>=0; + } + void Dump() + { + int x,y,z; + V.IPos(x,y,z,rpos,lpos); + printf("Iterator r %4i l %4i (%3i %3i %3i)\n",rpos,lpos,x,y,z); + } + + + +}; +#endif diff --git a/vcg/complex/algorithms/create/plymc/voxel.h b/vcg/complex/algorithms/create/plymc/voxel.h new file mode 100644 index 00000000..3dc18642 --- /dev/null +++ b/vcg/complex/algorithms/create/plymc/voxel.h @@ -0,0 +1,195 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004-2016 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ +#ifndef __VOXEL_H__ +#define __VOXEL_H__ +template +class Voxel +{ + public: + typedef SCALAR_TYPE scalar; + + Voxel(SCALAR_TYPE vv, bool bb, Point3 nn, short _cnt) {v=vv;b=bb;n=nn;cnt=_cnt;} + Voxel(SCALAR_TYPE vv, Point3 nn, scalar qq) {v=vv;b=true;n=nn;cnt=0;q=qq;} + + const scalar &N(const int i) const { return n[i]; } + + const Point3 &N() const { return n; } + + void SetN(const Point3 &nn) { n=nn; } + const scalar &V() const { return v; } + + void SetV(const scalar &vv) { v=vv; } + + const scalar &Q() const { return q; } + + void SetQ(const scalar &qq) { q=qq; } + + + bool B() const {return b;}; + void SetB(bool val) {b=val;} + int Cnt() const {return cnt;} + void SetCnt(int val) {cnt=val;} + inline void Blend( Voxel const & vx, scalar w) + { + float w1=1.0-w; + v=v*w1+vx.v*w; + q=q*w1+vx.q*w; + n=n*w1+vx.n*w; + //return *this; + } + + inline Voxel & operator += ( Voxel const & vx) + { + if(cnt==0) + { + assert(!b); + v=vx.v; + q=vx.q; + n=vx.n; + cnt=1; + b=false; + } + else + { + assert(!b); + v+=vx.v; + q+=vx.q; + n+=vx.n; + ++cnt; + } + return *this; + } + + inline bool Normalize(int thr) + { + assert(cnt>0); + assert(!B()); + if(cnt n; + +}; + + +class Voxelfc :public Voxel +{ +public: + + Voxelfc(float vv, bool bb, Point3f nn, short _cnt) :Voxel(vv,bb,nn,_cnt){} + Voxelfc(float vv, Point3f nn, scalar qq) :Voxel(vv,nn,qq) {} + Voxelfc(float vv, Point3f nn, scalar qq,Color4b cc) :Voxel(vv,nn,qq) + { + c[0]=cc[0]; + c[1]=cc[1]; + c[2]=cc[2]; + } + + inline bool Normalize(int thr) + { + if(cnt>=thr) c/=cnt; + return Voxel::Normalize(thr); + } + + static const Voxelfc &Zero() { + static Voxelfc tt(0,false,Point3f(0,0,0),0); + return tt; + } + + void Merge(const Voxelfc &VOX) + { + c=( c*q + VOX.C()*VOX.Q() )/(q+VOX.Q()); + Voxel::Merge(VOX); + } + + void Set(const Voxelfc &VOX) + { + Voxel::Set(VOX); + c=VOX.c; + } + + const float &C(const int i) const { return c[i]; } + const Point3f &C() const { return c; } + void SetC(const Point3f &cc) { c=cc; } + Color4b C4b() const + { + static Color4b cc; + cc=Color4b(c[0],c[1],c[2],255); + return cc; + } + inline void Blend( Voxelfc const & vx, scalar w) + { + float w1=1.0-w; + v=v*w1+vx.v*w; + q=q*w1+vx.q*w; + n=n*w1+vx.n*w; + c=c*w1+vx.c*w; + //return *this; + } + + inline Voxelfc & operator += ( Voxelfc const & vx) + { + Voxel::operator +=(vx); + if(cnt==1) c =vx.c; + else c+=vx.c; + return *this; + } + +private: + Point3f c; +}; + +#endif