diff --git a/apps/metro/defs.h b/apps/metro/defs.h new file mode 100644 index 00000000..e19f2487 --- /dev/null +++ b/apps/metro/defs.h @@ -0,0 +1,63 @@ +// ----------------------------------------------------------------------------------------------- +#ifndef _DEFS_H +#define _DEFS_H +// ----------------------------------------------------------------------------------------------- + +// command line parameters +#define CMD_LINE_ARG_HIST 'H' +#define CMD_LINE_ARG_VERTEX_SAMPLE 'V' +#define CMD_LINE_ARG_EDGE_SAMPLE 'E' +#define CMD_LINE_ARG_FACE_SAMPLE 'F' +#define CMD_LINE_ARG_SAMPLE_TYPE 'S' +#define CMD_LINE_ARG_MONTECARLO_SAMPLING 'M' +#define CMD_LINE_ARG_SUBDIVISION_SAMPLING 'S' +#define CMD_LINE_ARG_SIMILAR_TRIANGLES_SAMPLING 'T' +#define CMD_LINE_ARG_N_SAMPLES 'N' +#define CMD_LINE_ARG_SAMPLES_PER_AREA_UNIT 'A' +#define CMD_LINE_ARG_SAVE_DISPLACEMENT 'O' +#define CMD_LINE_ARG_SAVE_ERROR_AS_COLOUR 'C' + + +// error messages +#define MSG_ERR_N_ARGS "\nUsage: "\ + "metro file1 file2 [opt]\n"\ + "where opt can be:\n"\ + "-V disable vertex sampling\n"\ + "-E disable edge sampling\n"\ + "-F disable face sampling\n"\ + "-Sx set the face sampling mode\n"\ + " where x can be:\n"\ + " -M montecarlo sampling\n"\ + " -S subdivision sampling\n"\ + " -T similar triangles sampling\n"\ + "-N# set the required number of samples (overrides -A)\n"\ + "-A# set the required number of samples per area unit (overrides -N)\n"\ + "-Hxxx.yyy save histogram of error values in the file xxx.yyy\n"\ + "-O save error as vertex quality\n"\ + "-C save error as vertex colour"\ + "\n" + +#define MSG_ERR_MESH_LOAD "error loading the input meshes.\n" +#define MSG_ERR_INVALID_OPTION "unable to parse option '%s'\n" +#define MSG_ERR_FILE_OPEN "unable to open the output file.'n" +#define MSG_ERR_UNKNOWN_FORMAT "unknown file format '%s'.\n" + +// global constants +#define NO_SAMPLES_PER_FACE 10 +#define N_SAMPLES_EDGE_TO_FACE_RATIO 0.1 +#define BBOX_FACTOR 0.1 +#define INFLATE_PERCENTAGE 0.02 +#define MIN_SIZE 125 /* 125 = 5^3 */ +#define N_HIST_BINS 256 +#define PRINT_EVERY_N_ELEMENTS 1000 +#define FILE_EXT_SMF "smf" +#define FILE_EXT_PLY "ply" + +// strings +#define STR_HIST_FILENAME_DEFAULT "hist.txt" +#define STR_NEW_MESH_FILENAME_DEFAULT "error.ply" +#define STR_NEW_MESH_FILENAME_DEFAULT_2 "error_colour.ply" + +// ----------------------------------------------------------------------------------------------- +#endif +// ----------------------------------------------------------------------------------------------- diff --git a/apps/metro/error.h b/apps/metro/error.h new file mode 100644 index 00000000..d38ec347 --- /dev/null +++ b/apps/metro/error.h @@ -0,0 +1,49 @@ +#ifndef __METRO__VCG__HAUS +#define __METRO__VCG__HAUS + +#include "sampling.h" + +template +double Distance(MESH_TYPE & S1,MESH_TYPE & S2){ + int flags = 0x18; + int n_samples_target = 1.5 * __max(S1.fn,S2.fn); + + double dist12_max,dist21_max; + // compute face information + tri::UpdateEdges::Set(S1); + tri::UpdateEdges::Set(S2); + + // set bounding boxes for S1 and S2 + tri::UpdateBounding::Box(S1); + tri::UpdateBounding::Box(S2); + + // set Bounding Box. + Box3 bbox, tmp_bbox_M1=S1.bbox, tmp_bbox_M2=S2.bbox; + bbox.Add(S1.bbox); + bbox.Add(S2.bbox); + bbox.InflateFix(0.1); + S1.bbox = bbox; + S2.bbox = bbox; + + Sampling ForwardSampling(S1,S2); + Sampling BackwardSampling(S2,S1); + + ForwardSampling.SetFlags(flags); + BackwardSampling.SetFlags(flags); + + ForwardSampling.SetSamplesTarget(n_samples_target); + BackwardSampling.SetSamplesTarget(n_samples_target); + + ForwardSampling.Hausdorff(); + dist12_max = ForwardSampling.GetDistMax(); + + BackwardSampling.Hausdorff(); + dist21_max = BackwardSampling.GetDistMax(); + + return math::Max(dist12_max,dist12_max); + + } + +#endif + + diff --git a/apps/metro/mesh_type.h b/apps/metro/mesh_type.h new file mode 100644 index 00000000..9354d278 --- /dev/null +++ b/apps/metro/mesh_type.h @@ -0,0 +1,56 @@ +/***************************************************************************** + * VCGLib * + * * + * Visual Computing Group o> * + * IEI Institute, CNUCE Institute, CNR Pisa <| * + * / \ * + * Copyright(C) 1999 by Paolo Cignoni, Claudio Rocchini * + * All rights reserved. * + * * + * Permission to use, copy, modify, distribute and sell this software and * + * its documentation for any purpose is hereby granted without fee, provided * + * that the above copyright notice appear in all copies and that both that * + * copyright notice and this permission notice appear in supporting * + * documentation. the author makes no representations about the suitability * + * of this software for any purpose. It is provided "as is" without express * + * or implied warranty. * + * * + *****************************************************************************/ +/**************************************************************************** + History + 2003 Dic 17 modifiche per conversione alla versione template (gano) +****************************************************************************/ + +// ----------------------------------------------------------------------------------------------- +#ifndef _CMESH_H +#define _CMESH_H +// ----------------------------------------------------------------------------------------------- + +#pragma warning(disable:4786 4804 4666) +// standard libraries. +#include + +// VCG library. +//#include +#include +#include +#include +#include + + + +using namespace vcg; +using namespace std; + +// Vertex, Face, Mesh and Grid definitions. +class MyEdge; +class CFace; +class CVertex : public VertexCMNQ< double,MyEdge,CFace > {}; +class CFace : public FaceRTFMFN< CVertex > {public: + CFace*& F(const int j){ return (CFace*&) FaceRTFMFN::F(j);} + }; +class CMesh : public tri::TriMesh< vector, vector > {}; + +// ----------------------------------------------------------------------------------------------- +#endif +// ----------------------------------------------------------------------------------------------- diff --git a/apps/metro/metro.cpp b/apps/metro/metro.cpp new file mode 100644 index 00000000..0ec2ec57 --- /dev/null +++ b/apps/metro/metro.cpp @@ -0,0 +1,375 @@ +// ----------------------------------------------------------------------------------------------- + +// standard libraries +#include + +// project definitions. +#include "defs.h" +#include "sampling.h" +#include "mesh_type.h" +#include +#include +//#include +#include +#include + + + +// ----------------------------------------------------------------------------------------------- + + +////////////////// Command line Flags and parameters +bool ComputeHistFlag = false; +bool VertexSampleFlag = true; +bool EdgeSampleFlag = true; +bool FaceSampleFlag = true; +bool MontecarloSamplingFlag = false; +bool SubdivisionSamplingFlag = false; +bool SimilarTrianglesSamplingFlag = false; +bool NumberOfSamples = false; +bool SamplesPerAreaUnit = false; +bool SaveErrorDisplacement = false; +bool SaveErrorAsColour = false; + +// ----------------------------------------------------------------------------------------------- + + +inline char* GetExtension(char* filename) +{ + for(int i=strlen(filename)-1; i >= 0; i--) + if(filename[i] == '.') + break; + if(i > 0) + return &(filename[i+1]); + else + return NULL; +} + + +void main(int argc, char**argv) +{ + CMesh S1, S2; + double dist1_max, dist1_mean, dist1_RMS, volume_1; + double dist2_max, dist2_mean, dist2_RMS, volume_2; + double mesh_dist_max; + unsigned long n_samples_target, n_samples_output, elapsed_time; + double n_samples_per_area_unit; + int flags, flags_fwd, flags_back, n_samples_area, n_samples_edge, n_samples_vertex, err; + char *fmt, *hist_filename, *new_mesh_filename, *new_mesh_filename_2; + char fname_1[] = STR_NEW_MESH_FILENAME_DEFAULT, fname_2[] = STR_NEW_MESH_FILENAME_DEFAULT_2; + FILE *fd; + + // print program info + printf("-------------------------------\n" + " Metro\n" + " release date: "__DATE__"\n" + "-------------------------------\n\n"); + + // load input meshes. + if(argc <= 2) + { + printf(MSG_ERR_N_ARGS); + exit(-1); + } + + // load mesh M1. + if(!(fmt = GetExtension(argv[1]))) + { + printf(MSG_ERR_UNKNOWN_FORMAT, fmt); + exit(-1); + } + if(!_stricmp(FILE_EXT_PLY, fmt)) + { + printf("reading the mesh `%s'...", argv[1]); + err = tri::io::ImporterPLY::Open(S1,argv[1]); + } + /* else + if(!_stricmp(FILE_EXT_SMF, fmt)) + { + printf("reading the mesh `%s'...", argv[1]); + err = tri::io::ImporterSMF::Open(S1,argv[1]); + }*/ + else + { + printf(MSG_ERR_UNKNOWN_FORMAT, fmt); + exit(-1); + } + if(err < 0) + { + printf("\n"); + printf(MSG_ERR_MESH_LOAD); + exit(-1); + } + else + printf("done\n"); + + // load mesh M2. + if(!(fmt = GetExtension(argv[2]))) + { + printf(MSG_ERR_UNKNOWN_FORMAT, fmt); + exit(-1); + } + if(!_stricmp(FILE_EXT_PLY, fmt)) + { + printf("reading the mesh `%s'...", argv[2]); + err = tri::io::ImporterPLY::Open(S2,argv[2]); + } + /*else + if(!_stricmp(FILE_EXT_SMF, fmt)) + { + printf("reading the mesh `%s'...", argv[2]); + err = S2.Load_Smf(argv[2]); + }*/ + else + { + printf(MSG_ERR_UNKNOWN_FORMAT, fmt); + exit(-1); + } + if(err < 0) + { + printf(MSG_ERR_MESH_LOAD); + exit(-1); + } + else + printf("done\n"); + + // parse command line. + for(int i=3; i < argc;) + { + if(argv[i][0]=='-') + switch(argv[i][1]) + { + case CMD_LINE_ARG_HIST : ComputeHistFlag = true; hist_filename = &(argv[i][2]); + if(hist_filename[0] == '\0') + strcpy(hist_filename, STR_HIST_FILENAME_DEFAULT); + break; + case CMD_LINE_ARG_VERTEX_SAMPLE : VertexSampleFlag = false; break; + case CMD_LINE_ARG_EDGE_SAMPLE : EdgeSampleFlag = false; break; + case CMD_LINE_ARG_FACE_SAMPLE : FaceSampleFlag = false; break; + case CMD_LINE_ARG_SAMPLE_TYPE : + switch(argv[i][2]) + { + case CMD_LINE_ARG_MONTECARLO_SAMPLING : MontecarloSamplingFlag = true; break; + case CMD_LINE_ARG_SUBDIVISION_SAMPLING : SubdivisionSamplingFlag = true; break; + case CMD_LINE_ARG_SIMILAR_TRIANGLES_SAMPLING : SimilarTrianglesSamplingFlag = true; break; + default : printf(MSG_ERR_INVALID_OPTION, argv[i]); + exit(0); + } + break; + case CMD_LINE_ARG_N_SAMPLES : NumberOfSamples = true; n_samples_target = atoi(&(argv[i][2])); break; + case CMD_LINE_ARG_SAMPLES_PER_AREA_UNIT : SamplesPerAreaUnit = true; n_samples_per_area_unit = (double) atof(&(argv[i][2])); break; + case CMD_LINE_ARG_SAVE_DISPLACEMENT : SaveErrorDisplacement = true; new_mesh_filename = &(argv[i][2]); + if(new_mesh_filename[0] == '\0') + new_mesh_filename = fname_1; + break; + case CMD_LINE_ARG_SAVE_ERROR_AS_COLOUR : SaveErrorAsColour = true; new_mesh_filename_2 = &(argv[i][2]); + if(new_mesh_filename_2[0] == '\0') + new_mesh_filename_2 = fname_2; + break; + default : printf(MSG_ERR_INVALID_OPTION, argv[i]); + exit(0); + } + i++; + } + + // set sampling scheme + int sampling_method = MontecarloSamplingFlag + SubdivisionSamplingFlag + SimilarTrianglesSamplingFlag; + // defaults + if(!sampling_method) + SimilarTrianglesSamplingFlag = true; + if(sampling_method > 1) + { + printf("Cannot choose more than one sampling method. Similar Triangles sampling assumed.\n"); + SimilarTrianglesSamplingFlag = true; + MontecarloSamplingFlag = false; + SubdivisionSamplingFlag = false; + } + if(!NumberOfSamples && !SamplesPerAreaUnit) + { + NumberOfSamples = true; + n_samples_target = NO_SAMPLES_PER_FACE * max(S1.fn,S2.fn); + } + + // compute face information + tri::UpdateEdges::Set(S1); + tri::UpdateEdges::Set(S2); + + // set bounding boxes for S1 and S2 + tri::UpdateBounding::Box(S1); + tri::UpdateBounding::Box(S2); + + // set Bounding Box. + Box3d bbox, tmp_bbox_M1=S1.bbox, tmp_bbox_M2=S2.bbox; + bbox.Add(S1.bbox); + bbox.Add(S2.bbox); + bbox.InflateFix(INFLATE_PERCENTAGE); + S1.bbox = bbox; + S2.bbox = bbox; + + // set flags. + flags = 0; + if(ComputeHistFlag) + flags |= FLAG_HIST; + if(VertexSampleFlag) + flags |= FLAG_VERTEX_SAMPLING; + if(EdgeSampleFlag) + flags |= FLAG_EDGE_SAMPLING; + if(FaceSampleFlag) + flags |= FLAG_FACE_SAMPLING; + if(MontecarloSamplingFlag) + flags |= FLAG_MONTECARLO_SAMPLING; + if(SubdivisionSamplingFlag) + flags |= FLAG_SUBDIVISION_SAMPLING; + if(SimilarTrianglesSamplingFlag) + flags |= FLAG_SIMILAR_TRIANGLES_SAMPLING; + flags_fwd = flags; + flags_back = flags; + if(SaveErrorDisplacement) + { + if(S1.vn >= S2.vn) + flags_fwd |= FLAG_SAVE_ERROR_DISPLACEMENT; + else + flags_back |= FLAG_SAVE_ERROR_DISPLACEMENT; + } + + if(SaveErrorAsColour) + { + if(S1.vn >= S2.vn) + flags_fwd |= FLAG_SAVE_ERROR_AS_COLOUR; + else + flags_back |= FLAG_SAVE_ERROR_AS_COLOUR; + } + + // initialize time info. + int t0=clock(); + + // print mesh info. + Sampling ForwardSampling(S1,S2); + Sampling BackwardSampling(S2,S1); + + printf("Mesh info:\n"); + printf(" M1: '%s'\n\t%vertices %7i\n\tfaces %7i\n\tarea %12.4f\n", argv[1], S1.vn, S1.fn, ForwardSampling.GetArea()); + printf("\tbbox (%7.4f %7.4f %7.4f)-(%7.4f %7.4f %7.4f)\n", tmp_bbox_M1.min[0], tmp_bbox_M1.min[1], tmp_bbox_M1.min[2], tmp_bbox_M1.max[0], tmp_bbox_M1.max[1], tmp_bbox_M1.max[2]); + printf("\tbbox diagonal %f\n", (float)tmp_bbox_M1.Diag()); + printf(" M2: '%s'\n\t%vertices %7i\n\tfaces %7i\n\tarea %12.4f\n", argv[2], S2.vn, S2.fn, BackwardSampling.GetArea()); + printf("\tbbox (%7.4f %7.4f %7.4f)-(%7.4f %7.4f %7.4f)\n", tmp_bbox_M2.min[0], tmp_bbox_M2.min[1], tmp_bbox_M2.min[2], tmp_bbox_M2.max[0], tmp_bbox_M2.max[1], tmp_bbox_M2.max[2]); + printf("\tbbox diagonal %f\n", (float)tmp_bbox_M2.Diag()); + + // Forward distance. + printf("\nForward distance (M1 -> M2):\n"); + ForwardSampling.SetFlags(flags_fwd); + if(NumberOfSamples) + { + ForwardSampling.SetSamplesTarget(n_samples_target); + n_samples_per_area_unit = ForwardSampling.GetNSamplesPerAreaUnit(); + } + else + { + ForwardSampling.SetSamplesPerAreaUnit(n_samples_per_area_unit); + n_samples_target = ForwardSampling.GetNSamplesTarget(); + } + printf("target # samples : %u\ntarget # samples/area : %f\n", n_samples_target, n_samples_per_area_unit); + ForwardSampling.Hausdorff(); + dist1_max = ForwardSampling.GetDistMax(); + dist1_mean = ForwardSampling.GetDistMean(); + dist1_RMS = ForwardSampling.GetDistRMS(); + volume_1 = ForwardSampling.GetDistVolume(); + n_samples_output = ForwardSampling.GetNSamples(); + n_samples_area = ForwardSampling.GetNAreaSamples(); + n_samples_edge = ForwardSampling.GetNEdgeSamples(); + n_samples_vertex = ForwardSampling.GetNVertexSamples(); + printf("\ndistance:\n max : %f (%f with respect to bounding box diagonal)\n mean : %f\n RMS : %f\n", (float)dist1_max, (float)dist1_max/bbox.Diag(), (float)dist1_mean, (float)dist1_RMS); + if(VertexSampleFlag) + printf("# vertex samples %d\n", n_samples_vertex); + if(EdgeSampleFlag) + printf("# edge samples %d\n", n_samples_edge); + printf("# area samples %d\n# total samples %d\nsamples per area unit: %f\n\n", n_samples_area, n_samples_output, ForwardSampling.GetNSamplesPerAreaUnit()); + + // Backward distance. + printf("\nBackward distance (M2 -> M1):\n"); + BackwardSampling.SetFlags(flags_back); + if(NumberOfSamples) + { + BackwardSampling.SetSamplesTarget(n_samples_target); + n_samples_per_area_unit = BackwardSampling.GetNSamplesPerAreaUnit(); + } + else + { + BackwardSampling.SetSamplesPerAreaUnit(n_samples_per_area_unit); + n_samples_target = BackwardSampling.GetNSamplesTarget(); + } + printf("target # samples : %u\ntarget # samples/area : %f\n", n_samples_target, n_samples_per_area_unit); + BackwardSampling.Hausdorff(); + dist2_max = BackwardSampling.GetDistMax(); + dist2_mean = BackwardSampling.GetDistMean(); + dist2_RMS = BackwardSampling.GetDistRMS(); + volume_2 = BackwardSampling.GetDistVolume(); + n_samples_output = BackwardSampling.GetNSamples(); + n_samples_area = BackwardSampling.GetNAreaSamples(); + n_samples_edge = BackwardSampling.GetNEdgeSamples(); + n_samples_vertex = BackwardSampling.GetNVertexSamples(); + printf("\ndistance:\n max : %f (%f with respect to bounding box diagonal)\n mean : %f\n RMS : %f\n", (float)dist2_max, (float)dist2_max/bbox.Diag(), (float)dist2_mean, (float)dist2_RMS); + if(VertexSampleFlag) + printf("# vertex samples %d\n", n_samples_vertex); + if(EdgeSampleFlag) + printf("# edge samples %d\n", n_samples_edge); + printf("# area samples %d\n# total samples %d\nsamples per area unit: %f\n\n", n_samples_area, n_samples_output, BackwardSampling.GetNSamplesPerAreaUnit()); + + // compute time info. + elapsed_time = clock() - t0; + + // save error distribution histogram + /*if(ComputeHistFlag) + { + const Hist &hist1 = ForwardSampling.GetHist(); + const Hist &hist2 = BackwardSampling.GetHist(); + if(!(fd = fopen(hist_filename, "w"))) + { + printf(MSG_ERR_FILE_OPEN); + exit(-1); + } + vector::const_iterator ii; + vector::const_iterator fi; + + fprintf(fd, "error distribution histogram (forward distance)\n\n"); + for(ii=hist1.H.begin(), fi=hist1.R.begin(); ii != hist1.H.end(); ++fi,ii++) + fprintf(fd, "%6.4f\t%d\n", *fi, *ii); + + fprintf(fd, "\n\nerror distribution histogram (backward distance)\n"); + for(ii=hist2.H.begin(), fi=hist2.R.begin(); ii != hist2.H.end(); ++fi,ii++) + fprintf(fd, "%6.4f\t%d\n", *fi, *ii); + + fclose(fd); + }*/ + + // max distance. + mesh_dist_max = max(dist1_max , dist2_max); + printf("\nHausdorff distance: %f (%f with respect to bounding box diagonal)\nComputation time : %d ms\n# samples/second : %f\n\n", (float)mesh_dist_max, (float)mesh_dist_max/bbox.Diag(), (int)elapsed_time, (float)n_samples_output/(float)elapsed_time*2000.0F); + + // save error files. + if((flags_fwd & FLAG_SAVE_ERROR_DISPLACEMENT) && (flags_fwd & FLAG_SAVE_ERROR_AS_COLOUR)) + if(!strcmp(new_mesh_filename, new_mesh_filename_2)) + { + tri::io::ExporterPLY::Save( S1,new_mesh_filename); + exit(0); + } + if((flags_back & FLAG_SAVE_ERROR_DISPLACEMENT) && (flags_back & FLAG_SAVE_ERROR_AS_COLOUR)) + if(!strcmp(new_mesh_filename, new_mesh_filename_2)) + { + tri::io::ExporterPLY::Save( S2,new_mesh_filename_2); + exit(0); + } + + //if(flags_fwd & FLAG_SAVE_ERROR_DISPLACEMENT) + // S1.SavePly(new_mesh_filename, CMesh::SM_ALL & (CMesh::SM_ALL ^ CMesh::SM_VERTCOLOR)); + //else + // if(flags_back & FLAG_SAVE_ERROR_DISPLACEMENT) + // S2.SavePly(new_mesh_filename, CMesh::SM_ALL & (CMesh::SM_ALL ^ CMesh::SM_VERTCOLOR)); + //if(flags_fwd & FLAG_SAVE_ERROR_AS_COLOUR) + // S1.SavePly(new_mesh_filename_2, CMesh::SM_ALL & (CMesh::SM_ALL ^ CMesh::SM_VERTQUALITY)); + //else + // if(flags_back & FLAG_SAVE_ERROR_AS_COLOUR) + // S2.SavePly(new_mesh_filename_2, CMesh::SM_ALL & (CMesh::SM_ALL ^ CMesh::SM_VERTQUALITY)); +} + +// ----------------------------------------------------------------------------------------------- diff --git a/apps/metro/min_dist_point.h b/apps/metro/min_dist_point.h new file mode 100644 index 00000000..f50a9ea8 --- /dev/null +++ b/apps/metro/min_dist_point.h @@ -0,0 +1,206 @@ +/*#*************************************************************************** + * MinDistPoint.h o o * + * o o * + * Visual Computing Group _ O _ * + * IEI Institute, CNUCE Institute, CNR Pisa \/)\/ * + * /\/| * + * Copyright(C) 1999 by Paolo Cignoni, Paolo Pingi, Claudio Rocchini | * + * All rights reserved. \ * + * * + * Permission to use, copy, modify, distribute and sell this software and * + * its documentation for any purpose is hereby granted without fee, provided * + * that the above copyright notice appear in all copies and that both that * + * copyright notice and this permission notice appear in supporting * + * documentation. the author makes no representations about the suitability * + * of this software for any purpose. It is provided "as is" without express * + * or implied warranty. * + * * + * NOTE THAT THIS FILE SHOULD NOT DIRECTL BE INCLUDED * + * It is automatically included by Mesh.h * + * * + ***************************************************************************#*/ +/*#************************************************************************** + History + + 2000 Nov 06 First Working release (pc) + 08 Aggiunto if(gr.bbox.IsIn(p)) per evitare piantamenti se si chiede + un punto fuori. + Tolto un Normalize inutile + 2001 May 17 Aggiunta versione della Mindistpoint che da anche le coord + baricentriche del punto trovato (pc); aggiunto wrapper per la vecchia + versione. + Dec 10 Corretto prodotto scalare vettore nell'ordine giusto in un paio di posti + 2002 Mar 29 Templatata anche in funzione del tipo scalare. (pc) + Oct 24 Corretti warning (unsigned mismatch) del vc7, + 2003 Apr 15 Corretti mismatch iterator / pointer + Jun 08 Aggiornato UGridLink -> UGrid::Link + 2004 Gen 09 Aggiunte le inclusioni a Point3[4].h + 2004 Gen 19 Corretto qualche ->Normal() in ->cN() +****************************************************************************/ +#ifndef __VCG_MINDISTPOINT +#define __VCG_MINDISTPOINT +#include + +#include +#include +#include +using namespace vcg; + + +/* +aka MetroCore +data una mesh m e una ug sulle sue facce trova il punto di m piu' vicino ad +un punto dato. +*/ + +// input: mesh, punto, griglia, distanza limite +// output: normale alla faccia e punto piu' vicino su di essa + +// Nota che il parametro template GRID non ci dovrebbe essere, visto che deve essere +// UGrid, ma non sono riuscito a definirlo implicitamente + +template +void MinDistPoint( MESH & mesh, const Point3 & p, GRID & gr, SCALAR & mdist, + Point3 & normf, Point3 & bestq, typename MESH::FaceType * &f, Point3 &ip) +{ + typedef SCALAR scalar; + typedef Point3 Point3x; + typedef Box3 Box3x; + + if(!gr.bbox.IsIn(p)) return; + typedef GridStaticObj::Link A2UGridLink; + scalar ax = p[0] - gr.bbox.min[0]; // Real coodinate of point refer to + scalar ay = p[1] - gr.bbox.min[1]; + scalar az = p[2] - gr.bbox.min[2]; + + int gx = int( ax/gr.voxel[0] ); // Integer coordinate of the point + int gy = int( ay/gr.voxel[1] ); // voxel + int gz = int( az/gr.voxel[2] ); + + scalar vx = gr.bbox.min[0]+gx*gr.voxel[0]; // Real world coordinate of the Voxel + scalar vy = gr.bbox.min[1]+gy*gr.voxel[1]; // origin + scalar vz = gr.bbox.min[2]+gz*gr.voxel[2]; + + scalar dx = math::Min(p[0] - vx, vx+gr.voxel[0]-p[0]); // Dist from the voxel + scalar dy = math::Min(p[1] - vy, vy+gr.voxel[1]-p[1]); + scalar dz = math::Min(p[2] - vz, vz+gr.voxel[2]-p[2]); + + scalar vdist,vstep; + + if(dxElem())) ) + { + if( (*(l->Elem())).Dist( p, error, q) ) + { + bestq = q; + bestf = l->Elem(); + MESH::ScalarType alfa=1, beta=1, gamma=1; + + //bestf->InterpolationParameters(q, alfa, beta); + //calcolo normale con interpolazione trilineare + /*normf = (1-(alfa+beta))*(bestf->V(0)->Normal())+ + (alfa*(bestf->V(1)->Normal()))+ + (beta*(bestf->V(2)->Normal()));*/ + bool ret=bestf->InterpolationParameters(q, alfa, beta, gamma); + //assert(ret); + normf = (bestf->V(0)->cN())*alfa+ + (bestf->V(1)->cN())*beta+ + (bestf->V(2)->cN())*gamma; + normf.Normalize(); + ip[0]=alfa;ip[1]=beta;ip[2]=gamma; + + } + + mesh.Mark( &*(l->Elem()) ); + } + } + else + { + for(int ix=gx-s;ix<=gx+s;++ix) + if( ix>=0 && ix=0 && iy=0 && izElem())) ) + { + if( (*(l->Elem())).Dist( p, error, q) ) + { + bestq = q; + bestf = l->Elem(); + MESH::ScalarType alfa, beta, gamma; + //bestf->InterpolationParameters(q, alfa, beta); + //calcolo normale con interpolazione trilineare + bestf->InterpolationParameters(q, alfa, beta, gamma); + normf = (bestf->V(0)->cN())*alfa+ + (bestf->V(1)->cN())*beta+ + (bestf->V(2)->cN())*gamma ; + ip[0]=alfa;ip[1]=beta;ip[2]=gamma; + //normf.Normalize(); inutile si assume le normali ai vertici benfatte + } + mesh.Mark(&*l->Elem()); + } + } + } + } + } + + if( fabs(error) +void MinDistPoint( MESH & mesh, const Point3 & p, GRID & gr, SCALAR & mdist, + Point3 & normf, Point3 & bestq, typename MESH::face_type * &f) +{ + Point3 ip; + MinDistPoint(mesh,p,gr,mdist,normf,bestq,f,ip); +} +#endif \ No newline at end of file diff --git a/apps/metro/sampling.h b/apps/metro/sampling.h new file mode 100644 index 00000000..a89a28ab --- /dev/null +++ b/apps/metro/sampling.h @@ -0,0 +1,607 @@ +/***************************************************************************** + * VCGLib * + * * + * Visual Computing Group o> * + * IEI Institute, CNUCE Institute, CNR Pisa <| * + * / \ * + * Copyright(C) 1999 by Paolo Cignoni, Claudio Rocchini * + * All rights reserved. * + * * + * Permission to use, copy, modify, distribute and sell this software and * + * its documentation for any purpose is hereby granted without fee, provided * + * that the above copyright notice appear in all copies and that both that * + * copyright notice and this permission notice appear in supporting * + * documentation. the author makes no representations about the suitability * + * of this software for any purpose. It is provided "as is" without express * + * or implied warranty. * + * * + *****************************************************************************/ +/**************************************************************************** + History + 2003 Dic 17 modifiche per conversione alla versione template (gano) + 2004 Jan 19 qualche ->P() in ->cP() (Snarf) +****************************************************************************/ + +// ----------------------------------------------------------------------------------------------- +#ifndef _SAMPLING_H +#define _SAMPLING_H +// ----------------------------------------------------------------------------------------------- + +// standard libraries. +#include + +// VCG library. +//#include +#include "min_dist_point.h" +//#include +#include +#include +#include +using namespace vcg; +// ----------------------------------------------------------------------------------------------- + +// flags. +#define FLAG_HIST 0x0001 +#define FLAG_VERTEX_SAMPLING 0x0002 +#define FLAG_EDGE_SAMPLING 0x0004 +#define FLAG_FACE_SAMPLING 0x0008 +#define FLAG_MONTECARLO_SAMPLING 0x0010 +#define FLAG_SUBDIVISION_SAMPLING 0x0020 +#define FLAG_SIMILAR_TRIANGLES_SAMPLING 0x0040 +#define FLAG_SAVE_ERROR_DISPLACEMENT 0x0080 +#define FLAG_SAVE_ERROR_AS_COLOUR 0x0100 + +// global constants +#define NO_SAMPLES_PER_FACE 10 +#define N_SAMPLES_EDGE_TO_FACE_RATIO 0.1 +#define BBOX_FACTOR 0.1 +#define INFLATE_PERCENTAGE 0.02 +#define MIN_SIZE 125 /* 125 = 5^3 */ +#define N_HIST_BINS 256 +#define PRINT_EVERY_N_ELEMENTS 1000 +#define FILE_EXT_SMF "smf" +#define FILE_EXT_PLY "ply" + +// ----------------------------------------------------------------------------------------------- +template +class Sampling +{ +private: + typedef GridStaticObj< typename MetroMesh::FaceContainer > MetroMeshGrid; + typedef Point3 Point3x; + + // data structures + MetroMesh &S1; + MetroMesh &S2; + MetroMeshGrid gS2; + + + // parameters + double dist_upper_bound; + double n_samples_per_area_unit; + double n_samples_target; + int Flags; + + // results +// Hist hist; + unsigned long n_total_samples; + unsigned long n_total_area_samples; + unsigned long n_total_edge_samples; + unsigned long n_total_vertex_samples; + double max_dist; + double mean_dist; + double RMS_dist; + double volume; + double area_S1; + + // globals + int n_samples; + + // private methods + inline double ComputeMeshArea(MetroMesh & mesh); + float AddSample(const Point3x &p); + inline void AddRandomSample(typename MetroMesh::FaceIterator &T); + inline void SampleEdge(const Point3x & v0, const Point3x & v1, int n_samples_per_edge); + void VertexSampling(); + void EdgeSampling(); + void FaceSubdiv(const Point3x & v0, const Point3x &v1, const Point3x & v2, int maxdepth); + void SimilarTriangles(const Point3x &v0, const Point3x &v1, const Point3x &v2, int n_samples_per_edge); + void MontecarloFaceSampling(); + void SubdivFaceSampling(); + void SimilarFaceSampling(); + +public : + // public methods + Sampling(MetroMesh &_s1, MetroMesh &_s2); + void Hausdorff(); + double GetArea() {return area_S1;} + double GetDistMax() {return max_dist;} + double GetDistMean() {return mean_dist;} + double GetDistRMS() {return RMS_dist;} + double GetDistVolume() {return volume;} + double GetNSamples() {return n_total_samples;} + double GetNAreaSamples() {return n_total_area_samples;} + double GetNEdgeSamples() {return n_total_edge_samples;} + double GetNVertexSamples() {return n_total_vertex_samples;} + double GetNSamplesPerAreaUnit() {return n_samples_per_area_unit;} + double GetNSamplesTarget() {return n_samples_target;} +// Hist &GetHist() {return hist;} + void SetFlags(int flags) {Flags = flags;} + void ClearFlag(int flag) {Flags &= (flag ^ -1);} + void SetParam(double _n_samp) {n_samples_target = _n_samp;} + void SetSamplesTarget(int _n_samp); + void SetSamplesPerAreaUnit(double _n_samp); +}; + +// ----------------------------------------------------------------------------------------------- + +// constructor +template +Sampling::Sampling(MetroMesh &_s1, MetroMesh &_s2):S1(_s1),S2(_s2) +{ + Flags = 0; + area_S1 = ComputeMeshArea(_s1); +} + + +// set sampling parameters +template +void Sampling::SetSamplesTarget(int _n_samp) +{ + n_samples_target = _n_samp; + n_samples_per_area_unit = (double) n_samples_target / area_S1; +} + +template +void Sampling::SetSamplesPerAreaUnit(double _n_samp) +{ + n_samples_per_area_unit = _n_samp; + n_samples_target = (int)((double) n_samples_per_area_unit * area_S1); +} + + +// auxiliary functions +template +inline double Sampling::ComputeMeshArea(MetroMesh & mesh) +{ + MetroMesh::FaceIterator face; + double area = 0.0; + + for(face=mesh.face.begin(); face != mesh.face.end(); face++) + if(!(*face).IsD()) + area += face->Area(); + + return area; +} + +template +float Sampling::AddSample(const Point3x &p) +{ + MetroMesh::FaceType *f=0; + Point3x normf, bestq, ip; + MetroMesh::ScalarType dist; + + dist = dist_upper_bound; + + // compute distance between p_i and the mesh S2 + MinDistPoint(S2, p, gS2, dist, normf, bestq, f, ip); + + // update distance measures + if(dist == dist_upper_bound) + return -1.0; + + if(dist > max_dist) + max_dist = dist; // L_inf + mean_dist += dist; // L_1 + RMS_dist += dist*dist; // L_2 + n_total_samples++; + + //if(Flags & FLAG_HIST) + // hist.Add((float)fabs(dist)); + + return (float)dist; +} + + +// ----------------------------------------------------------------------------------------------- +// --- Vertex Sampling --------------------------------------------------------------------------- + +template +void Sampling::VertexSampling() +{ + // Vertex sampling. + int cnt = 0; + float error; + + printf("Vertex sampling\n"); + MetroMesh::VertexIterator vi; + for(vi=S1.vert.begin();vi!=S1.vert.end();++vi) + { + error = AddSample((*vi).cP()); + n_total_vertex_samples++; + + // save vertex quality + if(Flags & (FLAG_SAVE_ERROR_DISPLACEMENT | FLAG_SAVE_ERROR_AS_COLOUR)) + (*vi).Q() = error; + +/* + if(Flags & FLAG_SAVE_ERROR_AS_COLOUR) + { + ColorUB col = ColorUB(ColorUB::White); + + if(error < dist_upper_bound) + // colour mapped distance + col.ColorRamp(0, dist_upper_bound, dist_upper_bound-error); + //else + // no matching mesh patches -> white + + (*vi).C() = col; + } +*/ + + // print progress information + if(!(++cnt % PRINT_EVERY_N_ELEMENTS)) + printf("Sampling vertices %d%%\r", (100 * cnt/S1.vn)); + } + printf(" \r"); +} + + +// ----------------------------------------------------------------------------------------------- +// --- Edge Sampling ----------------------------------------------------------------------------- + +template +inline void Sampling::SampleEdge(const Point3x & v0, const Point3x & v1, int n_samples_per_edge) +{ + // uniform sampling of the segment v0v1. + Point3x e((v1-v0)/(double)(n_samples_per_edge+1)); + int i; + + for(i=1; i <= n_samples_per_edge; i++) + { + AddSample(v0 + e*i); + n_total_edge_samples++; + } +} + + +template +void Sampling::EdgeSampling() +{ + // Edge sampling. +typedef pair pvv; + vector< pvv > Edges; + + printf("Edge sampling\n"); + + // compute edge list. + MetroMesh::FaceIterator fi; + for(fi=S1.face.begin(); fi != S1.face.end(); fi++) + for(int i=0; i<3; ++i) + { + Edges.push_back(make_pair((*fi).V0(i),(*fi).V1(i))); + if(Edges.back().first > Edges.back().second) + swap(Edges.back().first, Edges.back().second); + } + sort(Edges.begin(), Edges.end()); + vector::iterator edgeend = unique(Edges.begin(), Edges.end()); + Edges.resize(edgeend-Edges.begin()); + + // sample edges. + vector::iterator ei; + double n_samples_per_length_unit; + double n_samples_decimal = 0.0; + int cnt=0; + if(Flags & FLAG_FACE_SAMPLING) + n_samples_per_length_unit = sqrt(n_samples_per_area_unit); + else + n_samples_per_length_unit = n_samples_per_area_unit; + for(ei=Edges.begin(); ei!=Edges.end(); ++ei) + { + n_samples_decimal += Distance((*ei).first->cP(),(*ei).second->cP()) * n_samples_per_length_unit; + n_samples = (int) n_samples_decimal; + SampleEdge((*ei).first->cP(), (*ei).second->cP(), (int) n_samples); + n_samples_decimal -= (double) n_samples; + + // print progress information + if(!(++cnt % PRINT_EVERY_N_ELEMENTS)) + printf("Sampling edge %d%%\r", (100 * cnt/Edges.size())); + } + printf(" \r"); +} + + +// ----------------------------------------------------------------------------------------------- +// --- Face Sampling ----------------------------------------------------------------------------- + +// Montecarlo sampling. +template +inline void Sampling::AddRandomSample(typename MetroMesh::FaceIterator &T) +{ + // random sampling over the input face. + double rnd_1, rnd_2; + + // vertices of the face T. + Point3x p0(T->V(0)->cP()); + Point3x p1(T->V(1)->cP()); + Point3x p2(T->V(2)->cP()); + // calculate two edges of T. + Point3x v1(p1 - p0); + Point3x v2(p2 - p0); + + // choose two random numbers. + rnd_1 = (double)rand() / (double)RAND_MAX; + rnd_2 = (double)rand() / (double)RAND_MAX; + if(rnd_1 + rnd_2 > 1.0) + { + rnd_1 = 1.0 - rnd_1; + rnd_2 = 1.0 - rnd_2; + } + + // add a random point on the face T. + AddSample (p0 + (v1 * rnd_1 + v2 * rnd_2)); + n_total_area_samples++; +} + +template +void Sampling::MontecarloFaceSampling() +{ + // Montecarlo sampling. + int cnt = 0; + double n_samples_decimal = 0.0; + MetroMesh::FaceIterator fi; + + srand(clock()); + // printf("Montecarlo face sampling\n"); + for(fi=S1.face.begin(); fi != S1.face.end(); fi++) + if(!(*fi).IsD()) + { + // compute # samples in the current face. + n_samples_decimal += fi->Area() * n_samples_per_area_unit; + n_samples = (int) n_samples_decimal; + + // for every sample p_i in T... + for(int i=0; i < n_samples; i++) + AddRandomSample(fi); + + n_samples_decimal -= (double) n_samples; + + // print progress information +// if(!(++cnt % PRINT_EVERY_N_ELEMENTS)) + // printf("Sampling face %d%%\r", (100 * cnt/S1.fn)); + } + // printf(" \r"); +} + + +// Subdivision sampling. +template +void Sampling::FaceSubdiv(const Point3x & v0, const Point3x & v1, const Point3x & v2, int maxdepth) +{ + // recursive face subdivision. + if(maxdepth == 0) + { + // ground case. + AddSample((v0+v1+v2)/3.0f); + n_total_area_samples++; + n_samples++; + return; + } + + // compute the longest edge. + double maxd01 = SquaredDistance(v0,v1); + double maxd12 = SquaredDistance(v1,v2); + double maxd20 = SquaredDistance(v2,v0); + int res; + if(maxd01 > maxd12) + if(maxd01 > maxd20) res = 0; + else res = 2; + else + if(maxd12 > maxd20) res = 1; + else res = 2; + + // break the input triangle along the median to the the longest edge. + Point3x pp; + switch(res) + { + case 0 : pp = (v0+v1)/2; + FaceSubdiv(v0,pp,v2,maxdepth-1); + FaceSubdiv(pp,v1,v2,maxdepth-1); + break; + case 1 : pp = (v1+v2)/2; + FaceSubdiv(v0,v1,pp,maxdepth-1); + FaceSubdiv(v0,pp,v2,maxdepth-1); + break; + case 2 : pp = (v2+v0)/2; + FaceSubdiv(v0,v1,pp,maxdepth-1); + FaceSubdiv(pp,v1,v2,maxdepth-1); + break; + } +} + +template +void Sampling::SubdivFaceSampling() +{ + // Subdivision sampling. + int cnt = 0, maxdepth; + double n_samples_decimal = 0.0; + MetroMesh::FaceIterator fi; + + printf("Subdivision face sampling\n"); + for(fi=S1.face.begin(); fi != S1.face.end(); fi++) + { + // compute # samples in the current face. + n_samples_decimal += fi->Area() * n_samples_per_area_unit; + n_samples = (int) n_samples_decimal; + if(n_samples) + { + // face sampling. + maxdepth = (int)log((double)n_samples)/log(2.0); + n_samples = 0; + FaceSubdiv((*fi).V(0)->cP(), (*fi).V(1)->cP(), (*fi).V(2)->cP(), maxdepth); + } + n_samples_decimal -= (double) n_samples; + + // print progress information + if(!(++cnt % PRINT_EVERY_N_ELEMENTS)) + printf("Sampling face %d%%\r", (100 * cnt/S1.fn)); + } + printf(" \r"); +} + + +// Similar Triangles sampling. +template +void Sampling::SimilarTriangles(const Point3x & v0, const Point3x & v1, const Point3x & v2, int n_samples_per_edge) +{ + Point3x V1((v1-v0)/(double)(n_samples_per_edge-1)); + Point3x V2((v2-v0)/(double)(n_samples_per_edge-1)); + int i, j; + + // face sampling. + for(i=1; i < n_samples_per_edge-1; i++) + for(j=1; j < n_samples_per_edge-1-i; j++) + { + AddSample( v0 + (V1*(double)i + V2*(double)j) ); + n_total_area_samples++; + n_samples++; + } +} + +template +void Sampling::SimilarFaceSampling() +{ + // Similar Triangles sampling. + int cnt = 0, n_samples_per_edge; + double n_samples_decimal = 0.0; + MetroMesh::FaceIterator fi; + + printf("Similar Triangles face sampling\n"); + for(fi=S1.face.begin(); fi != S1.face.end(); fi++) + { + // compute # samples in the current face. + n_samples_decimal += fi->Area() * n_samples_per_area_unit; + n_samples = (int) n_samples_decimal; + if(n_samples) + { + // face sampling. + n_samples_per_edge = (int)((sqrt(1.0+8.0*(double)n_samples) +5.0)/2.0); + n_samples = 0; + SimilarTriangles((*fi).V(0)->cP(), (*fi).V(1)->cP(), (*fi).V(2)->cP(), n_samples_per_edge); + } + n_samples_decimal -= (double) n_samples; + + // print progress information + if(!(++cnt % PRINT_EVERY_N_ELEMENTS)) + printf("Sampling face %d%%\r", (100 * cnt/S1.fn)); + } + printf(" \r"); +} + + +// ----------------------------------------------------------------------------------------------- +// --- Distance ---------------------------------------------------------------------------------- + +template +void Sampling::Hausdorff() +{ + Box3< MetroMesh::ScalarType> bbox; + + // set grid meshes. + gS2.SetBBox(S2.bbox); + if(S2.face.size() < MIN_SIZE) + gS2.Set(S2.face, MIN_SIZE); + else + gS2.Set(S2.face); + + // set bounding box + bbox = S2.bbox; + dist_upper_bound = /*BBOX_FACTOR * */bbox.Diag(); + //if(Flags & FLAG_HIST) + // hist.SetRange(0.0, dist_upper_bound, N_HIST_BINS); + + // initialize sampling statistics. + n_total_area_samples = n_total_edge_samples = n_total_vertex_samples = n_total_samples = n_samples = 0; + max_dist = -HUGE_VAL; + mean_dist = RMS_dist = 0; + + // Vertex sampling. + if(Flags & FLAG_VERTEX_SAMPLING) + VertexSampling(); + // Edge sammpling. + n_samples_target -= (int) n_total_samples; + if(n_samples_target > 0) + { + n_samples_per_area_unit = n_samples_target / area_S1; + if(Flags & FLAG_EDGE_SAMPLING) + { + EdgeSampling(); + n_samples_target -= (int) n_total_samples; + } + // Face sampling. + if((Flags & FLAG_FACE_SAMPLING) && (n_samples_target > 0)) + { + n_samples_per_area_unit = n_samples_target / area_S1; + if(Flags & FLAG_MONTECARLO_SAMPLING) MontecarloFaceSampling(); + if(Flags & FLAG_SUBDIVISION_SAMPLING) SubdivFaceSampling(); + if(Flags & FLAG_SIMILAR_TRIANGLES_SAMPLING) SimilarFaceSampling(); + } + } + + // compute vertex colour + if(Flags & FLAG_SAVE_ERROR_AS_COLOUR) + { + MetroMesh::VertexIterator vi; + float error; + int cnt = 0; + for(vi=S1.vert.begin();vi!=S1.vert.end();++vi) + { + Color4b col = Color4b(Color4b::White); + error = (*vi).Q(); + + if(error < dist_upper_bound) + // colour mapped distance + col.ColorRamp(0, (float)max_dist, (float)max_dist-error); + //else + // no matching mesh patches -> white + + (*vi).C() = col; + + // print progress information + if(!(++cnt % PRINT_EVERY_N_ELEMENTS)) + printf("Computing vertex colour %d%%\r", (100 * cnt/S1.vn)); + } + printf(" \r"); + } + + // compute statistics + n_samples_per_area_unit = (double) n_total_samples / area_S1; + volume = mean_dist / n_samples_per_area_unit / 2.0; + mean_dist /= n_total_samples; + RMS_dist = sqrt(RMS_dist / n_total_samples); +} +// ----------------------------------------------------------------------------------------------- +//#undef FLAG_HIST +//#undef FLAG_VERTEX_SAMPLING +//#undef FLAG_EDGE_SAMPLING +//#undef FLAG_FACE_SAMPLING +//#undef FLAG_MONTECARLO_SAMPLING +//#undef FLAG_SUBDIVISION_SAMPLING +//#undef FLAG_SIMILAR_TRIANGLES_SAMPLING +//#undef FLAG_SAVE_ERROR_DISPLACEMENT +//#undef FLAG_SAVE_ERROR_AS_COLOUR +// +//// global constants +//#undef NO_SAMPLES_PER_FACE +//#undef N_SAMPLES_EDGE_TO_FACE_RATIO +//#undef BBOX_FACTOR +//#undef INFLATE_PERCENTAGE +//#undef MIN_SIZE +//#undef N_HIST_BINS +//#undef PRINT_EVERY_N_ELEMENTS +//#undef FILE_EXT_SMF +//#undef FILE_EXT_PLY +// + +#endif +// -----------------------------------------------------------------------------------------------