mirror
/
OpenDroneMap-ODM
kopia lustrzana https://github.com/OpenDroneMap/ODM
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

Remove all backup(~) files 

Former-commit-id: d38ceca6a9
gh-pages
Jesse Crocker 2015-01-22 18:15:33 -07:00
rodzic 94a83e8e38
commit 1d70d0e6a4
13 zmienionych plików z 1 dodań i 4700 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

1
.gitignore vendored
Wyświetl plik

@ -0,0 +1 @@
*~

848
cmvs/program/base/pmvs/patchOrganizerS.cc~
Wyświetl plik

@ -1,848 +0,0 @@
#include <string>
#include "patchOrganizerS.h"
#include "findMatch.h"
using namespace PMVS3;
using namespace Patch;
using namespace std;
Ppatch CpatchOrganizerS::m_MAXDEPTH(new Cpatch());
Ppatch CpatchOrganizerS::m_BACKGROUND(new Cpatch());
CpatchOrganizerS::CpatchOrganizerS(CfindMatch& findMatch) : m_fm(findMatch) {
}
// change the contents of m_images from images to indexes
void CpatchOrganizerS::image2index(Cpatch& patch) {
// first image has to be target image
vector<int> newimages;
for (int i = 0; i < (int)patch.m_images.size(); ++i) {
const int index = m_fm.m_pss.image2index(patch.m_images[i]);
if (index != -1)
newimages.push_back(index);
}
patch.m_images.swap(newimages);
// make sure that the reference image is the tagetting image
int exist = -1;
for (int j = 0; j < (int)patch.m_images.size(); ++j) {
if (patch.m_images[j] < m_fm.m_tnum) {
exist = j;
break;
}
}
if (exist == -1)
patch.m_images.clear();
else if (exist != 0)
swap(patch.m_images[0], patch.m_images[exist]);
}
// change the contents of m_images from indexes to images
void CpatchOrganizerS::index2image(Cpatch& patch) {
for (int i = 0; i < (int)patch.m_images.size(); ++i)
patch.m_images[i] = m_fm.m_pss.m_images[patch.m_images[i]];
for (int i = 0; i < (int)patch.m_vimages.size(); ++i)
patch.m_vimages[i] = m_fm.m_pss.m_images[patch.m_vimages[i]];
}
void CpatchOrganizerS::init(void) {
m_pgrids.clear(); m_pgrids.resize(m_fm.m_tnum);
m_vpgrids.clear(); m_vpgrids.resize(m_fm.m_tnum);
m_dpgrids.clear(); m_dpgrids.resize(m_fm.m_tnum);
m_counts.clear(); m_counts.resize(m_fm.m_tnum);
m_gwidths.clear(); m_gwidths.resize(m_fm.m_num);
m_gheights.clear(); m_gheights.resize(m_fm.m_num);
for (int index = 0; index < m_fm.m_num; ++index) {
const int gwidth = (m_fm.m_pss.getWidth(index, m_fm.m_level)
+ m_fm.m_csize - 1) / m_fm.m_csize;
const int gheight = (m_fm.m_pss.getHeight(index, m_fm.m_level)
+ m_fm.m_csize - 1) / m_fm.m_csize;
m_gwidths[index] = gwidth;
m_gheights[index] = gheight;
if (index < m_fm.m_tnum) {
m_pgrids[index].resize(gwidth * gheight);
m_vpgrids[index].resize(gwidth * gheight);
m_dpgrids[index].resize(gwidth * gheight);
m_counts[index].resize(gwidth * gheight);
fill(m_dpgrids[index].begin(), m_dpgrids[index].end(), m_MAXDEPTH);
}
}
}
//----------------------------------------------------------------------
void CpatchOrganizerS::writePatches(void) {
for (int index = 0; index < m_fm.m_tnum; ++index) {
const int image = m_fm.m_timages[index];
vector<Ppatch> ppatches;
collectNonFixPatches(index, ppatches);
if (ppatches.empty())
continue;
char buffer[1024];
sprintf(buffer, "%smodels/%08d.patc%d", m_fm.m_prefix.c_str(), image, m_fm.m_level);
ofstream ofstr;
ofstr.open(buffer);
ofstr << "PATCHES" << endl
<< (int)ppatches.size() << endl;
for (int p = 0; p < (int)ppatches.size(); ++p) {
Cpatch patch = *ppatches[p];
index2image(patch);
ofstr << patch << endl;
}
ofstr.close();
sprintf(buffer, "%smodels/%08d-%d.ply", m_fm.m_prefix.c_str(), image, m_fm.m_level);
writePLY(ppatches, buffer);
}
collectPatches(1);
char buffer[1024];
sprintf(buffer, "%smodels/m-%08d-%d.ply", m_fm.m_prefix.c_str(),
m_fm.m_pss.m_images[0], m_fm.m_level);
writePLY(m_ppatches, buffer);
}
void CpatchOrganizerS::writePatches2(const std::string prefix) {
collectPatches(1);
{
char buffer[1024];
sprintf(buffer, "%s.ply", prefix.c_str());
writePLY(m_ppatches, buffer);
}
{
char buffer[1024];
sprintf(buffer, "%s.patch", prefix.c_str());
ofstream ofstr;
ofstr.open(buffer);
ofstr << "PATCHES" << endl
<< (int)m_ppatches.size() << endl;
for (int p = 0; p < (int)m_ppatches.size(); ++p) {
Cpatch patch = *ppatches[p];
index2image(patch);
ofstr << patch << endl;
}
ofstr.close();
}
{
char buffer[1024];
sprintf(buffer, "%s.pset", prefix.c_str());
ofstream ofstr;
ofstr.open(buffer);
for (int p = 0; p < (int)m_ppatches.size(); ++p)
ofstr << m_ppatches[p]->m_coord[0] << ' '
<< m_ppatches[p]->m_coord[1] << ' '
<< m_ppatches[p]->m_coord[2] << ' '
<< m_ppatches[p]->m_normal[0] << ' '
<< m_ppatches[p]->m_normal[1] << ' '
<< m_ppatches[p]->m_normal[2] << endl;
ofstr.close();
}
}
void CpatchOrganizerS::readPatches(void) {
// Read-in existing reconstructed points. set m_fix to one for non-targetting images
for (int i = 0; i < m_fm.m_tnum; ++i) {
const int image = m_fm.m_images[i];
char buffer[1024];
sprintf(buffer, "%smodels/%08d.patc%d",
m_fm.m_prefix.c_str(), image, m_fm.m_level);
ifstream ifstr;
ifstr.open(buffer);
if (!ifstr.is_open())
continue;
string header; int pnum;
ifstr >> header >> pnum;
cerr << image << ' ' << pnum << " patches" << endl;
for (int p = 0; p < pnum; ++p) {
Ppatch ppatch(new Cpatch());
ifstr >> *ppatch;
ppatch->m_fix = 0;
ppatch->m_vimages.clear();
image2index(*ppatch);
if (ppatch->m_images.empty())
continue;
// m_vimages must be targetting images
#ifdef DEBUG
for (int j = 0; j < (int)ppatch->m_vimages.size(); ++j)
if (m_fm.m_tnum <= ppatch->m_vimages[j]) {
cerr << "Impossible in readPatches. m_vimages must be targetting images" << endl
<< "for patches stored in targetting images, if visdata2 have been consistent" << endl;
exit (1);
}
#endif
setGrids(*ppatch);
addPatch(ppatch);
}
ifstr.close();
}
// For patches in non-targetting images
for (int i = m_fm.m_tnum; i < m_fm.m_num; ++i) {
const int image = m_fm.m_images[i];
char buffer[1024];
sprintf(buffer, "%smodels/%08d.patc%d", m_fm.m_prefix.c_str(), image, m_fm.m_level);
ifstream ifstr;
ifstr.open(buffer);
if (!ifstr.is_open())
continue;
string header; int pnum;
ifstr >> header >> pnum;
cerr << image << ' ' << pnum << " patches" << endl;
for (int p = 0; p < pnum; ++p) {
Ppatch ppatch(new Cpatch());
ifstr >> *ppatch;
ppatch->m_fix = 1;
ppatch->m_vimages.clear();
image2index(*ppatch);
if (ppatch->m_images.empty())
continue;
setGrids(*ppatch);
addPatch(ppatch);
}
ifstr.close();
}
}
void CpatchOrganizerS::collectPatches(const int target) {
m_ppatches.clear();
for (int index = 0; index < m_fm.m_tnum; ++index) {
for (int i = 0; i < (int)m_pgrids[index].size(); ++i) {
vector<Ppatch>::iterator begin = m_pgrids[index][i].begin();
while (begin != m_pgrids[index][i].end()) {
(*begin)->m_id = -1;
begin++;
}
}
}
int count = 0;
for (int index = 0; index < m_fm.m_tnum; ++index) {
for (int i = 0; i < (int)m_pgrids[index].size(); ++i) {
vector<Ppatch>::iterator begin = m_pgrids[index][i].begin();
while (begin != m_pgrids[index][i].end()) {
if ((*begin)->m_id == -1) {
(*begin)->m_id = count++;
if (target == 0 || (*begin)->m_fix == 0)
m_ppatches.push_back(*begin);
}
++begin;
}
}
}
}
void CpatchOrganizerS::collectPatches(std::priority_queue<Patch::Ppatch,
std::vector<Patch::Ppatch>,
P_compare>& pqpatches) {
for (int index = 0; index < m_fm.m_tnum; ++index) {
for (int i = 0; i < (int)m_pgrids[index].size(); ++i) {
vector<Ppatch>::iterator begin = m_pgrids[index][i].begin();
while (begin != m_pgrids[index][i].end()) {
if ((*begin)->m_flag == 0) {
(*begin)->m_flag = 1;
pqpatches.push(*begin);
}
++begin;
}
}
}
}
void CpatchOrganizerS::collectPatches(const int index,
std::priority_queue<Patch::Ppatch, std::vector<Patch::Ppatch>,
P_compare>& pqpatches) {
pthread_rwlock_wrlock(&m_fm.m_imageLocks[index]);
for (int i = 0; i < (int)m_pgrids[index].size(); ++i) {
vector<Ppatch>::iterator begin = m_pgrids[index][i].begin();
vector<Ppatch>::iterator end = m_pgrids[index][i].end();
while (begin != end) {
if ((*begin)->m_images[0] == index && (*begin)->m_flag == 0) {
(*begin)->m_flag = 1;
pqpatches.push(*begin);
}
++begin;
}
}
pthread_rwlock_unlock(&m_fm.m_imageLocks[index]);
}
// Should be used only for writing
void CpatchOrganizerS::collectNonFixPatches(const int index,
std::vector<Patch::Ppatch>& ppatches) {
pthread_rwlock_wrlock(&m_fm.m_imageLocks[index]);
for (int i = 0; i < (int)m_pgrids[index].size(); ++i) {
vector<Ppatch>::iterator begin = m_pgrids[index][i].begin();
vector<Ppatch>::iterator end = m_pgrids[index][i].end();
while (begin != end) {
if ((*begin)->m_images[0] == index && (*begin)->m_fix == 0) {
ppatches.push_back(*begin);
}
++begin;
}
}
pthread_rwlock_unlock(&m_fm.m_imageLocks[index]);
}
void CpatchOrganizerS::clearFlags(void) {
vector<Ppatch>::iterator bppatch = m_ppatches.begin();
vector<Ppatch>::iterator eppatch = m_ppatches.end();
while (bppatch != eppatch) {
(*bppatch)->m_flag = 0;
++bppatch;
}
}
void CpatchOrganizerS::clearCounts(void) {
for (int index = 0; index < m_fm.m_tnum; ++index) {
vector<unsigned char>::iterator begin = m_counts[index].begin();
vector<unsigned char>::iterator end = m_counts[index].end();
while (begin != end) {
*begin = (unsigned char)0;
++begin;
}
}
}
void CpatchOrganizerS::addPatch(Patch::Ppatch& ppatch) {
// First handle m_vimages
vector<int>::iterator bimage = ppatch->m_images.begin();
vector<int>::iterator eimage = ppatch->m_images.end();
vector<Vec2i>::iterator bgrid = ppatch->m_grids.begin();
while (bimage != eimage) {
const int index = *bimage;
if (m_fm.m_tnum <= index) {
++bimage; ++bgrid;
continue;
}
const int index2 = (*bgrid)[1] * m_gwidths[index] + (*bgrid)[0];
pthread_rwlock_wrlock(&m_fm.m_imageLocks[index]);
m_pgrids[index][index2].push_back(ppatch);
pthread_rwlock_unlock(&m_fm.m_imageLocks[index]);
++bimage;
++bgrid;
}
// If depth, set vimages
if (m_fm.m_depth == 0)
return;
bimage = ppatch->m_vimages.begin();
eimage = ppatch->m_vimages.end();
bgrid = ppatch->m_vgrids.begin();
while (bimage != eimage) {
const int index = *bimage;
const int index2 = (*bgrid)[1] * m_gwidths[index] + (*bgrid)[0];
pthread_rwlock_wrlock(&m_fm.m_imageLocks[index]);
m_vpgrids[index][index2].push_back(ppatch);
pthread_rwlock_unlock(&m_fm.m_imageLocks[index]);
++bimage;
++bgrid;
}
updateDepthMaps(ppatch);
}
void CpatchOrganizerS::updateDepthMaps(Ppatch& ppatch) {
for (int image = 0; image < m_fm.m_tnum; ++image) {
const Vec3f icoord = m_fm.m_pss.project(image, ppatch->m_coord, m_fm.m_level);
const float fx = icoord[0] / m_fm.m_csize;
const int xs[2] = {(int)floor(fx), (int)ceil(fx)};
const float fy = icoord[1] / m_fm.m_csize;
const int ys[2] = {(int)floor(fy), (int)ceil(fy)};
const float depth = m_fm.m_pss.m_photos[image].m_oaxis * ppatch->m_coord;
pthread_rwlock_wrlock(&m_fm.m_imageLocks[image]);
for (int j = 0; j < 2; ++j) {
for (int i = 0; i < 2; ++i) {
if (xs[i] < 0 || m_gwidths[image] <= xs[i] ||
ys[j] < 0 || m_gheights[image] <= ys[j])
continue;
const int index = ys[j] * m_gwidths[image] + xs[i];
if (m_dpgrids[image][index] == m_MAXDEPTH)
m_dpgrids[image][index] = ppatch;
else {
const float dtmp = m_fm.m_pss.m_photos[image].m_oaxis *
m_dpgrids[image][index]->m_coord;
if (depth < dtmp)
m_dpgrids[image][index] = ppatch;
}
}
}
pthread_rwlock_unlock(&m_fm.m_imageLocks[image]);
}
}
void CpatchOrganizerS::setGridsImages(Patch::Cpatch& patch,
const std::vector<int>& images) const {
patch.m_images.clear();
patch.m_grids.clear();
vector<int>::const_iterator bimage = images.begin();
vector<int>::const_iterator eimage = images.end();
while (bimage != eimage) {
const Vec3f icoord = m_fm.m_pss.project(*bimage, patch.m_coord, m_fm.m_level);
const int ix = ((int)floor(icoord[0] + 0.5f)) / m_fm.m_csize;
const int iy = ((int)floor(icoord[1] + 0.5f)) / m_fm.m_csize;
if (0 <= ix && ix < m_gwidths[*bimage] &&
0 <= iy && iy < m_gheights[*bimage]) {
patch.m_images.push_back(*bimage);
patch.m_grids.push_back(Vec2i(ix, iy));
}
++bimage;
}
}
void CpatchOrganizerS::setGrids(Ppatch& ppatch) const{
setGrids(*ppatch);
}
void CpatchOrganizerS::setGrids(Cpatch& patch) const{
patch.m_grids.clear();
for (int i = 0; i < (int)patch.m_images.size(); ++i) {
const int image = patch.m_images[i];
Vec3f icoord = m_fm.m_pss.project(image, patch.m_coord, m_fm.m_level);
const int ix = ((int)floor(icoord[0] + 0.5f)) / m_fm.m_csize;
const int iy = ((int)floor(icoord[1] + 0.5f)) / m_fm.m_csize;
patch.m_grids.push_back(TVec2<int>(ix, iy));
}
}
void CpatchOrganizerS::setVImagesVGrids(Ppatch& ppatch) {
setVImagesVGrids(*ppatch);
}
void CpatchOrganizerS::setVImagesVGrids(Cpatch& patch) {
vector<int> used;
used.resize(m_fm.m_tnum);
fill(used.begin(), used.end(), 0);
vector<int>::iterator bimage = patch.m_images.begin();
vector<int>::iterator eimage = patch.m_images.end();
while (bimage != eimage) {
if ((*bimage) < m_fm.m_tnum)
used[*(bimage)] = 1;
++bimage;
}
bimage = patch.m_vimages.begin();
eimage = patch.m_vimages.end();
while (bimage != eimage)
used[*(bimage++)] = 1;
for (int image = 0; image < m_fm.m_tnum; ++image) {
if (used[image])
continue;
int ix, iy;
if (isVisible0(patch, image, ix, iy,
m_fm.m_neighborThreshold, 1) == 0) {
continue;
}
if (m_fm.m_pss.getEdge(patch.m_coord, image, m_fm.m_level) == 0)
continue;
patch.m_vimages.push_back(image);
patch.m_vgrids.push_back(TVec2<int>(ix, iy));
}
}
void CpatchOrganizerS::removePatch(const Ppatch& ppatch) {
for (int i = 0; i < (int)ppatch->m_images.size(); ++i) {
const int image = ppatch->m_images[i];
if (m_fm.m_tnum <= image)
continue;
const int& ix = ppatch->m_grids[i][0];
const int& iy = ppatch->m_grids[i][1];
const int index = iy * m_gwidths[image] + ix;
m_pgrids[image][index].erase(remove(m_pgrids[image][index].begin(),
m_pgrids[image][index].end(),
ppatch),
m_pgrids[image][index].end());
}
for (int i = 0; i < (int)ppatch->m_vimages.size(); ++i) {
const int image = ppatch->m_vimages[i];
#ifdef DEBUG
if (m_fm.m_tnum <= image) {
cerr << "Impossible in removePatch. m_vimages must be targetting images" << endl;
exit (1);
}
#endif
const int& ix = ppatch->m_vgrids[i][0];
const int& iy = ppatch->m_vgrids[i][1];
const int index = iy * m_gwidths[image] + ix;
m_vpgrids[image][index].erase(remove(m_vpgrids[image][index].begin(),
m_vpgrids[image][index].end(),
ppatch),
m_vpgrids[image][index].end());
}
}
int CpatchOrganizerS::isVisible0(const Cpatch& patch, const int image,
int& ix, int& iy,
const float strict, const int lock) {
const Vec3f icoord =
m_fm.m_pss.project(image, patch.m_coord, m_fm.m_level);
ix = ((int)floor(icoord[0] + 0.5f)) / m_fm.m_csize;
iy = ((int)floor(icoord[1] + 0.5f)) / m_fm.m_csize;
return isVisible(patch, image, ix, iy, strict, lock);
}
int CpatchOrganizerS::isVisible(const Cpatch& patch, const int image,
const int& ix, const int& iy,
const float strict, const int lock) {
const int& gwidth = m_gwidths[image];
const int& gheight = m_gheights[image];
if (ix < 0 || gwidth <= ix || iy < 0 || gheight <= iy)
return 0;
if (m_fm.m_depth == 0)
return 1;
int ans = 0;
Ppatch dppatch = m_MAXDEPTH;
const int index = iy * gwidth + ix;
if (lock)
pthread_rwlock_rdlock(&m_fm.m_imageLocks[image]);
if (m_dpgrids[image][index] == m_MAXDEPTH)
ans = 1;
else
dppatch = m_dpgrids[image][index];
if (lock)
pthread_rwlock_unlock(&m_fm.m_imageLocks[image]);
if (ans == 1)
return 1;
Vec4f ray = patch.m_coord - m_fm.m_pss.m_photos[image].m_center;
unitize(ray);
const float diff = ray * (patch.m_coord - dppatch->m_coord);
const float factor = min(2.0, 2.0 + ray * patch.m_normal);
if (diff < m_fm.m_optim.getUnit(image, patch.m_coord) * m_fm.m_csize * strict * factor)
return 1;
else
return 0;
}
void CpatchOrganizerS::findNeighbors(const Patch::Cpatch& patch,
std::vector<Patch::Ppatch>& neighbors,
const int lock, const float scale,
const int margin,
const int skipvis) {
const float radius = 1.5 * margin * m_fm.m_expand.computeRadius(patch);
vector<int>::const_iterator bimage = patch.m_images.begin();
vector<int>::const_iterator eimage = patch.m_images.end();
vector<Vec2i>::const_iterator bgrid = patch.m_grids.begin();
#ifdef DEBUG
if (patch.m_images.empty()) {
cerr << "Empty patches in findCloses" << endl;
exit (1);
}
#endif
float unit = 0.0f;
for (int i = 0; i < (int)patch.m_images.size(); ++i)
unit += m_fm.m_optim.getUnit(patch.m_images[i], patch.m_coord);
unit /= (int)patch.m_images.size();
unit *= m_fm.m_csize;
while (bimage != eimage) {
if (m_fm.m_tnum <= *bimage) {
++bimage;
++bgrid;
continue;
}
const int image = *bimage;
const int& ix = (*bgrid)[0];
const int& iy = (*bgrid)[1];
if (lock)
pthread_rwlock_rdlock(&m_fm.m_imageLocks[image]);
for (int j = -margin; j <= margin; ++j) {
const int ytmp = iy + j;
if (ytmp < 0 || m_fm.m_pos.m_gheights[image] <= ytmp)
continue;
for (int i = -margin; i <= margin; ++i) {
const int xtmp = ix + i;
if (xtmp < 0 || m_fm.m_pos.m_gwidths[image] <= xtmp)
continue;
const int index = ytmp * m_fm.m_pos.m_gwidths[image] + xtmp;
vector<Ppatch>::const_iterator bpatch =
m_fm.m_pos.m_pgrids[image][index].begin();
vector<Ppatch>::const_iterator epatch =
m_fm.m_pos.m_pgrids[image][index].end();
while (bpatch != epatch) {
if (m_fm.isNeighborRadius(patch, **bpatch, unit,
m_fm.m_neighborThreshold * scale,
radius))
neighbors.push_back(*bpatch);
++bpatch;
}
bpatch = m_fm.m_pos.m_vpgrids[image][index].begin();
epatch = m_fm.m_pos.m_vpgrids[image][index].end();
while (bpatch != epatch) {
if (m_fm.isNeighborRadius(patch, **bpatch, unit,
m_fm.m_neighborThreshold * scale,
radius))
neighbors.push_back(*bpatch);
++bpatch;
}
}
}
if (lock)
pthread_rwlock_unlock(&m_fm.m_imageLocks[image]);
++bimage;
++bgrid;
}
if (skipvis == 0) {
bimage = patch.m_vimages.begin();
eimage = patch.m_vimages.end();
bgrid = patch.m_vgrids.begin();
while (bimage != eimage) {
const int image = *bimage;
const int& ix = (*bgrid)[0];
const int& iy = (*bgrid)[1];
if (lock)
pthread_rwlock_rdlock(&m_fm.m_imageLocks[image]);
for (int j = -margin; j <= margin; ++j) {
const int ytmp = iy + j;
if (ytmp < 0 || m_fm.m_pos.m_gheights[image] <= ytmp)
continue;
for (int i = -margin; i <= margin; ++i) {
const int xtmp = ix + i;
if (xtmp < 0 || m_fm.m_pos.m_gwidths[image] <= xtmp)
continue;
const int index = ytmp * m_fm.m_pos.m_gwidths[image] + xtmp;
vector<Ppatch>::const_iterator bpatch =
m_fm.m_pos.m_pgrids[image][index].begin();
vector<Ppatch>::const_iterator epatch =
m_fm.m_pos.m_pgrids[image][index].end();
while (bpatch != epatch) {
if (m_fm.isNeighborRadius(patch, **bpatch, unit,
m_fm.m_neighborThreshold * scale,
radius))
neighbors.push_back(*bpatch);
++bpatch;
}
bpatch = m_fm.m_pos.m_vpgrids[image][index].begin();
epatch = m_fm.m_pos.m_vpgrids[image][index].end();
while (bpatch != epatch) {
if (m_fm.isNeighborRadius(patch, **bpatch, unit,
m_fm.m_neighborThreshold * scale,
radius))
neighbors.push_back(*bpatch);
++bpatch;
}
}
}
if (lock)
pthread_rwlock_unlock(&m_fm.m_imageLocks[image]);
++bimage;
++bgrid;
}
}
sort(neighbors.begin(), neighbors.end());
neighbors.erase(unique(neighbors.begin(), neighbors.end()), neighbors.end());
}
float CpatchOrganizerS::computeUnit(const Patch::Cpatch& patch) const{
float unit = 0.0f;
for (int i = 0; i < (int)patch.m_images.size(); ++i)
unit += m_fm.m_optim.getUnit(patch.m_images[i], patch.m_coord);
unit /= (int)patch.m_images.size();
unit *= m_fm.m_csize;
return unit;
}
void CpatchOrganizerS::setScales(Patch::Cpatch& patch) const {
const float unit = m_fm.m_optim.getUnit(patch.m_images[0], patch.m_coord);
const float unit2 = 2.0f * unit;
Vec4f ray = patch.m_coord - m_fm.m_pss.m_photos[patch.m_images[0]].m_center;
unitize(ray);
const int inum = min(m_fm.m_tau, (int)patch.m_images.size());
// First compute, how many pixel difference per unit along vertical
//for (int i = 1; i < (int)patch.m_images.size(); ++i) {
for (int i = 1; i < inum; ++i) {
Vec3f diff = m_fm.m_pss.project(patch.m_images[i], patch.m_coord, m_fm.m_level) -
m_fm.m_pss.project(patch.m_images[i], patch.m_coord - ray * unit2, m_fm.m_level);
patch.m_dscale += norm(diff);
}
// set m_dscale to the vertical distance where average pixel move is half pixel
//patch.m_dscale /= (int)patch.m_images.size() - 1;
patch.m_dscale /= inum - 1;
patch.m_dscale = unit2 / patch.m_dscale;
patch.m_ascale = atan(patch.m_dscale / (unit * m_fm.m_wsize / 2.0f));
}
// write out results
void CpatchOrganizerS::writePLY(const std::vector<Ppatch>& patches,
const std::string filename) {
ofstream ofstr;
ofstr.open(filename.c_str());
ofstr << "ply" << endl
<< "format ascii 1.0" << endl
<< "element vertex " << (int)patches.size() << endl
<< "property float x" << endl
<< "property float y" << endl
<< "property float z" << endl
<< "property float nx" << endl
<< "property float ny" << endl
<< "property float nz" << endl
<< "property uchar diffuse_red" << endl
<< "property uchar diffuse_green" << endl
<< "property uchar diffuse_blue" << endl
<< "end_header" << endl;
vector<Ppatch>::const_iterator bpatch = patches.begin();
vector<Ppatch>::const_iterator bend = patches.end();
while (bpatch != bend) {
// Get color
Vec3i color;
const int mode = 0;
// 0: color from images
// 1: fix
// 2: angle
if (mode == 0) {
int denom = 0;
Vec3f colorf;
for (int i = 0; i < (int)(*bpatch)->m_images.size(); ++i) {
const int image = (*bpatch)->m_images[i];
colorf += m_fm.m_pss.getColor((*bpatch)->m_coord, image, m_fm.m_level);
denom++;
}
colorf /= denom;
color[0] = min(255,(int)floor(colorf[0] + 0.5f));
color[1] = min(255,(int)floor(colorf[1] + 0.5f));
color[2] = min(255,(int)floor(colorf[2] + 0.5f));
}
else if (mode == 1) {
if ((*bpatch)->m_tmp == 1.0f) {
color[0] = 255;
color[1] = 0;
color[2] = 0;
}
else {
color[0] = 255;
color[1] = 255;
color[2] = 255;
}
}
else if (mode == 2) {
float angle = 0.0f;
vector<int>::iterator bimage = (*bpatch)->m_images.begin();
vector<int>::iterator eimage = (*bpatch)->m_images.end();
while (bimage != eimage) {
const int index = *bimage;
Vec4f ray = m_fm.m_pss.m_photos[index].m_center - (*bpatch)->m_coord;
ray[3] = 0.0f;
unitize(ray);
angle += acos(ray * (*bpatch)->m_normal);
++bimage;
}
angle = angle / (M_PI / 2.0f);
float r, g, b;
Image::Cimage::gray2rgb(angle, r, g, b);
color[0] = (int)(r * 255.0f);
color[1] = (int)(g * 255.0f);
color[2] = (int)(b * 255.0f);
}
ofstr << (*bpatch)->m_coord[0] << ' '
<< (*bpatch)->m_coord[1] << ' '
<< (*bpatch)->m_coord[2] << ' '
<< (*bpatch)->m_normal[0] << ' '
<< (*bpatch)->m_normal[1] << ' '
<< (*bpatch)->m_normal[2] << ' '
<< color[0] << ' ' << color[1] << ' ' << color[2] << endl;
++bpatch;
}
ofstr.close();
}
void CpatchOrganizerS::writePLY(const std::vector<Ppatch>& patches,
const std::string filename,
const std::vector<Vec3i>& colors) {
ofstream ofstr;
ofstr.open(filename.c_str());
ofstr << "ply" << endl
<< "format ascii 1.0" << endl
<< "element vertex " << (int)patches.size() << endl
<< "property float x" << endl
<< "property float y" << endl
<< "property float z" << endl
<< "property float nx" << endl
<< "property float ny" << endl
<< "property float nz" << endl
<< "property uchar diffuse_red" << endl
<< "property uchar diffuse_green" << endl
<< "property uchar diffuse_blue" << endl
<< "end_header" << endl;
vector<Ppatch>::const_iterator bpatch = patches.begin();
vector<Ppatch>::const_iterator bend = patches.end();
vector<Vec3i>::const_iterator colorb = colors.begin();
while (bpatch != bend) {
ofstr << (*bpatch)->m_coord[0] << ' '
<< (*bpatch)->m_coord[1] << ' '
<< (*bpatch)->m_coord[2] << ' '
<< (*bpatch)->m_normal[0] << ' '
<< (*bpatch)->m_normal[1] << ' '
<< (*bpatch)->m_normal[2] << ' '
<< *colorb << endl;
++bpatch;
++colorb;
}
ofstr.close();
}

50
cmvs/program/base/stann/assert.hpp~
Wyświetl plik

@ -1,50 +0,0 @@
/*****************************************************************************/
/* */
/* Header: assert.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
#ifndef REVIVER_ASSERT_HPP
#define REVIVER_ASSERT_HPP
#include <iostream>
/*!
\file
\brief Implementation of assert functions
This file contains an assert function implementation
*/
/*!
\brief Assert function
This function prints error messages along with file and line numbers
\param b bool to be checked
\param desc Error message
\param line Line number
\param file File name
\return True is there is no error
*/
bool MyAssertFunction( bool b, char* desc, int line, char* file){
if (b) return true;
std::cerr << "\n\nAssertion Failure\n";
std::cerr << "Description : " << desc << std::endl;
std::cerr << "Filename : " << file << std::endl;
std::cerr << "Line No : " << line << std::endl;
exit(1);
}
#if defined( _DEBUG )
#define Assert( exp, description ) MyAssertFunction( (int)(exp), description, __LINE__, __FILE__ )
#else
#define Assert( exp, description )
#endif
#endif

126
cmvs/program/base/stann/bruteNN.hpp~
Wyświetl plik

@ -1,126 +0,0 @@
/*****************************************************************************/
/* */
/* Header: bruteNN.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
/*!
\file
\brief Brute force nearest neighbor search implementation
This file contains a simple brute force NN algorithm implementation.
It is very innefficient and should be used only for accuracy tests.
*/
#include <vector>
#include <qknn.hpp>
using namespace std;
/*!
\brief A Brute force NN search class
This class is a simple brute force algorithm for computing nearest neighbors.
It is very innefficient and should be used only for accuracy tests.
*/
template<typename Point>
class bruteNN
{
public:
/*!
\brief Constructor
Constructs a brute force nearest neighbor search structure out
of the given data points.
\param *points Pointer to list of input points
\param N number of input points
*/
bruteNN(Point *points, long int N);
/*!
\brief Destructor
*/
~bruteNN();
/*!
\brief Nearest neighbor search function
This function returns the nearest k neighbors from the data
structure to the given query point. The return vector contains
the indexes to the answer points in the original data array passed
at construction time.
\param q Query point
\param k number of neighbors to search for
\param nn_idx Vector in which answer is written
*/
void ksearch(Point q, int k,
vector<long unsigned int> &nn_idx);
/*!
\brief Nearest neighbor search function
This function returns the nearest k neighbors from the data
structure to the given query point, as well as the squared
distance to the query point. The return vector contains
the indexes to the answer points in the original data array
passed at construction time.
\param q Query point
\param k number of neighbors to search for
\param nn_idx Vector in which answer is written
\param d_index Vector in which square distances are written
*/
void ksearch(Point q, int k,
vector<long unsigned int> &nn_idx,
vector<double> &d_index);
private:
vector<Point> points;
};
template<typename Point>
bruteNN<Point>::bruteNN(Point *p, long int N)
{
points.resize(N);
for(int i=0;i < N;++i)
{
points[i] = p[i];
}
}
template<typename Point>
bruteNN<Point>::~bruteNN()
{
}
template<typename Point>
void bruteNN<Point>::ksearch(Point q, int k,
vector<long unsigned int> &nn_idx)
{
qknn que;
double distance;
que.set_size(k);
for(unsigned int i=0;i < points.size();++i)
{
distance = q.sqr_dist(points[i]);
que.update(distance, i);
}
que.answer(nn_idx);
}
template<typename Point>
void bruteNN<Point>::ksearch(Point q, int k,
vector<long unsigned int> &nn_idx, vector<double> &d_indx)
{
qknn que;
double distance;
for(int i=0;i < points.size();++i)
{
distance = q.sqr_dist(points[i]);
que.update(i, distance);
}
que.answer(nn_idx, d_indx);
}

149
cmvs/program/base/stann/bsearch.hpp~
Wyświetl plik

@ -1,149 +0,0 @@
/*****************************************************************************/
/* */
/* Header: bsearch.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
#ifndef __SFCNN_BSEARCH__
#define __SFCNN_BSEARCH__
#include <vector>
#include <zorder_lt.hpp>
/*! \file
\brief Binary search functions
This file contains binary search functions for z-order operations.
*/
//! Binary search function
/*!
This function executes a binary search on a vector of pointers to points
\param A Vector of pointers to search
\param *q Pointer to query point
\param lt A less_than comparetor
\return If found: index of point. Otherwise: index of first smaller point
*/
template<typename Point>
long int BinarySearch(vector<Point *> &A, Point *q, zorder_lt<Point> lt)
{
long int low = 0;
long int high = A.size()-1;
long int middle;
while(low <= high)
{
middle = (low+high)/2;
if(q == A[middle])
return middle;
else if(lt(q, A[middle]))
high = middle-1;
else
low = middle+1;
}
return middle;
}
//! A Binary Search function
/*!
This function conducts a binary search for two points at the same time.
\param A Reference to the vector of points being searched
\param q1 pointer to first point to be searched for
\param q2 pointer to second point to be searched for
\param lt less than comparetor
\param p1 reference to return value for q1 location
\param p2 reference to return value for q2 location
*/
template<typename Point>
void PairBinarySearch(vector<Point> &A, Point q1,
Point q2, zorder_lt<Point> lt, int &p1, int &p2)
{
int low_q1=0;
int low_q2=0;
int high_q1 = A.size()-1;
int high_q2 = A.size()-1;
int middle = 0;
int middle_store;
p1 = -2;
p2 = -2;
while((low_q1 == low_q2) && (high_q1 == high_q2) && (p1 == -2) && (p2 == -2))
{
middle = (low_q1+high_q1)/2;
if(q1 == A[middle])
p1 = middle;
else if(lt(q1, A[middle]))
high_q1 = middle-1;
else
low_q1 = middle+1;
if(q2 == A[middle])
p2 = middle;
else if(lt(q2, A[middle]))
high_q2 = middle-1;
else
low_q2 = middle+1;
}
middle_store = middle;
while(low_q1 <= high_q1)
{
middle = (low_q1+high_q1)/2;
if(q1 == A[middle])
break;
else if(lt(q1, A[middle]))
high_q1 = middle-1;
else
low_q1 = middle+1;
}
p1 = middle;
middle = middle_store;
while(low_q2 <= high_q2)
{
middle = (low_q2+high_q2)/2;
if(q2 == A[middle])
break;
else if(lt(q2, A[middle]))
high_q2 = middle-1;
else
low_q2 = middle+1;
}
p2 = middle;
}
//! A binary search Function.
/*
This function executes a binary search on a vector of points
\param A Vector of points to search
\param q Query point
\param lt A less_than comparetor
\return If found: index of point. Otherwise: index of first smaller point
*/
template<typename Point>
long int BinarySearch(vector<Point> &A, Point q, zorder_lt<Point> lt)
{
long int low = 0;
long int high = A.size()-1;
long int middle = 0;
while(low <= high)
{
middle = (low+high)/2;
if(q == A[middle])
return middle;
else if(lt(q, A[middle]))
high = middle-1;
else
low = middle+1;
}
return middle;
}
#endif

498
cmvs/program/base/stann/sfcnn.hpp~
Wyświetl plik

@ -1,498 +0,0 @@
/*****************************************************************************/
/* */
/* Header: sfcnn.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
#ifndef __SFCNN___
#define __SFCNN___
#include <cstdlib>
#include <cmath>
#include <climits>
#include <vector>
#include <queue>
#include <algorithm>
#include <pthread.h>
#include <pair_iter.hpp>
#include <qknn.hpp>
#include <zorder_lt.hpp>
#include <bsearch.hpp>
/*!
\mainpage STANN Doxygen Index Page
This Doxygen API is intended for users interested in modifying
the STANN library. For instructions on using STANN, see the
STANN webpage at http://www.compgeom.com/~stann or the README.txt
file located in the main directory of the STANN distribution.
The API is still being updated. For requests/comments/errors,
please contact us at stann @ compgeom dot com
*/
/*!
\file
\brief Space filling curve nearest neighbor search
This file contains the implementation of a space filling curve
nearest neighbor search data structure
*/
using namespace std;
/*!
\brief A space filling curve nearest neighbor class.
This is the workhorse class for the sfcnn class.
The Space Filling Curve Nearest Neighbor (SFCNN) algorithm sorts the
input data set into 2-order Morton ordering. Nearest neighbors are then
calculated based on that curve. The algorithm has a runtime of O(ln(N)),
a construction time of O(Nlog(N)), and a space requirement of
O(N). The query functions of the algorithm are thread-safe.
*/
template <typename Point>
class sfcnn_work
{
template<typename P, unsigned int D, typename N> friend class sfcnn;
public:
/*!
\brief Constructor
*/
sfcnn_work();
/*!
\brief Destructor
*/
~sfcnn_work();
private:
/*!
\brief Nearest Neighbor search function
Searches for the k nearest neighbors to the point q. The answer
vector returned will contain the indexes to the answer points
This function is thread-safe.
\param q The query point
\param k The number of neighbors to return
\param nn_idx Answer vector
\param eps Error tolerence, default of 0.0.
\param sl Unused, for backwards compatibility
*/
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, float eps = 0);
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, vector<double> &dist, float eps = 0);
vector<Point> points;
vector<long unsigned int> pointers;
zorder_lt<Point> lt;
float eps;
typename Point::__NumType max, min;
void compute_bounding_box(Point q, Point &q1, Point &q2, double r);
void sfcnn_work_init(int num_threads);
void ksearch_common(Point q, unsigned int k, long unsigned int j, qknn &que, float Eps);
inline void recurse(int s, int n, Point q,
qknn &ans, Point &q1, Point &q2, int bb, int bt);
};
template<typename Point>
sfcnn_work<Point>::sfcnn_work()
{
}
template<typename Point>
void sfcnn_work<Point>::sfcnn_work_init(int num_threads)
{
max = numeric_limits<typename Point::__NumType>::max();
min = numeric_limits<typename Point::__NumType>::min();
pair_iter<typename vector<Point>::iterator,
typename vector<long unsigned int>::iterator> a(points.begin(), pointers.begin()),
b(points.end(), pointers.end());
sort(a,b,lt);
}
template<typename Point>
sfcnn_work<Point>::~sfcnn_work()
{
}
template<typename Point>
void sfcnn_work<Point>::recurse(int s, // Starting index
int n, // Number of points
Point q, // Query point
qknn &ans, // Answer que
Point &bound_box_lower_corner,
Point &bound_box_upper_corner,
int initial_scan_lower_range,
int initial_scan_upper_range)
{
if(n < 4)
{
if(n == 0) return;
bool update=false;
for(int i=0;i < n;++i)
{
if((s+i >= initial_scan_lower_range)
&& (s+i < initial_scan_upper_range))
continue;
update = ans.update(points[s+i].sqr_dist(q), pointers[s+i]) || update;
}
if(update)
compute_bounding_box(q, bound_box_lower_corner, bound_box_upper_corner, sqrt(ans.topdist()));
return;
}
if((s+n/2 >= initial_scan_lower_range) && (s+n/2 < initial_scan_upper_range))
{
}
else if(ans.update(points[s+n/2].sqr_dist(q), pointers[s+n/2]))
compute_bounding_box(q, bound_box_lower_corner, bound_box_upper_corner, sqrt(ans.topdist()));
double dsqb = lt.dist_sq_to_quad_box(q,points[s], points[s+n-1]);
//cout << "dsqb: " << dsqb << endl;
//cout << "dist: " << ans.topdist() << endl;
//cout << "p1 : " << points[s] << endl;
//cout << "p2 : " << points[s+n-1] << endl;
if(dsqb > ans.topdist())
return;
if(lt(q,points[s+n/2]))
{
//search_queue.push(pair<int, int>(s, n/2));
recurse(s, n/2, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range);
if(lt(points[s+n/2],bound_box_upper_corner))
//search_queue.push(pair<int, int>(s+n/2+1, n-n/2-1));
recurse(s+n/2+1,n-n/2-1, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range);
}
else
{
recurse(s+n/2+1, n-n/2-1, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range);
//search_queue.push(pair<int, int>(s+n/2+1, n-n/2-1));
if(lt(bound_box_lower_corner,points[s+n/2]))
//search_queue.push(pair<int, int>(s, n/2));
recurse(s, n/2, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range);
}
}
template<typename Point>
void sfcnn_work<Point>::compute_bounding_box(Point q, Point &q1, Point &q2, double R)
{
typename Point::__NumType radius;
radius = (typename Point::__NumType) ceil(R);
for(unsigned int i=0;i<Point::__DIM;++i)
{
if(q[i] < (min+radius)) q1[i] = min;
else q1[i] = q[i]-radius;
if(q[i] > (max-radius)) q2[i] = max;
else q2[i] = q[i]+radius;
}
}
template<typename Point>
void sfcnn_work<Point>::ksearch_common(Point q, unsigned int k, long unsigned int query_point_index, qknn &que, float Eps)
{
Point bound_box_lower_corner, bound_box_upper_corner;
Point low, high;
que.set_size(k);
eps=1.0+Eps;
if(query_point_index >= (k)) query_point_index -= (k);
else query_point_index=0;
int initial_scan_upper_range=query_point_index+2*k+1;
if(initial_scan_upper_range > (int)points.size())
initial_scan_upper_range = points.size();
low = points[query_point_index];
high = points[initial_scan_upper_range-1];
for(int i=query_point_index;i<initial_scan_upper_range;++i)
{
que.update(points[i].sqr_dist(q), pointers[i]);
}
compute_bounding_box(q, bound_box_lower_corner, bound_box_upper_corner, sqrt(que.topdist()));
if(lt(bound_box_upper_corner, high) && lt(low,bound_box_lower_corner))
{
//cout << "Inital search break!" << endl;
//cout << "Bb1: " << bound_box_lower_corner << endl;
//cout << "Bb2: " << bound_box_upper_corner << endl;
return;
}
//Recurse through the entire set
recurse(0, points.size(), q, que,
bound_box_lower_corner,
bound_box_upper_corner,
query_point_index,
initial_scan_upper_range);
//cout << "Bb1: " << bound_box_lower_corner << endl;
//cout << "Bb2: " << bound_box_upper_corner << endl;
}
template<typename Point>
void sfcnn_work<Point>::ksearch(Point q, unsigned int k,
vector<long unsigned int> &nn_idx, float Eps)
{
long unsigned int query_point_index;
qknn que;
query_point_index = BinarySearch(points, q, lt);
ksearch_common(q, k, query_point_index, que, Eps);
que.answer(nn_idx);
}
template<typename Point>
void sfcnn_work<Point>::ksearch(Point q, unsigned int k,
vector<long unsigned int> &nn_idx, vector<double> &dist, float Eps)
{
long unsigned int query_point_index;
qknn que;
query_point_index = BinarySearch(points, q, lt);
ksearch_common(q, k, query_point_index, que, Eps);
que.answer(nn_idx, dist);
}
/*!
\brief Nearest Neighbor search class
This class is a wrapper class for the sfcnn_work class.
There are several template specializations for this class
that select the appropriate seperated float type to be used
for various floating point coordinate types.
\param Point Data type that stores user points
\param Dim Dimension of user points
\param NumType Data type that stores coordinates of user points
*/
template <typename Point, unsigned int Dim, typename NumType>
class sfcnn
{
public:
/*!
\brief Default Constructor
*/
sfcnn(){};
/*!
\brief Constructor
Constructs a nearest neighbor data structure using the
given data points.
\param *points Pointer to the first data point
\param N number of data points
\param num_threads Currently unused. (multiple processor construction soon!)
*/
sfcnn(Point *points, long int N, int num_threads=1)
{
sfcnn_init(points,N,num_threads);
};
~sfcnn(){};
/*!
\brief Nearest Neighbor search function
Searches for the k nearest neighbors to the point q. The answer
vector returned will contain the indexes to the answer points
This function is thread-safe.
\param q The query point
\param k The number of neighbors to return
\param nn_idx Answer vector
\param eps Error tolerence, default of 0.0.
*/
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, float eps = 0)
{
NN.ksearch(q,k,nn_idx,eps);
};
/*!
\brief Nearest Neighbor search function
Searches for the k nearest neighbors to the point q. The answer
vector returned will contain the indexes to the answer points.
The distance vector contains the square distances to the answer points.
This function is thread-safe.
\param q The query point
\param k The number of neighbors to return
\param nn_idx Answer vector
\param dist Distance vector
\param eps Error tolerence, default of 0.0.
*/
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, vector<double> &dist, float eps=0)
{
NN.ksearch(q,k,nn_idx,dist,eps);
};
private:
sfcnn_work<reviver::dpoint<NumType, Dim> > NN;
void sfcnn_init(Point *points, long int N, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
for(long int i=0;i < N;++i)
{
NN.pointers[i] = (i);
for(unsigned int j=0;j < Dim;++j)
NN.points[i][j] = points[i][j];
}
NN.sfcnn_work_init(num_threads);
};
};
template <typename Point, unsigned int Dim>
class sfcnn <Point, Dim, float>
{
public:
sfcnn(){};
sfcnn(Point *points, long int N, int num_threads=1)
{
sfcnn_init(points,N,num_threads);
};
~sfcnn(){};
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, float eps = 0)
{
reviver::dpoint<sep_float<float>, Dim> Q;
for(unsigned int i=0;i < Dim;++i)
Q[i] = q[i];
NN.ksearch(Q,k,nn_idx,eps);
};
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, vector<double> &dist, float eps=0)
{
reviver::dpoint<sep_float<float>, Dim> Q;
for(unsigned int i=0;i < Dim;++i)
Q[i] = q[i];
NN.ksearch(Q,k,nn_idx,dist,eps);
};
;
private:
sfcnn_work<reviver::dpoint<sep_float<float>, Dim> > NN;
void sfcnn_init(Point *points, long int N, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
for(long int i=0;i < N;++i)
{
NN.pointers[i] = i;
for(unsigned int j=0;j < Dim;++j)
{
NN.points[i][j] = points[i][j];
}
}
NN.sfcnn_work_init(num_threads);
};
};
template <typename Point, unsigned int Dim>
class sfcnn <Point, Dim, double>
{
public:
sfcnn(){};
sfcnn(Point *points, long int N, int num_threads=1)
{
sfcnn_init(points,N,num_threads);
};
~sfcnn(){};
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, float eps = 0)
{
reviver::dpoint<sep_float<double>, Dim> Q;
for(unsigned int i=0;i < Dim;++i)
Q[i] = q[i];
NN.ksearch(Q,k,nn_idx,eps);
};
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, vector<double> &dist, float eps=0)
{
reviver::dpoint<sep_float<double>, Dim> Q;
for(unsigned int i=0;i < Dim;++i)
Q[i] = q[i];
NN.ksearch(Q,k,nn_idx,dist,eps);
};
private:
void sfcnn_init(Point *points, long int N, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
for(long int i=0;i < N;++i)
{
NN.pointers[i] = i;
for(unsigned int j=0;j < Dim;++j)
NN.points[i][j] = points[i][j];
}
NN.sfcnn_work_init(num_threads);
};
sfcnn_work<reviver::dpoint<sep_float<double>, Dim> > NN;
};
template <typename Point, unsigned int Dim>
class sfcnn <Point, Dim, long double>
{
public:
sfcnn(){};
sfcnn(Point *points, long int N, int num_threads=1)
{
sfcnn_init(points,N,num_threads);
};
~sfcnn(){};
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, float eps = 0)
{
reviver::dpoint<sep_float<long double>, Dim> Q;
for(unsigned int i=0;i < Dim;++i)
Q[i] = q[i];
NN.ksearch(Q,k,nn_idx,eps);
};
void ksearch(Point q, unsigned int k, vector<long unsigned int> &nn_idx, vector<double> &dist, float eps=0)
{
reviver::dpoint<sep_float<long double>, Dim> Q;
for(unsigned int i=0;i < Dim;++i)
Q[i] = q[i];
NN.ksearch(Q,k,nn_idx,dist,eps);
};
private:
void sfcnn_init(Point *points, long int N, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
for(long int i=0;i < N;++i)
{
NN.pointers[i] = i;
for(unsigned int j=0;j < Dim;++j)
NN.points[i][j] = points[i][j];
}
NN.sfcnn_work_init(num_threads);
};
sfcnn_work<reviver::dpoint<sep_float<long double>, Dim> > NN;
};
#endif

396
cmvs/program/base/stann/sfcnn_knng.hpp~
Wyświetl plik

@ -1,396 +0,0 @@
/*****************************************************************************/
/* */
/* Header: sfcnn_knng.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
#ifndef __KNNGRAPH__
#define __KNNGRAPH__
#include <cstdlib>
#include <cmath>
#include <climits>
#include <vector>
#include <queue>
#include <algorithm>
#include <pthread.h>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <pair_iter.hpp>
#include <qknn.hpp>
#include <zorder_lt.hpp>
#include <bsearch.hpp>
/*!
\file sfcnn_knng.hpp
\brief Implements K-Nearest Neighbor Graph construction using SFC nearest neighbor algorithm. Construction is done in parallel using OpenMP.
*/
using namespace std;
template <typename Point>
class sfcnn_knng_work
{
public:
sfcnn_knng_work(){};
~sfcnn_knng_work(){};
void sfcnn_knng_work_init(long int N, unsigned int k, int num_threads);
vector<vector<long unsigned int> > answer;
vector<Point> points;
vector<long unsigned int> pointers;
private:
void compute_bounding_box(Point q, Point &q1, Point &q2, double r);
void recurse(int s, // Starting index
int n, // Number of points
long int q, // Query point
qknn &ans, // Answer que
Point &bound_box_lower_corner,
Point &bound_box_upper_corner,
int initial_scan_lower_range,
int initial_scan_upper_range,
zorder_lt<Point> &lt);
typename Point::__NumType max, min;
};
template<typename Point>
void sfcnn_knng_work<Point>::sfcnn_knng_work_init(long int N, unsigned int k, int num_threads)
{
zorder_lt<Point> lt;
Point bound_box_lower_corner,
bound_box_upper_corner;
double distance;
long int range_b;
long int range_e;
max = numeric_limits<typename Point::__NumType>::max();
min = numeric_limits<typename Point::__NumType>::min();
answer.resize(N);
#ifdef _OPENMP
long int chunk = N/num_threads;
omp_set_num_threads(num_threads);
#endif
pair_iter<typename vector<Point>::iterator,
typename vector<long unsigned int>::iterator> a(points.begin(), pointers.begin()),
b(points.end(), pointers.end());
sort(a,b,lt);
vector<qknn> que;
que.resize(N);
#ifdef _OPENMP
#pragma omp parallel private(distance, range_b, range_e)
#endif
{
#ifdef _OPENMP
#pragma omp for schedule(static, chunk)
#endif
for(long int i=0;i < N;++i)
{
range_b = i-k;
if(range_b < 0) range_b = 0;
range_e = i+k+1;
if(range_e > N) range_e = N;
que[i].set_size(k);
for(long int j=range_b;j < i;++j)
{
distance = points[i].sqr_dist(points[j]);
que[i].update(distance, pointers[j]);
}
for(long int j=i+1;j < range_e;++j)
{
distance = points[i].sqr_dist(points[j]);
que[i].update(distance, pointers[j]);
}
}
}
#ifdef _OPENMP
#pragma omp parallel private(distance, range_b, range_e, bound_box_lower_corner, bound_box_upper_corner)
#endif
{
#ifdef _OPENMP
#pragma omp for schedule(static, chunk)
#endif
for(long int i=0;i < N;++i)
{
range_b = i-k;
if(range_b < 0) range_b = 0;
range_e = i+k+1;
if(range_e > N) range_e = N;
compute_bounding_box(points[i], bound_box_lower_corner, bound_box_upper_corner, sqrt(que[i].topdist()));
if(!lt(bound_box_upper_corner, points[range_e]) || !lt(points[range_b], bound_box_lower_corner))
{
recurse(0, points.size(), i, que[i],
bound_box_lower_corner,
bound_box_upper_corner,
range_b,
range_e,
lt);
}
que[i].answer(answer[pointers[i]]);
}
}
points.clear();
pointers.clear();
}
template<typename Point>
void sfcnn_knng_work<Point>::compute_bounding_box(Point q, Point &q1, Point &q2, double R)
{
typename Point::__NumType radius;
radius = (typename Point::__NumType) ceil(R);
for(unsigned int i=0;i<Point::__DIM;++i)
{
if(q[i] < (min+radius)) q1[i] = min;
else q1[i] = q[i]-radius;
if(q[i] > (max-radius)) q2[i] = max;
else q2[i] = q[i]+radius;
}
}
template<typename Point>
void sfcnn_knng_work<Point>::recurse(int s, // Starting index
int n, // Number of points
long int q,
qknn &ans, // Answer que
Point &bound_box_lower_corner,
Point &bound_box_upper_corner,
int initial_scan_lower_range,
int initial_scan_upper_range,
zorder_lt<Point> &lt)
{
double distance;
if(n < 4)
{
if(n == 0) return;
bool update=false;
for(int i=0;i < n;++i)
{
if((s+i >= initial_scan_lower_range)
&& (s+i < initial_scan_upper_range))
continue;
distance = points[q].sqr_dist(points[s+i]);
update = ans.update(distance, pointers[s+i]) || update;
}
if(update)
compute_bounding_box(points[q], bound_box_lower_corner, bound_box_upper_corner, sqrt(ans.topdist()));
return;
}
if((s+n/2 >= initial_scan_lower_range) && (s+n/2 < initial_scan_upper_range))
{
}
else
{
distance = points[q].sqr_dist(points[s+n/2]);
if(ans.update(distance, pointers[s+n/2]))
compute_bounding_box(points[q], bound_box_lower_corner, bound_box_upper_corner, sqrt(ans.topdist()));
}
if((lt.dist_sq_to_quad_box(points[q],points[s], points[s+n-1])) > ans.topdist())
return;
if(lt(points[q],points[s+n/2]))
{
recurse(s, n/2, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range,
lt);
if(lt(points[s+n/2],bound_box_upper_corner))
recurse(s+n/2+1,n-n/2-1, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range,
lt);
}
else
{
recurse(s+n/2+1, n-n/2-1, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range,
lt);
if(lt(bound_box_lower_corner,points[s+n/2]))
recurse(s, n/2, q, ans,
bound_box_lower_corner,
bound_box_upper_corner,
initial_scan_lower_range,
initial_scan_upper_range,
lt);
}
};
template <typename Point, unsigned int Dim, typename NumType>
class sfcnn_knng
{
public:
sfcnn_knng(){};
sfcnn_knng(Point *points, long int N, unsigned int k, int num_threads=1)
{
sfcnn_knng_init(points,N,k,num_threads);
};
~sfcnn_knng(){};
vector<long unsigned int>& operator[](long unsigned int i)
{
return NN.answer[i];
};
private:
void sfcnn_knng_init(Point *points, long int N, unsigned int k, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
#ifdef _OPENMP
long int chunk = N/num_threads;
omp_set_num_threads(num_threads);
#pragma omp parallel for schedule(static, chunk)
#endif
for(int i=0;i < N;++i)
{
NN.pointers[i] = (i);
for(unsigned int j=0;j < Dim;++j)
NN.points[i][j] = points[i][j];
}
NN.sfcnn_knng_work_init(N, k, num_threads);
};
sfcnn_knng_work<reviver::dpoint<NumType, Dim> > NN;
};
template <typename Point, unsigned int Dim>
class sfcnn_knng <Point, Dim, float>
{
public:
sfcnn_knng(){};
sfcnn_knng(Point *points, long int N, unsigned int k, int num_threads=1)
{
sfcnn_knng_init(points,N,k,num_threads);
};
~sfcnn_knng(){};
vector<long unsigned int>& operator[](long unsigned int i)
{
return NN.answer[i];
};
private:
void sfcnn_knng_init(Point *points, long int N, unsigned int k, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
#ifdef _OPENMP
long int chunk = N/num_threads;
omp_set_num_threads(num_threads);
#pragma omp parallel for schedule(static, chunk)
#endif
for(int i=0;i < N;++i)
{
NN.pointers[i] = i;
for(unsigned int j=0;j < Dim;++j)
{
NN.points[i][j] = points[i][j];
}
}
NN.sfcnn_knng_work_init(N,k,num_threads);
};
sfcnn_knng_work<reviver::dpoint<sep_float<float>, Dim> > NN;
};
template <typename Point, unsigned int Dim>
class sfcnn_knng <Point, Dim, double>
{
public:
sfcnn_knng(){};
sfcnn_knng(Point *points, long int N, unsigned int k, int num_threads=1)
{
sfcnn_knng_init(points,N,k,num_threads);
};
~sfcnn_knng(){};
vector<long unsigned int>& operator[](long unsigned int i)
{
return NN.answer[i];
};
private:
void sfcnn_knng_init(Point *points, long int N, unsigned int k, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
#ifdef _OPENMP
long int chunk = N/num_threads;
omp_set_num_threads(num_threads);
#pragma omp parallel for schedule(static, chunk)
#endif
for(int i=0;i < N;++i)
{
NN.pointers[i] = i;
for(unsigned int j=0;j < Dim;++j)
NN.points[i][j] = points[i][j];
}
NN.sfcnn_knng_work_init(N,k,num_threads);
};
sfcnn_knng_work<reviver::dpoint<sep_float<double>, Dim> > NN;
};
template <typename Point, unsigned int Dim>
class sfcnn_knng <Point, Dim, long double>
{
public:
sfcnn_knng(){};
sfcnn_knng(Point *points, long int N, unsigned int k, int num_threads=1)
{
sfcnn_knng_init(points,N,k,num_threads);
};
~sfcnn_knng(){};
vector<long unsigned int>& operator[](long unsigned int i)
{
return NN.answer[i];
};
private:
void sfcnn_knng_init(Point *points, long int N, unsigned int k, int num_threads=1)
{
NN.points.resize(N);
NN.pointers.resize(N);
#ifdef _OPENMP
long int chunk = N/num_threads;
omp_set_num_threads(num_threads);
#pragma omp parallel for schedule(static, chunk)
#endif
for(int i=0;i < N;++i)
{
NN.pointers[i] = i;
for(unsigned int j=0;j < Dim;++j)
NN.points[i][j] = points[i][j];
}
NN.sfcnn_knng_work_init(N,k,num_threads);
};
sfcnn_knng_work<reviver::dpoint<sep_float<long double>, Dim> > NN;
};
#endif

143
cmvs/program/base/stann/test.hpp~
Wyświetl plik

@ -1,143 +0,0 @@
/*****************************************************************************/
/* */
/* Header: test.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
#ifndef STANN_TEST
#define STANN_TEST
#include <iostream>
#include <vector>
#include <bruteNN.hpp>
#include <dpoint.hpp>
#include <rand.hpp>
#include <sfcnn.hpp>
#include <sfcnn_knng.hpp>
using namespace std;
/*! \file test.hpp
\brief Header file for the test.cpp program.*/
template<typename Point, typename T>
Point newRandomPoint(T Min, T Max)
{
double d;
Point a;
double max, min;
max = (double) Max / (double) numeric_limits<T>::max();
min = (double) Min / (double) numeric_limits<T>::max();
for(unsigned int i=0;i < Point::__DIM;++i)
{
d = __drand48__();
d = d*(max-min)-max;
d *= (double) numeric_limits<T>::max();
d *= -1;
a[i] = (T) d;
}
return a;
}
template<typename T, unsigned DIM>
bool testNN(unsigned int Size, unsigned int k, T min, T max)
{
typedef reviver::dpoint<T, DIM> Point;
vector<Point> data;
vector<Point> query;
vector<long unsigned int> sfcnn_ans;
vector<long unsigned int> bf_ans;
data.resize(Size);
query.resize(Size);
for(unsigned int i=0;i < data.size();++i)
{
data[i] = newRandomPoint<Point, T>(min, max);
query[i] = newRandomPoint<Point, T>(min, max);
}
bruteNN<Point> BF(&data[0], data.size());
sfcnn<Point, DIM, T> SFC(&data[0], data.size());
for(unsigned int i=0;i < data.size();++i)
{
BF.ksearch(query[i], k, bf_ans);
SFC.ksearch(query[i], k, sfcnn_ans);
for(unsigned int j=0;j < Point::__DIM;++j)
{
if(bf_ans[j] != sfcnn_ans[j])
{
/*
cerr << "SFCNN:" << endl;
for(unsigned int q=0;q < sfcnn_ans.size();++q)
{
cerr << sfcnn_ans[q] << endl;
}
cerr << "BF:" << endl;
for(unsigned int q=0;q < bf_ans.size();++q)
{
cerr << bf_ans[q] << endl;
}
*/
return false;
}
}
}
return true;
}
template<typename T, unsigned DIM>
bool testKNNG(unsigned int Size, unsigned int k, T min, T max, int num_threads)
{
typedef reviver::dpoint<T, DIM> Point;
vector<Point> data;
vector<long unsigned int> bf_ans;
data.resize(Size);
for(unsigned int i=0;i < data.size();++i)
{
data[i] = newRandomPoint<Point, T>(min, max);
}
bruteNN<Point> BF(&data[0], data.size());
sfcnn_knng<Point, DIM, T> SFC(&data[0], data.size(), k, num_threads);
for(unsigned int i=0;i < data.size();++i)
{
BF.ksearch(data[i], k+1, bf_ans);
for(unsigned int j=1;j < k+1;++j)
{
if(bf_ans[j] != SFC[i][j-1])
{
cerr << "SFCNN:" << endl;
for(unsigned int q=0;q < SFC[i].size();++q)
{
cerr << SFC[i][q] << endl;
}
cerr << "BF:" << endl;
for(unsigned int q=0;q < bf_ans.size();++q)
{
cerr << bf_ans[q] << endl;
}
return false;
}
}
}
return true;
}
#endif

275
cmvs/program/base/stann/zorder_lt.hpp~
Wyświetl plik

@ -1,275 +0,0 @@
/*****************************************************************************/
/* */
/* Header: zorder_lt.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
#ifndef __SFCNN_ZORDER_LT__
#define __SFCNN_ZORDER_LT__
#include <dpoint.hpp>
#include <cmath>
#include <climits>
#include <iostream>
#include <zorder_type_traits.hpp>
#include <sep_float.hpp>
#include <pair_iter.hpp>
using namespace std;
template<typename Point>
class zorder_lt;
template<typename Point, typename CType, typename sign_trait, typename integral_trait, typename sep_trait>
class zorder_lt_worker;
//Z Order spec for unsigned integral types
template<typename Point, typename CType>
class zorder_lt_worker <Point, CType, zorder_f, zorder_t, zorder_f>
{
public:
zorder_lt_worker(){;}
~zorder_lt_worker(){;}
bool operator()(const Point &p, const Point &q)
{
return lt_func(p,q);
}
double dist_sq_to_quad_box(Point &q, Point &p1, Point &p2)
{
unsigned int j;
int i;
CType x,y;
double z;
z=0;
for(j=x=0;j < Point::__DIM;++j)
{
y = p1[j]^p2[j];
if(less_msb(x, y))
{
x = y;
}
}
frexp((float)x, &i);
for(j=0;j < Point::__DIM;++j)
{
x = (((p1)[j])>>i)<<i;
y = x+(1 << i);
if(q[j] < x)
z+= pow(((double) q[j]-(double) x), 2.0);
else if(q[j] > y)
z+= pow(((double) q[j]-(double) y), 2.0);
}
//cout << q << endl << p1 << endl << p2 << endl << z << endl;
return z;
}
private:
bool lt_func(const Point &p, const Point &q)
{
CType j,x,y,k;
for(j=k=x=0;k < Point::__DIM;++k)
{
y = (p[k])^(q[k]);
if(less_msb(x,y))
{
j=k;
x=y;
}
}
return p[j] < q[j];
}
bool less_msb(CType x, CType y)
{
return (x < y) && (x < (x^y));
}
};
//Z Order spec for signed integral types
template<typename Point, typename CType>
class zorder_lt_worker <Point, CType, zorder_t, zorder_t, zorder_f>
{
public:
zorder_lt_worker(){;}
~zorder_lt_worker(){;}
bool operator()(const Point &p, const Point &q)
{
return lt_func(p,q);
}
double dist_sq_to_quad_box(Point &q, Point &p1, Point &p2)
{
unsigned int j;
int i;
CType x,y;
double z, X, Y;
x = numeric_limits<CType>::min();
for(j=0;j < Point::__DIM;++j)
{
if((p1[j] < 0) != (p2[j] < 0))
{
//cout << "Break out" << endl;
//cout << "P1: " << p1 << endl << "P2: " << p2 << endl;
return 0;
}
y = p1[j]^p2[j];
if(less_msb(x,y)) x = y;
}
frexp((double)x, &i);
//cout << "i: " << i << endl;
//int exp;
for(z=j=0;j < Point::__DIM;++j)
{
X = (p1[j]>>i)<<i;
//frexp(p1[j], &exp);
//cout << "p1[j]" << (p1[j] >> i) << endl;
//cout << "exp: " << exp << endl;
//cout << "X: " << X << endl;
Y = X+(1 << i);
if(q[j] < X)
z+=pow(((double) X - (double) q[j]), 2.0);
else if(q[j] > Y)
z+=pow(((double) q[j] - (double) Y), 2.0);
}
return z;
}
private:
bool lt_func(const Point &p, const Point &q)
{
CType j,x,y;
unsigned int k;
x = numeric_limits<CType>::min();
for(j=k=0;k < Point::__DIM;++k)
{
if((p[k] < 0) != (q[k] < 0))
return p[k] < q[k];
y = (p[k])^(q[k]);
if(less_msb(x,y))
{
j=k;
x=y;
}
}
return p[j] < q[j];
}
bool less_msb(CType x, CType y)
{
return (x < y) && (x < (x^y));
}
};
//Z Order spec for seperated floating point types
template<typename Point, typename CType>
class zorder_lt_worker <Point, CType, zorder_t, zorder_f, zorder_t>
{
public:
zorder_lt_worker()
{
}
~zorder_lt_worker(){;}
bool operator()(const Point &p, const Point &q)
{
return lt_func(p,q);
}
double dist_sq_to_quad_box(Point &q, Point &p1, Point &p2)
{
unsigned int j;
int x,y;
double box_edge_1, box_edge_2;
double z;
typename CType::flt_type box_dist;
z = 0;
x = 0;
//cout << "X: " << x << endl;
for(j=0;j < Point::__DIM;++j)
{
if((p1[j].get_flt() < 0) != (p2[j].get_flt() < 0))
return 0;
y = msdb(p1[j], p2[j]);
if(y > x)
{
x = y;
}
}
box_dist = pow(2.0,x);
for(j=0;j < Point::__DIM;++j)
{
//cout << "p1[j]/2**i: " << floor(p1[j] / box_dist) << endl;
//box_edge_1 = p1[j].get_flt() - fmod(p1[j].get_flt(),box_dist);
box_edge_1 = floor(p1[j] / box_dist) * box_dist;
//cout << "Box1: " << box_edge_1 << endl;
box_edge_2 = box_edge_1+box_dist;
if(q[j].get_flt() < box_edge_1)
z+=(q[j].get_flt()-box_edge_1)*(q[j].get_flt()-box_edge_1);
else if(q[j].get_flt() > box_edge_2)
z+= (q[j].get_flt()-box_edge_2)*(q[j].get_flt()-box_edge_2);
}
return z;
}
private:
bool lt_func(const Point &p, const Point &q)
{
int y,x;
unsigned int k,j;
j = 0;
x = 0;
for(k=0;k < Point::__DIM;++k)
{
if((p[k].get_flt() < 0) != (q[k].get_flt() < 0))
return p[k].get_flt() < q[k].get_flt();
y = msdb(p[k], q[k]);
if(x < y)
{
j = k;
x = y;
}
}
return p[j] < q[j];
}
};
template<typename Point>
class zorder_lt
{
public:
zorder_lt()
{
zorder_traits<typename Point::__NumType>::check_type();
}
~zorder_lt(){;}
bool operator()(const Point &p, const Point &q)
{
return lt(p, q);
}
bool operator()(const Mypair<typename vector<Point>::iterator,typename vector<long unsigned int>::iterator> &p, const Mypair<typename vector<Point>::iterator, typename vector<long unsigned int>::iterator> &q)
{
return lt(p.val1, q.val1);
}
double dist_sq_to_quad_box(Point &q, Point &p1, Point &p2)
{
return lt.dist_sq_to_quad_box(q,p1,p2);
}
private:
zorder_lt_worker<Point, typename Point::__NumType, typename zorder_traits<typename Point::__NumType>::is_signed, typename zorder_traits<typename Point::__NumType>::is_integral, typename zorder_traits<typename Point::__NumType>::is_seperated> lt;
};
#endif

208
cmvs/program/base/stann/zorder_type_traits.hpp~
Wyświetl plik

@ -1,208 +0,0 @@
/*****************************************************************************/
/* */
/* Header: zorder_type_traits.hpp */
/* */
/* Accompanies STANN Version 0.5 Beta */
/* Aug 05, 2008 */
/* */
/* Copyright 2007, 2008 */
/* Michael Connor and Piyush Kumar */
/* Florida State University */
/* Tallahassee FL, 32306-4532 */
/* */
/*****************************************************************************/
#ifndef __ZORDER_TYPE_TRAITS__
#define __ZORDER_TYPE_TRAITS__
#include <iostream>
#include <sep_float.hpp>
using namespace std;
struct zorder_t {};
struct zorder_f {};
template<typename TYPE>
class zorder_traits
{
public:
static void check_type()
{
//cerr << "Error: Type traits not defined." << endl;
//cerr << "Please ensure the appropriate type is added to zorder_traits.hpp" << endl;
//cerr << "and contact author for update." << endl;
//exit(1);
}
typedef zorder_t is_signed;
typedef zorder_t is_integral;
typedef int unsigned_type;
typedef unsigned int signed_type;
};
template<>
class zorder_traits <int>
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef unsigned int unsigned_type;
typedef int signed_type;
};
template<>
class zorder_traits <unsigned int>
{
public:
static void check_type(){};
typedef zorder_f is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef unsigned int unsigned_type;
typedef int signed_type;
};
template<>
class zorder_traits <char>
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef unsigned char unsigned_type;
typedef char signed_type;
};
template<>
class zorder_traits <unsigned char>
{
public:
static void check_type(){};
typedef zorder_f is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef unsigned char unsigned_type;
typedef char signed_type;
};
template<>
class zorder_traits <short>
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef unsigned short unsigned_type;
typedef short signed_type;
};
template<>
class zorder_traits <unsigned short>
{
public:
static void check_type(){};
typedef zorder_f is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef unsigned short unsigned_type;
typedef unsigned short signed_type;
};
template<>
class zorder_traits <long>
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef long unsigned int unsigned_type;
typedef long int signed_type;
};
template<>
class zorder_traits <unsigned long>
{
public:
static void check_type(){};
typedef zorder_f is_signed;
typedef zorder_t is_integral;
typedef zorder_f is_seperated;
typedef long unsigned int unsigned_type;
typedef long unsigned int signed_type;
};
template<>
class zorder_traits <float>
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_f is_integral;
typedef zorder_f is_seperated;
typedef float unsigned_type;
typedef float signed_type;
};
template<>
class zorder_traits <double>
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_f is_integral;
typedef zorder_f is_seperated;
typedef double unsigned_type;
typedef double signed_type;
};
template<>
class zorder_traits <long double>
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_f is_integral;
typedef zorder_f is_seperated;
typedef long double unsigned_type;
typedef long double signed_type;
};
template<>
class zorder_traits <sep_float<float> >
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_f is_integral;
typedef zorder_t is_seperated;
};
template<>
class zorder_traits <sep_float<double> >
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_f is_integral;
typedef zorder_t is_seperated;
};
template<>
class zorder_traits <sep_float<long double> >
{
public:
static void check_type(){};
typedef zorder_t is_signed;
typedef zorder_f is_integral;
typedef zorder_t is_seperated;
};
#endif

530
odm_georef/src/Georef.cpp~
Wyświetl plik

@ -1,530 +0,0 @@
// PCL
#include <pcl/io/obj_io.h>
#include <pcl/common/transforms.h>
// Modified PCL
#include "modifiedPclFunctions.hpp"
// This
#include "Georef.hpp"
std::ostream& operator<<(std::ostream &os, const GeorefSystem &geo)
{
return os << geo.system_ << " " << static_cast<int>(geo.falseEasting_) << " " << static_cast<int>(geo.falseNorthing_);
}
GeorefCamera::GeorefCamera()
:focalLength_(0.0), k1_(0.0), k2_(0.0), transform_(NULL), position_(NULL)
{
}
GeorefCamera::GeorefCamera(const GeorefCamera &other)
: focalLength_(other.focalLength_), k1_(other.k1_), k2_(other.k2_),
easting_(other.easting_), northing_(other.northing_), altitude_(other.altitude_),
transform_(NULL), position_(NULL)
{
if(NULL != other.transform_)
{
transform_ = new Eigen::Affine3f(*other.transform_);
}
if(NULL != other.position_)
{
position_ = new Eigen::Vector3f(*other.position_);
}
}
GeorefCamera::~GeorefCamera()
{
if(NULL != transform_)
{
delete transform_;
transform_ = NULL;
}
if(NULL != position_)
{
delete position_;
position_ = NULL;
}
}
void GeorefCamera::extractCamera(std::ifstream &bundleStream)
{
// Extract intrinsic parameters.
bundleStream >> focalLength_ >> k1_ >> k2_;
Eigen::Vector3f t;
Eigen::Matrix3f rot;
Eigen::Affine3f transform;
bundleStream >> transform(0,0); // Read rotation (0,0) from bundle file
bundleStream >> transform(0,1); // Read rotation (0,1) from bundle file
bundleStream >> transform(0,2); // Read rotation (0,2) from bundle file
bundleStream >> transform(1,0); // Read rotation (1,0) from bundle file
bundleStream >> transform(1,1); // Read rotation (1,1) from bundle file
bundleStream >> transform(1,2); // Read rotation (1,2) from bundle file
bundleStream >> transform(2,0); // Read rotation (2,0) from bundle file
bundleStream >> transform(2,1); // Read rotation (2,1) from bundle file
bundleStream >> transform(2,2); // Read rotation (2,2) from bundle file
bundleStream >> t(0); // Read translation (0,3) from bundle file
bundleStream >> t(1); // Read translation (1,3) from bundle file
bundleStream >> t(2); // Read translation (2,3) from bundle file
rot = transform.matrix().topLeftCorner<3,3>();
// Calculate translation according to -R't and store in vector.
t = -rot.transpose()*t;
transform(0,3) = t(0);
transform(1,3) = t(1);
transform(2,3) = t(2);
// Set transform and position.
if(NULL != transform_)
{
delete transform_;
transform_ = NULL;
}
transform_ = new Eigen::Affine3f(transform);
if(NULL != position_)
{
delete position_;
position_ = NULL;
}
position_ = new Eigen::Vector3f(t);
}
void GeorefCamera::extractCameraGeoref(std::istringstream &coordStream)
{
coordStream >> easting_ >> northing_ >> altitude_;
}
Vec3 GeorefCamera::getPos()
{
return Vec3((*position_)(0),(*position_)(1),(*position_)(2));
}
Vec3 GeorefCamera::getReferencedPos()
{
return Vec3(easting_,northing_,altitude_);
}
std::ostream& operator<<(std::ostream &os, const GeorefCamera &cam)
{
os << "Focal, k1, k2 : " << cam.focalLength_ << ", " << cam.k1_ << ", " << cam.k2_ << "\n";
if(NULL != cam.transform_)
{
os << "Transform :\n" << cam.transform_->matrix() << "\n";
}
else
{
os << "Transform :\nNULL\n";
}
if(NULL != cam.position_)
{
os << "Position :\n" << cam.position_->matrix() << "\n";
}
else
{
os << "Position :\nNULL\n";
}
os << "east, north, alt : " << cam.easting_ << ", " << cam.northing_ << ", " << cam.altitude_ << '\n';
return os;
}
Georef::Georef()
{
log_.setIsPrintingInCout(true);
bundleFilename_ = "";
coordFilename_ = "";
inputObjFilename_ = "";
outputObjFilename_ = "";
}
Georef::~Georef()
{
}
int Georef::run(int argc, char *argv[])
{
try
{
parseArguments(argc, argv);
makeGeoreferencedModel();
}
catch (const GeorefException& e)
{
log_ << e.what() << "\n";
log_.print(logFile_);
return EXIT_FAILURE;
}
catch (const std::exception& e)
{
log_ << "Error in Georef:\n";
log_ << e.what() << "\n";
log_.print(logFile_);
return EXIT_FAILURE;
}
catch (...)
{
log_ << "Unknown error, terminating:\n";
log_.print(logFile_);
return EXIT_FAILURE;
}
log_.print(logFile_);
return EXIT_SUCCESS;
}
void Georef::parseArguments(int argc, char *argv[])
{
bool outputSpecified = false;
logFile_ = std::string(argv[0]) + "_log.txt";
log_ << logFile_ << "\n";
// If no arguments were passed, print help.
if (argc == 1)
{
printHelp();
}
log_ << "Arguments given\n";
for(int argIndex = 1; argIndex < argc; ++argIndex)
{
log_ << argv[argIndex] << '\n';
}
log_ << '\n';
for(int argIndex = 1; argIndex < argc; ++argIndex)
{
// The argument to be parsed.
std::string argument = std::string(argv[argIndex]);
if(argument == "-help")
{
printHelp();
}
else if(argument == "-verbose")
{
log_.setIsPrintingInCout(true);
}
else if(argument == "-bundleFile" && argIndex < argc)
{
argIndex++;
if (argIndex >= argc)
{
throw GeorefException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
bundleFilename_ = std::string(argv[argIndex]);
log_ << "Reading cameras from: " << bundleFilename_ << "\n";
}
else if(argument == "-coordFile" && argIndex < argc)
{
argIndex++;
if (argIndex >= argc)
{
throw GeorefException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
coordFilename_ = std::string(argv[argIndex]);
log_ << "Reading cameras georeferenced positions from: " << coordFilename_ << "\n";
}
else if(argument == "-inputFile" && argIndex < argc)
{
argIndex++;
if (argIndex >= argc)
{
throw GeorefException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
inputObjFilename_ = std::string(argv[argIndex]);
log_ << "Reading textured mesh from: " << inputObjFilename_ << "\n";
}
else if(argument == "-outputFile" && argIndex < argc)
{
argIndex++;
if (argIndex >= argc)
{
throw GeorefException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
outputObjFilename_ = std::string(argv[argIndex]);
log_ << "Writing output to: " << outputObjFilename_ << "\n";
outputSpecified = true;
}
else
{
printHelp();
throw GeorefException("Unrecognised argument '" + argument + "'");
}
}
if(!outputSpecified)
{
makeDefaultOutput();
}
}
void Georef::printHelp()
{
bool printInCoutPop = log_.isPrintingInCout();
log_.setIsPrintingInCout(true);
log_ << "Georef.exe\n\n";
log_ << "Purpose:" << "\n";
log_ << "Create an orthograpical photo from an oriented textured mesh." << "\n";
log_ << "Usage:" << "\n";
log_ << "The program requires a path to a camera bundle file, a camera georeference coords file, and an input OBJ mesh file. All other input parameters are optional." << "\n\n";
log_ << "The following flags are available\n";
log_ << "Call the program with flag \"-help\", or without parameters to print this message, or check any generated log file.\n";
log_ << "Call the program with flag \"-verbose\", to print log messages in the standard output stream as well as in the log file.\n\n";
log_ << "Parameters are specified as: \"-<argument name> <argument>\", (without <>), and the following parameters are configureable: " << "\n";
log_ << "\"-bundleFile <path>\" (mandatory)" << "\n";
log_ << "\"Input cameras bundle file.\n\n";
log_ << "\"-coordFile <path>\" (mandatory)" << "\n";
log_ << "\"Input cameras geroreferenced coords file.\n\n";
log_ << "\"-inputFile <path>\" (mandatory)" << "\n";
log_ << "\"Input obj file that must contain a textured mesh.\n\n";
log_ << "\"-outputFile <path>\" (optional, default <inputFile>_geo)" << "\n";
log_ << "\"Output obj file that will contain the georeferenced texture mesh.\n\n";
log_.setIsPrintingInCout(printInCoutPop);
}
void Georef::makeDefaultOutput()
{
if(inputObjFilename_.empty())
{
throw GeorefException("Tried to generate default ouptut file without having an input file.");
}
std::string tmp = inputObjFilename_;
size_t findPos = tmp.find_last_of(".");
if(std::string::npos == findPos)
{
throw GeorefException("Tried to generate default ouptut file, could not find .obj in the input file:\n\'"+inputObjFilename_+"\'");
}
tmp = tmp.substr(0, findPos);
outputObjFilename_ = tmp + "_geo.obj";
log_ << "Writing output to: " << outputObjFilename_ << "\n";
}
void Georef::makeGeoreferencedModel()
{
// Read translations from bundle file
std::ifstream bundleStream(bundleFilename_.c_str());
if (!bundleStream.good())
{
throw GeorefException("Failed opening " + bundleFilename_ + " for reading." + '\n');
}
// Read Cameras.
std::string bundleLine;
std::getline(bundleStream, bundleLine); // Read past bundle version comment
int numCameras, numPoints;
bundleStream >> numCameras >> numPoints;
for (int i=0; i<numCameras; ++i)
{
cameras_.push_back(GeorefCamera());
cameras_.back().extractCamera(bundleStream);
}
// Read coords from coord file generated by extract_utm tool
std::ifstream coordStream(coordFilename_.c_str());
if (!coordStream.good())
{
throw GeorefException("Failed opening " + coordFilename_ + " for reading." + '\n');
}
std::string coordString;
std::getline(coordStream, georefSystem_.system_); // System
{
std::getline(coordStream, coordString); // Flase easting & northing.
std::stringstream ss(coordString);
ss >> georefSystem_.falseEasting_ >> georefSystem_.falseNorthing_;
}
log_ << '\n';
log_ << "Geographical reference system\n";
log_ << georefSystem_ << '\n';
// The number of cameras in the coords file.
size_t nGeorefCameras = 0;
// Read the georefernced position for all cameras.
while (std::getline(coordStream, coordString))
{
if(nGeorefCameras >= cameras_.size())
{
throw GeorefException("Error to many cameras in \'" + coordFilename_ + "\' coord file.\n");
}
std::istringstream istr(coordString);
cameras_[nGeorefCameras].extractCameraGeoref(istr);
++nGeorefCameras;
}
coordStream.close();
if(nGeorefCameras < cameras_.size())
{
throw GeorefException("Not enough cameras in \'" + coordFilename_ + "\' coord file.\n");
}
// The optimal camera triplet.
size_t cam0, cam1, cam2;
log_ << '\n';
log_ << "Choosing optimal camera triplet...\n";
chooseBestCameraTriplet(cam0, cam1, cam2);
log_ << "... optimal camera triplet chosen:\n";
log_ << cam0 << ", " << cam1 << ", " << cam2 << '\n';
log_ << '\n';
FindTransform transFinal;
transFinal.findTransform(cameras_[cam0].getPos(), cameras_[cam1].getPos(), cameras_[cam2].getPos(),
cameras_[cam0].getReferencedPos(), cameras_[cam1].getReferencedPos(), cameras_[cam2].getReferencedPos());
log_ << "Final transform:\n";
log_ << transFinal.transform_ << '\n';
// The tranform used to move the chosen area into the ortho photo.
Eigen::Transform<float, 3, Eigen::Affine> transform;
transform(0, 0) = static_cast<float>(transFinal.transform_.r1c1_);
transform(1, 0) = static_cast<float>(transFinal.transform_.r2c1_);
transform(2, 0) = static_cast<float>(transFinal.transform_.r3c1_);
transform(3, 0) = static_cast<float>(transFinal.transform_.r4c1_);
transform(0, 1) = static_cast<float>(transFinal.transform_.r1c2_);
transform(1, 1) = static_cast<float>(transFinal.transform_.r2c2_);
transform(2, 1) = static_cast<float>(transFinal.transform_.r3c2_);
transform(3, 1) = static_cast<float>(transFinal.transform_.r4c2_);
transform(0, 2) = static_cast<float>(transFinal.transform_.r1c3_);
transform(1, 2) = static_cast<float>(transFinal.transform_.r2c3_);
transform(2, 2) = static_cast<float>(transFinal.transform_.r3c3_);
transform(3, 2) = static_cast<float>(transFinal.transform_.r4c3_);
transform(0, 3) = static_cast<float>(transFinal.transform_.r1c4_);
transform(1, 3) = static_cast<float>(transFinal.transform_.r2c4_);
transform(2, 3) = static_cast<float>(transFinal.transform_.r3c4_);
transform(3, 3) = static_cast<float>(transFinal.transform_.r4c4_);
log_ << '\n';
log_ << "Reading mesh file...\n";
// The textureds mesh.e
pcl::TextureMesh mesh;
pcl::io::loadOBJFile(inputObjFilename_, mesh);
log_ << ".. mesh file read.\n";
// Contains the vertices of the mesh.
pcl::PointCloud<pcl::PointXYZ>::Ptr meshCloud (new pcl::PointCloud<pcl::PointXYZ>);
pcl::fromPCLPointCloud2 (mesh.cloud, *meshCloud);
log_ << '\n';
log_ << "Applying transform to mesh...\n";
// Move the mesh into position.
pcl::transformPointCloud(*meshCloud, *meshCloud, transform);
log_ << ".. mesh transformed.\n";
// Update the mesh.
pcl::toPCLPointCloud2 (*meshCloud, mesh.cloud);
// Iterate over each part of the mesh (one per material), to make texture file paths relative the .mtl file.
for(size_t t = 0; t < mesh.tex_materials.size(); ++t)
{
// The material of the current submesh.
pcl::TexMaterial& material = mesh.tex_materials[t];
size_t find = material.tex_file.find_last_of("/\\");
if(std::string::npos != find)
{
material.tex_file = material.tex_file.substr(find + 1);
}
}
log_ << '\n';
log_ << "Saving mesh file to \'" << outputObjFilename_ << "\'...\n";
saveOBJFile(outputObjFilename_, mesh, 8);
log_ << ".. mesh file saved.\n";
printGeorefSystem();
}
void Georef::chooseBestCameraTriplet(size_t &cam0, size_t &cam1, size_t &cam2)
{
double minTotError = std::numeric_limits<double>::infinity();
for(size_t t = 0; t < cameras_.size(); ++t)
{
for(size_t s = t; s < cameras_.size(); ++s)
{
for(size_t p = s; p < cameras_.size(); ++p)
{
FindTransform trans;
trans.findTransform(cameras_[t].getPos(), cameras_[s].getPos(), cameras_[p].getPos(),
cameras_[t].getReferencedPos(), cameras_[s].getReferencedPos(), cameras_[p].getReferencedPos());
// The total error for the curren camera triplet.
double totError = 0.0;
for(size_t r = 0; r < cameras_.size(); ++r)
{
totError += trans.error(cameras_[r].getPos(), cameras_[r].getReferencedPos());
}
if(minTotError > totError)
{
minTotError = totError;
cam0 = t;
cam1 = s;
cam2 = p;
}
}
}
}
log_ << "Mean georeference error " << minTotError / static_cast<double>(cameras_.size()) << '\n';
}
void Georef::printGeorefSystem()
{
if(outputObjFilename_.empty())
{
throw GeorefException("Output file path empty!.");
}
std::string tmp = outputObjFilename_;
size_t findPos = tmp.find_last_of(".");
if(std::string::npos == findPos)
{
throw GeorefException("Tried to generate default ouptut file, could not find .obj in the output file:\n\'"+outputObjFilename_+"\'");
}
tmp = tmp.substr(0, findPos);
tmp = tmp + "_georef_system.txt";
log_ << '\n';
log_ << "Saving georeference system file to \'" << tmp << "\'...\n";
std::ofstream geoStream(tmp.c_str());
geoStream << georefSystem_ << std::endl;
geoStream.close();
log_ << "... georeference system saved.\n";
}

361
odm_meshing/src/OdmMeshing.cpp~
Wyświetl plik

@ -1,361 +0,0 @@
#include "OdmMeshing.hpp"
OdmMeshing::OdmMeshing() : log_(false)
{
meshCreator_ = pcl::Poisson<pcl::PointNormal>::Ptr(new pcl::Poisson<pcl::PointNormal>());
points_ = pcl::PointCloud<pcl::PointNormal>::Ptr(new pcl::PointCloud<pcl::PointNormal>());
mesh_ = pcl::PolygonMeshPtr(new pcl::PolygonMesh);
decimatedMesh_ = pcl::PolygonMeshPtr(new pcl::PolygonMesh);
// Set default values
outputFile_ = "";
logFilePath_ = "";
maxVertexCount_ = 0;
treeDepth_ = 0;
solverDivide_ = 9.0;
samplesPerNode_ = 1.0;
decimationFactor_ = 0.0;
logFilePath_ = "odm_meshing_log.txt";
log_ << logFilePath_ << "\n";
}
OdmMeshing::~OdmMeshing()
{
}
int OdmMeshing::run(int argc, char **argv)
{
// If no arguments were passed, print help and return early.
if (argc <= 1)
{
printHelp();
return EXIT_SUCCESS;
}
try
{
parseArguments(argc, argv);
loadPoints();
createMesh();
decimateMesh();
writePlyFile();
}
catch (const OdmMeshingException& e)
{
log_.setIsPrintingInCout(true);
log_ << e.what() << "\n";
log_.printToFile(logFilePath_);
log_ << "For more detailed information, see log file." << "\n";
return EXIT_FAILURE;
}
catch (const std::exception& e)
{
log_.setIsPrintingInCout(true);
log_ << "Error in OdmMeshing:\n";
log_ << e.what() << "\n";
log_.printToFile(logFilePath_);
log_ << "For more detailed information, see log file." << "\n";
return EXIT_FAILURE;
}
catch (...)
{
log_.setIsPrintingInCout(true);
log_ << "Unknwon error in OdmMeshing:\n";
log_.printToFile(logFilePath_);
log_ << "For more detailed information, see log file." << "\n";
return EXIT_FAILURE;
}
log_.printToFile(logFilePath_);
return EXIT_SUCCESS;
}
void OdmMeshing::parseArguments(int argc, char **argv)
{
for(int argIndex = 1; argIndex < argc; ++argIndex)
{
// The argument to be parsed.
std::string argument = std::string(argv[argIndex]);
if(argument == "-help")
{
printHelp();
}
else if(argument == "-verbose")
{
log_.setIsPrintingInCout(true);
}
else if(argument == "-maxVertexCount" && argIndex < argc)
{
++argIndex;
if (argIndex >= argc)
{
throw OdmMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
std::stringstream ss(argv[argIndex]);
ss >> maxVertexCount_;
if (ss.bad())
{
throw OdmMeshingException("Argument '" + argument + "' has a bad value (wrong type).");
}
log_ << "Vertex count was manually set to: " << maxVertexCount_ << "\n";
}
else if(argument == "-octreeDepth" && argIndex < argc)
{
++argIndex;
if (argIndex >= argc)
{
throw OdmMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
std::stringstream ss(argv[argIndex]);
ss >> treeDepth_;
if (ss.bad())
{
throw OdmMeshingException("Argument '" + argument + "' has a bad value (wrong type).");
}
log_ << "Octree depth was manually set to: " << treeDepth_ << "\n";
}
else if(argument == "-solverDivide" && argIndex < argc)
{
++argIndex;
if (argIndex >= argc)
{
throw OdmMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
std::stringstream ss(argv[argIndex]);
ss >> solverDivide_;
if (ss.bad())
{
throw OdmMeshingException("Argument '" + argument + "' has a bad value (wrong type).");
}
log_ << "Numerical solver divisions was manually set to: " << treeDepth_ << "\n";
}
else if(argument == "-samplesPerNode" && argIndex < argc)
{
++argIndex;
if (argIndex >= argc)
{
throw OdmMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
std::stringstream ss(argv[argIndex]);
ss >> samplesPerNode_;
if (ss.bad())
{
throw OdmMeshingException("Argument '" + argument + "' has a bad value (wrong type).");
}
log_ << "The number of samples per octree node was manually set to: " << samplesPerNode_ << "\n";
}
else if(argument == "-inputFile" && argIndex < argc)
{
++argIndex;
if (argIndex >= argc)
{
throw OdmMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
inputFile_ = std::string(argv[argIndex]);
std::ifstream testFile(inputFile_.c_str(), std::ios::binary);
if (!testFile.is_open())
{
throw OdmMeshingException("Argument '" + argument + "' has a bad value. (file not accessible)");
}
testFile.close();
log_ << "Reading point cloud at: " << inputFile_ << "\n";
}
else if(argument == "-outputFile" && argIndex < argc)
{
++argIndex;
if (argIndex >= argc)
{
throw OdmMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
outputFile_ = std::string(argv[argIndex]);
std::ofstream testFile(outputFile_.c_str());
if (!testFile.is_open())
{
throw OdmMeshingException("Argument '" + argument + "' has a bad value.");
}
testFile.close();
log_ << "Writing output to: " << outputFile_ << "\n";
}
else if(argument == "-logFile" && argIndex < argc)
{
++argIndex;
if (argIndex >= argc)
{
throw OdmMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
logFilePath_ = std::string(argv[argIndex]);
std::ofstream testFile(outputFile_.c_str());
if (!testFile.is_open())
{
throw OdmMeshingException("Argument '" + argument + "' has a bad value.");
}
testFile.close();
log_ << "Writing log information to: " << logFilePath_ << "\n";
}
else
{
printHelp();
throw OdmMeshingException("Unrecognised argument '" + argument + "'");
}
}
}
void OdmMeshing::loadPoints()
{
if(pcl::io::loadPLYFile<pcl::PointNormal> (inputFile_.c_str(), *points_.get()) == -1) {
throw OdmMeshingException("Error when reading points and normals from:\n" + inputFile_ + "\n");
}
else
{
log_ << "Successfully loaded " << points_->size() << " points with corresponding normals from file.\n";
}
}
void OdmMeshing::printHelp()
{
bool printInCoutPop = log_.isPrintingInCout();
log_.setIsPrintingInCout(true);
log_ << "OpenDroneMapMeshing.exe\n\n";
log_ << "Purpose:" << "\n";
log_ << "Create a mesh from an oriented point cloud (points with normals) using the Poisson surface reconstruction method." << "\n";
log_ << "Usage:" << "\n";
log_ << "The program requires a path to an input PLY point cloud file, all other input parameters are optional." << "\n\n";
log_ << "The following flags are available\n";
log_ << "Call the program with flag \"-help\", or without parameters to print this message, or check any generated log file.\n";
log_ << "Call the program with flag \"-verbose\", to print log messages in the standard output stream as well as in the log file.\n\n";
log_ << "Parameters are specified as: \"-<argument name> <argument>\", (without <>), and the following parameters are configureable: " << "\n";
log_ << "\"-inputFile <path>\" (mandatory)" << "\n";
log_ << "\"Input ascii ply file that must contain a point cloud with normals.\n\n";
log_ << "\"-outputFile <path>\" (optional, default: odm_mesh.ply)" << "\n";
log_ << "\"Target file in which the mesh is saved.\n\n";
log_ << "\"-logFile <path>\" (optional, default: odm_meshing_log.txt)" << "\n";
log_ << "\"Target file in which the mesh is saved.\n\n";
log_ << "\"-maxVertexCount <integer>\" (optional, default: 100,000)" << "\n";
log_ << "Desired final vertex count (after decimation), set to 0 to disable decimation.\n\n";
log_ << "\"-treeDepth <integer>\" (optional, default: 0 (automatic))" << "\n";
log_ << "Controls octree depth used for poisson reconstruction. Recommended values (9-11).\n"
<< "Increasing the value on this parameter will raise initial vertex count."
<< "If omitted or zero, the depth is calculated automatically from the input point count.\n\n";
log_ << "\"-samplesPerNode <float>\" (optional, default: 1)" << "\n";
log_ << "Average number of samples (points) per octree node. Increasing this value might help if data is very noisy.\n\n";
log_ << "\"-solverDivide <integer>\" (optional, default: 9)" << "\n";
log_ << "Ocree depth at which the Laplacian equation is solved in the surface reconstruction step.\n";
log_ << "Increasing this value increases computation times slightly but helps reduce memory usage.\n\n";
log_.setIsPrintingInCout(printInCoutPop);
}
void OdmMeshing::createMesh()
{
// Attempt to calculate the depth of the tree if unspecified
if (treeDepth_ == 0)
{
treeDepth_ = calcTreeDepth(points_->size());
}
log_ << "Octree depth used for reconstruction is: " << treeDepth_ << "\n";
log_ << "Estimated initial vertex count: " << pow(4, treeDepth_) << "\n\n";
meshCreator_->setDepth(treeDepth_);
meshCreator_->setSamplesPerNode(samplesPerNode_);
meshCreator_->setInputCloud(points_);
// Guarantee manifold mesh.
meshCreator_->setManifold(true);
// Begin reconstruction
meshCreator_->reconstruct(*mesh_.get());
log_ << "Reconstruction complete:\n";
log_ << "Vertex count: " << mesh_->cloud.width << "\n";
log_ << "Triangle count: " << mesh_->polygons.size() << "\n\n";
}
void OdmMeshing::decimateMesh()
{
if (maxVertexCount_ <= 0)
{
log_ << "Vertex count not specified, decimation cancelled.\n";
return;
}
if (maxVertexCount_ > mesh_->cloud.height*mesh_->cloud.width)
{
log_ << "Vertex count in mesh lower than initially generated mesh, unable to decimate.\n";
return;
}
else
{
decimatedMesh_ = pcl::PolygonMeshPtr(new pcl::PolygonMesh);
double reductionFactor = 1.0 - double(maxVertexCount_)/double(mesh_->cloud.height*mesh_->cloud.width);
log_ << "Decimating mesh, removing " << reductionFactor*100 << " percent of vertices.\n";
pcl::MeshQuadricDecimationVTK decimator;
decimator.setInputMesh(mesh_);
decimator.setTargetReductionFactor(reductionFactor);
decimator.process(*decimatedMesh_.get());
log_ << "Decimation complete.\n";
log_ << "Decimated vertex count: " << decimatedMesh_->cloud.width << "\n";
log_ << "Decimated triangle count: " << decimatedMesh_->polygons.size() << "\n\n";
mesh_ = decimatedMesh_;
}
}
int OdmMeshing::calcTreeDepth(size_t nPoints)
{
// Assume points are located (roughly) in a plane.
double squareSide = std::sqrt(double(nPoints));
// Calculate octree depth such that if points were equally distributed in
// a quadratic plane, there would be at least 1 point per octree node.
int depth = 0;
while(std::pow<double>(2,depth) < squareSide/2)
{
depth++;
}
return depth;
}
void OdmMeshing::writePlyFile()
{
log_ << "Saving mesh to file.\n";
if (pcl::io::savePLYFile(outputFile_.c_str(), *mesh_.get()) == -1) {
throw OdmMeshingException("Error when saving mesh to file:\n" + outputFile_ + "\n");
}
else
{
log_ << "Successfully wrote mesh to:\n"
<< outputFile_ << "\n";
}
}

1116
odm_texturing/src/OdmTexturing.cpp~
Wyświetl plik

Plik diff jest za duży Load Diff