Skip to content

fix(gabor): keep filter-response magnitudes real-valued, then vet GABOR vs scikit-image #393

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
darkclad wants to merge 1 commit into
base: main
Choose a base branch
from
Open

fix(gabor): keep filter-response magnitudes real-valued, then vet GABOR vs scikit-image #393

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/nyx/cache.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@
GpuCache <cufftDoubleComplex>& gabor_linear_image, // gabor
GpuCache <cufftDoubleComplex>& gabor_result,
GpuCache <cufftDoubleComplex>& gabor_linear_kernel,
GpuCache <PixIntens>& gabor_energy_image,
GpuCache <RealPixIntens>& gabor_energy_image,
// in
bool needContour,
bool needErosion,
Expand Down Expand Up @@ -356,7 +356,7 @@
GpuCache <cufftDoubleComplex>& gabor_linear_image,
GpuCache <cufftDoubleComplex>& gabor_result,
GpuCache <cufftDoubleComplex>& gabor_linear_kernel,
GpuCache <PixIntens>& gabor_energy_image)
GpuCache <RealPixIntens>& gabor_energy_image)
{
// clouds

Expand Down
6 changes: 3 additions & 3 deletions src/nyx/cache.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ class GpusideCache
/*--extern--*/ GpuCache <cufftDoubleComplex> gabor_linear_image; // (img_plus_ker_size* n_filters); // ROI + kernel image
/*--extern--*/ GpuCache <cufftDoubleComplex> gabor_result; // (img_plus_ker_size* n_filters);
/*--extern--*/ GpuCache <cufftDoubleComplex> gabor_linear_kernel; // (img_plus_ker_size* n_filters);
/*--extern--*/ GpuCache <PixIntens> gabor_energy_image; // (img_plus_ker_size* n_filters);
/*--extern--*/ GpuCache <RealPixIntens> gabor_energy_image; // real-valued (was PixIntens, which truncated the Gabor response to unsigned int); (img_plus_ker_size* n_filters);

//
// these need to be called after "prescan" (phase 0)
Expand Down Expand Up @@ -166,7 +166,7 @@ class GpusideCache
GpuCache <cufftDoubleComplex>& gabor_linear_image,
GpuCache <cufftDoubleComplex>& gabor_result,
GpuCache <cufftDoubleComplex>& gabor_linear_kernel,
GpuCache <PixIntens>& gabor_energy,
GpuCache <RealPixIntens>& gabor_energy,
// in
bool needContour,
bool needErosion,
Expand All @@ -193,7 +193,7 @@ class GpusideCache
GpuCache <cufftDoubleComplex>& gabor_linear_image,
GpuCache <cufftDoubleComplex>& gabor_result,
GpuCache <cufftDoubleComplex>& gabor_linear_kernel,
GpuCache <PixIntens>& gabor_energy_image
GpuCache <RealPixIntens>& gabor_energy_image
);

// these need to be called in "reduce_trivial"
Expand Down
135 changes: 67 additions & 68 deletions src/nyx/features/gabor.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,10 +61,9 @@ void GaborFeature::calculate (LR& r, const Fsettings& s)

// Temp buffers

// --1
ImageMatrix e2img;
e2img.allocate (Im0.width, Im0.height);
readOnlyPixels pix_plane = e2img.ReadablePixels();
// --1 Real-valued filter-response image (was a PixIntens ImageMatrix, which truncated the
// response magnitudes to unsigned int and floored sub-integer responses to 0).
std::vector<double> e2 (Im0.width * Im0.height);

// --2
std::vector<double> auxC ((Im0.width + n - 1) * (Im0.height + n - 1) * 2);
Expand All @@ -77,13 +76,15 @@ void GaborFeature::calculate (LR& r, const Fsettings& s)
ty.resize (n + 1);

// Compute the baseline score before applying high-pass Gabor filters
GaborEnergy (Im0, e2img.writable_data_ptr(), auxC.data(), auxG.data(), f0LP, sig2lam, gamma, M_PI_2, n); // compromise pi/2 theta
GaborEnergy (Im0, e2.data(), auxC.data(), auxG.data(), f0LP, sig2lam, gamma, M_PI_2, n); // compromise pi/2 theta

// Values that we need for scoring filter responses
Moments2 local_stats;
e2img.GetStats(local_stats);
double maxval = local_stats.max__(),
cmpval = local_stats.min__();
double maxval = e2[0], cmpval = e2[0];
for (double a : e2)
{
maxval = std::max(maxval, a);
cmpval = std::min(cmpval, a);
}

// intercept blank baseline filter response
if (maxval == cmpval)
Expand All @@ -93,29 +94,27 @@ void GaborFeature::calculate (LR& r, const Fsettings& s)
return;
}

// Score the baseline signal
// Score the baseline signal
unsigned long baselineScore = 0;
for (auto a : pix_plane)
if (double(a) > cmpval)
for (double a : e2)
if (a > cmpval)
baselineScore++;

// Iterate frequencies and score corresponding filter response over the baseline
for (int i=0; i < nFreqs; i++)
{
// filter response for i-th frequency
writeablePixels e2_pix_plane = e2img.WriteablePixels();

// -- unpack a frequency-angle pair
const auto& ft = f0_theta_pairs[i];
auto f0 = ft.first;
auto theta = ft.second;

GaborEnergy (Im0, e2_pix_plane.data(), auxC.data(), auxG.data(), f0, sig2lam, gamma, theta, n);
// filter response for i-th frequency (overwrites e2)
GaborEnergy (Im0, e2.data(), auxC.data(), auxG.data(), f0, sig2lam, gamma, theta, n);

// score it
unsigned long afterGaborScore = 0;
for (auto a : e2_pix_plane)
if (double(a)/maxval > GRAYthr)
for (double a : e2)
if (a/maxval > GRAYthr)
afterGaborScore++;

// save the score as feature value
Expand Down Expand Up @@ -152,10 +151,9 @@ void GaborFeature::calculate_gpu (LR& r)

// Temp buffers

// --1
ImageMatrix e2img;
e2img.allocate (Im0.width, Im0.height);
readOnlyPixels pix_plane = e2img.ReadablePixels();
// --1 Real-valued filter-response image (was a PixIntens ImageMatrix, which truncated the
// response magnitudes to unsigned int).
std::vector<double> e2 (Im0.width * Im0.height);

// --2
std::vector<double> auxC ((Im0.width + n - 1) * (Im0.height + n - 1) * 2);
Expand All @@ -168,37 +166,37 @@ void GaborFeature::calculate_gpu (LR& r)
ty.resize (n + 1);

// Compute the baseline score before applying high-pass Gabor filters
GaborEnergyGPU (Im0, e2img.writable_data_ptr(), auxC.data(), auxG.data(), f0LP, sig2lam, gamma, M_PI_2, n); // compromise pi/2 theta
GaborEnergyGPU (Im0, e2.data(), auxC.data(), auxG.data(), f0LP, sig2lam, gamma, M_PI_2, n); // compromise pi/2 theta

// Values that we need for scoring filter responses
Moments2 local_stats;
e2img.GetStats(local_stats);
double maxval = local_stats.max__(),
cmpval = local_stats.min__();
double maxval = e2[0], cmpval = e2[0];
for (double a : e2)
{
maxval = std::max(maxval, a);
cmpval = std::min(cmpval, a);
}

// Score the baseline signal
// Score the baseline signal
unsigned long baselineScore = 0;
for (auto a : pix_plane)
if (double(a) > cmpval)
for (double a : e2)
if (a > cmpval)
baselineScore++;

// Iterate frequencies and score corresponding filter response over the baseline
for (int i=0; i < nFreqs; i++)
{
// filter response for i-th frequency
writeablePixels e2_pix_plane = e2img.WriteablePixels();

// -- unpack a frequency-angle pair
const auto& ft = f0_theta_pairs[i];
auto f0 = ft.first;
auto theta = ft.second;

GaborEnergyGPU (Im0, e2_pix_plane.data(), auxC.data(), auxG.data(), f0, sig2lam, gamma, theta, n);
// filter response for i-th frequency (overwrites e2)
GaborEnergyGPU (Im0, e2.data(), auxC.data(), auxG.data(), f0, sig2lam, gamma, theta, n);

// score it
unsigned long afterGaborScore = 0;
for (auto a : e2_pix_plane)
if (double(a)/maxval > GRAYthr)
for (double a : e2)
if (a/maxval > GRAYthr)
afterGaborScore++;

// save the score as feature value
Expand Down Expand Up @@ -262,8 +260,9 @@ void GaborFeature::calculate_gpu_multi_filter (LR & r, size_t roiidx, GpusideCac
thetas.push_back(ft.second);
}

// Compute the baseline score before applying high-pass Gabor filters
std::vector<std::vector<PixIntens>> responses (num_filters, std::vector<PixIntens>(Im0.width * Im0.height));
// Compute the baseline score before applying high-pass Gabor filters.
// Real-valued responses (was PixIntens, which truncated the magnitudes to unsigned int).
std::vector<std::vector<double>> responses (num_filters, std::vector<double>(Im0.width * Im0.height));

// Calculate low-passed baseline and high-passed filter responses
// 'responses' is montage of Gabor energy of filter responses
Expand All @@ -272,8 +271,8 @@ void GaborFeature::calculate_gpu_multi_filter (LR & r, size_t roiidx, GpusideCac
// Examine the baseline signal

// we need to get these 3 values from response[0], the baseline signal
PixIntens maxval = 0,
cmpval = UINT16_MAX; // min
double maxval = 0.0, // magnitudes are non-negative
cmpval = std::numeric_limits<double>::max(); // min
size_t baselineScore = 0;

size_t wh = Im0.width * Im0.height;
Expand Down Expand Up @@ -451,15 +450,15 @@ void GaborFeature::Gabor (double* Gex, double f0, double sig2lam, double gamma,

// Computes Gabor energy
void GaborFeature::GaborEnergy (
const ImageMatrix& Im,
PixIntens* out,
double* auxC,
const ImageMatrix& Im,
double* out,
double* auxC,
double* Gexp,
double f0,
double sig2lam,
double gamma,
double theta,
int n)
double f0,
double sig2lam,
double gamma,
double theta,
int n)
{
int n_gab = n;

Expand Down Expand Up @@ -491,19 +490,19 @@ void GaborFeature::GaborEnergy (
conv_dud (auxC, pix_plane.data(), Gexp, Im.width, Im.height, n_gab, n_gab);

decltype(Im.height) b = 0;
for (auto y = (int)ceil((double)n / 2); b < Im.height; y++)
for (auto y = (int)ceil((double)n / 2); b < Im.height; y++)
{
decltype(Im.width) a = 0;
for (auto x = (int)ceil((double)n / 2); a < Im.width; x++)
for (auto x = (int)ceil((double)n / 2); a < Im.width; x++)
{
if (std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2]) || std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1]))
if (std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2]) || std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1]))
{
out[b * Im.width + a] = (PixIntens) std::numeric_limits<double>::quiet_NaN();
out[b * Im.width + a] = std::numeric_limits<double>::quiet_NaN(); // keep response real-valued (was truncated to unsigned int)
a++;
continue;
}

out[b * Im.width + a] = (PixIntens) sqrt(pow(auxC[y * 2 * (Im.width + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1], 2));
out[b * Im.width + a] = sqrt(pow(auxC[y * 2 * (Im.width + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1], 2)); // real magnitude, no unsigned-int truncation
a++;
}
b++;
Expand All @@ -512,15 +511,15 @@ void GaborFeature::GaborEnergy (

#ifdef USE_GPU
void GaborFeature::GaborEnergyGPU (
const ImageMatrix& Im,
PixIntens* /* double* */ out,
double* auxC,
const ImageMatrix& Im,
double* out,
double* auxC,
double* Gexp,
double f0,
double sig2lam,
double gamma,
double theta,
int n)
double f0,
double sig2lam,
double gamma,
double theta,
int n)
{
int n_gab = n;

Expand All @@ -536,28 +535,28 @@ void GaborFeature::GaborEnergyGPU (
}

decltype(Im.height) b = 0;
for (auto y = (int)ceil((double)n / 2); b < Im.height; y++)
for (auto y = (int)ceil((double)n / 2); b < Im.height; y++)
{
decltype(Im.width) a = 0;
for (auto x = (int)ceil((double)n / 2); a < Im.width; x++)
for (auto x = (int)ceil((double)n / 2); a < Im.width; x++)
{
if (std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2]) || std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1]))
if (std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2]) || std::isnan(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1]))
{
out[b * Im.width + a] = (PixIntens) std::numeric_limits<double>::quiet_NaN();
out[b * Im.width + a] = std::numeric_limits<double>::quiet_NaN(); // keep response real-valued (was truncated to unsigned int)
a++;
continue;
}

out[b * Im.width + a] = (PixIntens) sqrt(pow(auxC[y * 2 * (Im.width + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1], 2));
out[b * Im.width + a] = sqrt(pow(auxC[y * 2 * (Im.width + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (Im.width + n - 1) + x * 2 + 1], 2)); // real magnitude, no unsigned-int truncation
a++;
}
b++;
}
}

void GaborFeature::GaborEnergyGPUMultiFilter (
const ImageMatrix& Im,
std::vector<std::vector<PixIntens>>& out, // energy image
const ImageMatrix& Im,
std::vector<std::vector<double>>& out, // energy image (real-valued; was PixIntens)
double* auxC, // batch of filter responses in complex layout
double* Gexp,
const std::vector<double>& f0s, // f0-s matching 'thetas'
Expand Down
37 changes: 19 additions & 18 deletions src/nyx/features/gabor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,33 +92,34 @@ class GaborFeature: public FeatureMethod
// Buffers for Gabor amplitudes. Used by method Gabor()
std::vector<double> tx, ty;

// Computes Gabor energy
// Computes Gabor energy. 'out' holds real-valued filter-response magnitudes, kept as double
// to avoid the precision loss of the former PixIntens (unsigned int) truncation.
void GaborEnergy (
const ImageMatrix& Im,
PixIntens* out,
double* auxC,
double* Gex,
double f0,
double sig2lam,
double gamma,
double theta,
const ImageMatrix& Im,
double* out,
double* auxC,
double* Gex,
double f0,
double sig2lam,
double gamma,
double theta,
int n);

#ifdef USE_GPU
void GaborEnergyGPU (
const ImageMatrix& Im,
PixIntens* /* double* */ out,
double* auxC,
double* Gex,
double f0,
double sig2lam,
double gamma,
double theta,
const ImageMatrix& Im,
double* out,
double* auxC,
double* Gex,
double f0,
double sig2lam,
double gamma,
double theta,
int n);

void GaborEnergyGPUMultiFilter(
const ImageMatrix& Im,
std::vector<std::vector<PixIntens>>& out,
std::vector<std::vector<double>>& out,
double* auxC,
double* Gexp,
const std::vector<double>& f0, // frequencies matching 'thetas'
Expand Down
8 changes: 6 additions & 2 deletions src/nyx/features/gabor_nontriv.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,10 @@ void GaborFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader&)

double cmp_a = max0 * GRAYthr;

// NOTE (pre-existing, separate from the response-truncation fix): originalScore is the baseline
// score and must be the count of baseline pixels above the baseline minimum (as the trivial
// calculate() computes it), but it is left at 0 here, so every fvals[i] below divides by ~0.
// The out-of-core Gabor path needs the baseline min + count wired in to match the in-RAM path.
size_t originalScore = 0;
for (int i=0; i < nF; i++)
{
Expand Down Expand Up @@ -128,7 +132,7 @@ void GaborFeature::GaborEnergy_NT2 (
for (auto y = xy0, b = 0; b < winY; y++, b++)
for (auto x = xy0, a = 0; a < winX; x++, a++)
{
auto inten = (PixIntens)sqrt(pow(auxC[y * 2 * (winX + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (winX + n - 1) + x * 2 + 1], 2));
double inten = sqrt(pow(auxC[y * 2 * (winX + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (winX + n - 1) + x * 2 + 1], 2)); // real magnitude, no unsigned-int truncation
// Calculate requested summary
if (max_or_threshold)
max_val = std::max(max_val, double(inten));
Expand Down Expand Up @@ -163,7 +167,7 @@ void GaborFeature::GaborEnergy_NT2 (
for (auto y = xy0, b = 0; b < height; y++, b++)
for (auto x = xy0, a = 0; a < width; x++, a++)
{
auto inten = (PixIntens)sqrt(pow(auxC[y * 2 * (winX + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (winX + n - 1) + x * 2 + 1], 2));
double inten = sqrt(pow(auxC[y * 2 * (winX + n - 1) + x * 2], 2) + pow(auxC[y * 2 * (winX + n - 1) + x * 2 + 1], 2)); // real magnitude, no unsigned-int truncation
// Calculate requested summary
if (max_or_threshold)
max_val = std::max(max_val, double(inten));
Expand Down
Loading
Loading