Skip to content

Commit

Permalink
UPDATE: faster initial guess
Browse files Browse the repository at this point in the history
  • Loading branch information
@ilwoolyu
ilwoolyu committed Jun 6, 2021
1 parent 96bada3 commit 772b8a0
Show file tree
Hide file tree
Showing 6 changed files with 115 additions and 87 deletions.
4 changes: 2 additions & 2 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ project(HSD)

set(HSD_VERSION_MAJOR 1)
set(HSD_VERSION_MINOR 3)
set(HSD_VERSION_PATCH 8)
set(HSD_VERSION_PATCH 9)
set(HSD_VERSION
${HSD_VERSION_MAJOR}.${HSD_VERSION_MINOR}.${HSD_VERSION_PATCH})

Expand All @@ -30,7 +30,7 @@ if (CUDA_FOUND)
endif()

if (ENABLE_CUDA_BLAS)
set(CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS}; -O3 -gencode arch=compute_61,code=sm_61)
set(CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS}; -O3 -gencode=arch=compute_52,code=sm_52 -gencode=arch=compute_60,code=sm_60 -gencode=arch=compute_61,code=sm_61 -gencode=arch=compute_70,code=sm_70 -gencode=arch=compute_75,code=sm_75 -gencode=arch=compute_80,code=sm_80 -gencode=arch=compute_86,code=sm_86 -gencode=arch=compute_86,code=compute_86)
else()
find_package(BLAS REQUIRED)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -lopenblas")
Expand Down
173 changes: 99 additions & 74 deletions src/HSD.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@
* HSD.cpp
*
* Release: Sep 2016
* Update: July 2020
* Update: June 2021
*
* University of North Carolina at Chapel Hill
* Department of Computer Science
* Ulsan National Institute of Science and Technology
* Department of Computer Science and Engineering
*
* Ilwoo Lyu, ilwoolyu@cs.unc.edu
* Ilwoo Lyu, ilwoolyu@unist.ac.kr
*************************************************/

#include <cstring>
Expand Down Expand Up @@ -42,9 +42,10 @@ HSD::HSD(void)
m_resampling = false;
m_multi_res = true;
m_patch_range = 2;
m_guess_res = 3;
}

HSD::HSD(const char **sphere, int nSubj, const char **property, int nProperties, const char **output, const char **outputcoeff, const float *weight, int deg, const char **landmark, float weightMap, float weightLoc, float idprior, const char **coeff, const char **surf, int maxIter, const bool *fixedSubj, int icosahedron, bool realtimeCoeff, const char *tmpVariance, bool guess, const char *ico_mesh, int nCThreads, bool resampling)
HSD::HSD(const char **sphere, int nSubj, const char **property, int nProperties, const char **output, const char **outputcoeff, const float *weight, int deg, const char **landmark, float weightMap, float weightLoc, float idprior, const char **coeff, const char **surf, int maxIter, const bool *fixedSubj, int icosahedron, bool realtimeCoeff, const char *tmpVariance, bool guess, int guessRes, const char *ico_mesh, int nCThreads, bool resampling)
{
m_maxIter = maxIter;
m_nSubj = nSubj;
Expand All @@ -68,6 +69,7 @@ HSD::HSD(const char **sphere, int nSubj, const char **property, int nProperties,
m_resampling = resampling;
m_multi_res = true;
m_patch_range = 2;
m_guess_res = guessRes;
init(sphere, property, weight, landmark, weightLoc, coeff, surf, icosahedron, fixedSubj, tmpVariance, ico_mesh, nCThreads);
}

Expand Down Expand Up @@ -1273,7 +1275,7 @@ bool HSD::updateCoordinate(const float *v0, float *v1, const double *Y, const do
memcpy(v1, rv, sizeof(float) * 3);
if (delta[1] == 0 && delta[2] == 0)
return true;
else
else
return false;
}

Expand Down Expand Up @@ -2707,97 +2709,120 @@ void HSD::guessInitCoeff(void)

double delta[3];
float *feature_res = new float[nQuerySamples * (m_nProperties + m_nSurfaceProperties)];
for (double c1 = 0; c1 <= PI / 4; c1 += PI / 16)
double range_c1[3] = {0, PI / 2, PI / 8};
double range_c2[3] = {0, 2 * PI, PI / 8};
int intv = 3;

double opt_c1 = 0, opt_c2 = 0;
for (int res = 0; res < m_guess_res; res++)
{
for (double c2 = 0; c2 < 2 * PI; c2 += PI / 16)
for (double c1 = range_c1[0]; c1 <= range_c1[1]; c1 += range_c1[2])
{
if (c1 == 0 && c2 > 0) continue;
delta[1] = c1 * cos(c2);
delta[2] = c1 * sin(c2);
for (double c2 = range_c2[0]; c2 < range_c2[1]; c2 += range_c2[2])
{
if (c1 == 0 && c2 > 0) continue;
delta[1] = c1 * cos(c2);
delta[2] = c1 * sin(c2);

// exponential map (linear)
const float *axis = (P + Vector(m_spharm[subj].tan1) * delta[1] + Vector(m_spharm[subj].tan2) * delta[2]).unit().fv();
axis1[0] = axis[0]; axis1[1] = axis[1]; axis1[2] = axis[2];
// exponential map (linear)
const float *axis = (P + Vector(m_spharm[subj].tan1) * delta[1] + Vector(m_spharm[subj].tan2) * delta[2]).unit().fv();
axis1[0] = axis[0]; axis1[1] = axis[1]; axis1[2] = axis[2];

// standard pole
Vector Q = axis1;
float angle = (float)c1;
Vector A = P.cross(Q); A.unit();
// standard pole
Vector Q = axis1;
float angle = (float)c1;
Vector A = P.cross(Q); A.unit();

float rot[9];
Coordinate::rotation(A.fv(), -angle, rot);
float rot[9];
Coordinate::rotation(A.fv(), -angle, rot);

for (int i = 0; i < nQuerySamples; i++)
{
const float *v0 = &propertySamples[i * 3];
float v1[3];
Coordinate::rotPoint(v0, rot, v1);
float bary[3];
//int fid = m_spharm[subj].tree->closestFace(v1, bary);
int fid = tree->closestFace(v1, bary);
for (int k = 0; k < m_nProperties + m_nSurfaceProperties; k++)
//feature[nQuerySamples * k + i] = propertyInterpolation(&m_spharm[subj].property[nVertex * k], fid, bary, m_spharm[subj].sphere);
feature[nQuerySamples * k + i] = propertyInterpolation(&feature0[nQuerySamples * k], fid, bary, ico_mesh);
}
for (int i = 0; i < nQuerySamples; i++)
{
const float *v0 = &propertySamples[i * 3];
float v1[3];
Coordinate::rotPoint(v0, rot, v1);
float bary[3];
//int fid = m_spharm[subj].tree->closestFace(v1, bary);
int fid = tree->closestFace(v1, bary);
for (int k = 0; k < m_nProperties + m_nSurfaceProperties; k++)
//feature[nQuerySamples * k + i] = propertyInterpolation(&m_spharm[subj].property[nVertex * k], fid, bary, m_spharm[subj].sphere);
feature[nQuerySamples * k + i] = propertyInterpolation(&feature0[nQuerySamples * k], fid, bary, ico_mesh);
}

int intv = 6;
for (int w = 0; w < intv; w++)
{
for (int c3 = 0; c3 < 5; c3++)
for (int w = 0; w < intv; w++)
{
delta[0] = PI * 72.0 / 180.0 * c3 + PI * 72.0 / 180.0 * (double)w / (double)intv;
double cost = 0;
if (nLandmark > m_spharm[subj].landmark.size())
{
updateLandmark(subj);
cost += varLandmarks(subj);
}
if (nQuerySamples > 0)
for (int c3 = 0; c3 < 5; c3++)
{
if (c3 == 0)
delta[0] = 2 * PI / 5 * c3 + 2 * PI / 5 * (double)w / (double)intv;
double cost = 0;
if (nLandmark > m_spharm[subj].landmark.size())
{
Coordinate::rotation(P.fv(), (float)-delta[0], rot);
updateLandmark(subj);
cost += varLandmarks(subj);
}
for (int i = 0; i < nQuerySamples; i++)
if (nQuerySamples > 0)
{
int id = i;
if (c3 == 0)
{
float v1[3];
const float *v0 = &propertySamples[i * 3];
Coordinate::rotPoint(v0, rot, v1);
float bary[3];
int fid = tree->closestFace(v1, bary);
for (int k = 0; k < m_nProperties + m_nSurfaceProperties; k++)
feature_res[nQuerySamples * k + i] = propertyInterpolation(&feature[nQuerySamples * k], fid, bary, ico_mesh);
}
else
{
for (int j = 0; j < c3; j++)
id = cand[id];
Coordinate::rotation(P.fv(), (float)-delta[0], rot);
}
for (int k = 0; k < m_nProperties + m_nSurfaceProperties; k++)
for (int i = 0; i < nQuerySamples; i++)
{
float p = feature_res[nQuerySamples * k + id];
double m = mean[nQuerySamples * k + i];
double pm = (p - m);
cost += pm * pm / ((m_nProperties + m_nSurfaceProperties) * nQuerySamples);
int id = i;
if (c3 == 0)
{
float v1[3];
const float *v0 = &propertySamples[i * 3];
Coordinate::rotPoint(v0, rot, v1);
float bary[3];
int fid = tree->closestFace(v1, bary);
for (int k = 0; k < m_nProperties + m_nSurfaceProperties; k++)
feature_res[nQuerySamples * k + i] = propertyInterpolation(&feature[nQuerySamples * k], fid, bary, ico_mesh);
}
else
{
for (int j = 0; j < c3; j++)
id = cand[id];
}
for (int k = 0; k < m_nProperties + m_nSurfaceProperties; k++)
{
float p = feature_res[nQuerySamples * k + id];
double m = mean[nQuerySamples * k + i];
double pm = (p - m);
cost += pm * pm / ((m_nProperties + m_nSurfaceProperties) * nQuerySamples);
}
}
}
}

if (cost < mincost)
{
coeff[0] = delta[0] / m_spharm[subj].Y[0];
coeff[n] = delta[1] / m_spharm[subj].Y[0];
coeff[2 * n] = delta[2] / m_spharm[subj].Y[0];
mincost = cost;
// cout << mincost << ": ";
// cout << c1 << " " << c2 << " " << c3 << endl;
if (cost < mincost)
{
coeff[0] = delta[0] / m_spharm[subj].Y[0];
coeff[n] = delta[1] / m_spharm[subj].Y[0];
coeff[2 * n] = delta[2] / m_spharm[subj].Y[0];
mincost = cost;
// cout << mincost << ": ";
// cout << c1 << " " << c2 << " " << c3 << endl;
opt_c1 = c1;
opt_c2 = c2;
}
}
}
}
}
range_c1[0] = opt_c1 - range_c1[2];
range_c1[1] = opt_c1 + range_c1[2];
range_c1[2] /= 2;
range_c1[0] += range_c1[2];
range_c1[1] -= range_c1[2];
range_c1[0] = (range_c1[0] < 0) ? range_c1[2]: range_c1[0];
if (opt_c1 != 0)
{
range_c2[0] = opt_c2 - range_c2[2];
range_c2[1] = opt_c2 + range_c2[2];
}
range_c2[2] /= 2;
if (opt_c1 != 0) range_c2[0] += range_c2[2];
intv *= 2;
}
delete [] feature_res;
delete [] feature;
Expand Down
11 changes: 6 additions & 5 deletions src/HSD.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@
* HSD.h
*
* Release: Sep 2016
* Update: July 2020
* Update: June 2021
*
* University of North Carolina at Chapel Hill
* Department of Computer Science
* Ulsan National Institute of Science and Technology
* Department of Computer Science and Engineering
*
* Ilwoo Lyu, ilwoolyu@cs.unc.edu
* Ilwoo Lyu, ilwoolyu@unist.ac.kr
*************************************************/

#pragma once
Expand All @@ -29,7 +29,7 @@ class HSD
{
public:
HSD(void);
HSD(const char **sphere, int nSubj, const char **property, int nProperties, const char **output, const char **outputcoeff, const float *weight, int deg = 5, const char **landmark = NULL, float weightMap = 0.1, float weightLoc = 0, float idprior = 200, const char **coeff = NULL, const char **surf = NULL, int maxIter = 50, const bool *fixedSubj = NULL, int icosahedron = 7, bool realtimeCoeff = false, const char *tmpVariance = NULL, bool guess = true, const char *ico_mesh = NULL, int nCThreads = 1, bool resampling = false);
HSD(const char **sphere, int nSubj, const char **property, int nProperties, const char **output, const char **outputcoeff, const float *weight, int deg = 5, const char **landmark = NULL, float weightMap = 0.1, float weightLoc = 0, float idprior = 200, const char **coeff = NULL, const char **surf = NULL, int maxIter = 50, const bool *fixedSubj = NULL, int icosahedron = 7, bool realtimeCoeff = false, const char *tmpVariance = NULL, bool guess = true, int guessRes = 3, const char *ico_mesh = NULL, int nCThreads = 1, bool resampling = false);
~HSD(void);
void run(void);
void saveCoeff(const char *filename, int id);
Expand Down Expand Up @@ -156,6 +156,7 @@ class HSD
int m_nCThreads;
int m_patch_range;
int m_pivot_cache;
int m_guess_res;

double *m_coeff;
double *m_coeff_prev_step; // previous coefficients
Expand Down
2 changes: 1 addition & 1 deletion wrapper/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
set(CLI11_VER "1.9.0" CACHE STRING "CLI11 Version: see https://github.com/CLIUtils/CLI11/releases")
set(CLI11_VER "1.9.1" CACHE STRING "CLI11 Version: see https://github.com/CLIUtils/CLI11/releases")
mark_as_advanced(CLI11_VER)
set(CLI11_DIRS ${CMAKE_BINARY_DIR}/CMakeFiles/CLI11)
file(DOWNLOAD https://github.com/CLIUtils/CLI11/releases/download/v${CLI11_VER}/CLI11.hpp ${CLI11_DIRS}/CLI11.hpp)
Expand Down
2 changes: 2 additions & 0 deletions wrapper/PARSE_ARGS.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ bool realtimeCoeff = false;
bool noguess = false;
bool resampling = false;
int nThreads = 0;
int guessRes = 3;
std::vector<int> listFixedSubj;
int nCThreads = 0;

Expand Down Expand Up @@ -67,6 +68,7 @@ void PARSE_ARGS(int argc, char **argv)
app.add_option("--maxIter", maxIter, "Specify the maxmum number of iterations at the final phase", true)->check(CLI::NonNegativeNumber)->group("Optimization");
//app.add_option("--locationWeight", weightLoc, "Specify a weighting factor of location information", true);
app.add_option("--icomesh", icoMesh, "Specify a pre-defined icosahedral mesh, which overrides --icosaherdon")->check(CLI::ExistingFile)->group("Optimization");
app.add_option("--gres", guessRes, "Specify a search resolution for initial guess", true)->check(CLI::PositiveNumber)->group("Optimization");
app.add_flag("--noguess", noguess, "Do not execute an initial guess for rigid alignment")->group("Optimization");
app.add_flag("--resample", resampling, "Resample geometric properties using the current icosahedral level")->group("Optimization");
app.add_flag("--writecoeff", realtimeCoeff, "Write coefficient whenever the energy is minimized, which may lead to significant IO overhead")->group("Optimization");
Expand Down
10 changes: 5 additions & 5 deletions wrapper/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@
* main.cpp
*
* Release: Sep 2016
* Update: Dec 2019
* Update: June 2021
*
* University of North Carolina at Chapel Hill
* Department of Computer Science
* Ulsan National Institute of Science and Technology
* Department of Computer Science and Engineering
*
* Ilwoo Lyu, ilwoolyu@cs.unc.edu
* Ilwoo Lyu, ilwoolyu@unist.ac.kr
*************************************************/

#include <cstdlib>
Expand Down Expand Up @@ -198,7 +198,7 @@ int main(int argc, char **argv)
openblas_set_num_threads(nThreads);
#endif

HSD hsd(sphere, nSubj, property, nProperties / nSubj, output, outputcoeff, weight, degree, landmark, weightMap, weightLoc, idprior, coeff, surf, maxIter, fixedSubj, icosa, realtimeCoeff, prior, !noguess, icomesh, nCThreads, resampling);
HSD hsd(sphere, nSubj, property, nProperties / nSubj, output, outputcoeff, weight, degree, landmark, weightMap, weightLoc, idprior, coeff, surf, maxIter, fixedSubj, icosa, realtimeCoeff, prior, !noguess, guessRes, icomesh, nCThreads, resampling);
hsd.run();

// delete memory allocation
Expand Down

0 comments on commit 772b8a0

Please sign in to comment.