/*
* Written 2003 Lukas Kunc <Lukas.Kunc@seznam.cz>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*/
#include <stdlib.h>
#include "xmalloc.h"
#include "lambda.h"
static void get_variance_mirror(image_t* variance, image_t* img, real_t* pmin, real_t* pmax, int winsize)
{
int i, j, k, l;
real_t sum, sum2;
real_t c;
real_t num_points;
real_t var, minvar, maxvar;
minvar = R(1e20);
maxvar = R(0.0);
num_points = (real_t)((2*winsize+1)*(2*winsize+1));
for (i = 0; i < variance->x; i++) {
for (j = 0; j < variance->y; j++) {
sum = R(0.0);
sum2 = R(0.0);
for (k = -winsize; k <= winsize; k++) {
for (l = -winsize; l <= winsize; l++) {
c = image_get_mirror(img, i+k, j+l);
sum += c;
sum2 += c*c;
}
}
sum /= num_points;
sum2 /= num_points;
sum *= sum;
var = sum2 - sum;
image_set(variance, i, j, var);
if (var > maxvar) {
maxvar = var;
}
else if (var < minvar) {
minvar = var;
}
}
}
*pmax = maxvar;
*pmin = minvar;
}
static void get_variance_period(image_t* variance, image_t* img, real_t* pmin, real_t* pmax, int winsize)
{
int i, j, k, l;
real_t sum, sum2;
real_t c;
real_t num_points;
real_t var, minvar, maxvar;
minvar = R(1e20);
maxvar = R(0.0);
num_points = (real_t)((2*winsize+1)*(2*winsize+1));
for (i = 0; i < variance->x; i++) {
for (j = 0; j < variance->y; j++) {
sum = R(0.0);
sum2 = R(0.0);
for (k = -winsize; k <= winsize; k++) {
for (l = -winsize; l <= winsize; l++) {
c = image_get_period(img, i+k, j+l);
sum += c;
sum2 += c*c;
}
}
sum /= num_points;
sum2 /= num_points;
sum *= sum;
var = sum2 - sum;
image_set(variance, i, j, var);
if (var > maxvar) {
maxvar = var;
}
else if (var < minvar) {
minvar = var;
}
}
}
*pmax = maxvar;
*pmin = minvar;
}
lambda_t* lambda_create(lambda_t* lambda, int x, int y, real_t minlambda, int winsize, convmask_t* filter)
{
lambda->x = x;
lambda->y = y;
lambda->minlambda = minlambda;
lambda->winsize = winsize;
lambda->filter = filter;
lambda->lambda = (real_t*)xmalloc(sizeof(real_t) * x * y);
return lambda;
}
void lambda_destroy(lambda_t* lambda)
{
xfree(lambda->lambda);
}
void lambda_set_mirror(lambda_t* lambda, int mirror)
{
lambda->mirror = mirror;
}
void lambda_set_nl(lambda_t* lambda, int nl)
{
lambda->nl = nl;
}
lambda_t* lambda_calculate_period(lambda_t* lambda, image_t* image)
{
image_t imgenh, *imgcal;
image_t variance;
real_t minvar, maxvar;
real_t akoef, bkoef;
int i, size;
if (lambda->filter)
{
imgcal = image_create_copyparam(&imgenh, image);
image_convolve_period(imgcal, image, lambda->filter);
}
else
{
imgcal = image;
}
image_create_copyparam(&variance, imgcal);
get_variance_period(&variance, imgcal, &minvar, &maxvar, lambda->winsize);
bkoef = (R(1.0) - lambda->minlambda)/(minvar - maxvar);
akoef = R(1.0) - (minvar*(R(1.0) - lambda->minlambda))/(minvar - maxvar);
size = lambda->x * lambda->y;
for (i = 0; i < size; i++) {
lambda->lambda[i] = akoef + bkoef*variance.data[i];
}
if (lambda->filter)
{
image_destroy(imgcal);
}
image_destroy(&variance);
return lambda;
}
lambda_t* lambda_calculate_period_nl(lambda_t* lambda, image_t* image)
{
image_t imgenh, *imgcal;
image_t variance;
real_t minvar, maxvar;
real_t alpha;
int i, size;
if (lambda->filter)
{
imgcal = image_create_copyparam(&imgenh, image);
image_convolve_period(imgcal, image, lambda->filter);
}
else
{
imgcal = image;
}
image_create_copyparam(&variance, imgcal);
get_variance_period(&variance, imgcal, &minvar, &maxvar, lambda->winsize);
alpha = (R(1.0)-lambda->minlambda)/(lambda->minlambda*(maxvar-minvar));
size = lambda->x * lambda->y;
for (i = 0; i < size; i++) {
lambda->lambda[i] = R(1.0)/(R(1.0)+alpha*(variance.data[i]-minvar));
}
if (lambda->filter)
{
image_destroy(imgcal);
}
image_destroy(&variance);
return lambda;
}
lambda_t* lambda_calculate_mirror(lambda_t* lambda, image_t* image)
{
image_t imgenh, *imgcal;
image_t variance;
real_t minvar, maxvar;
real_t akoef, bkoef;
int i, size;
if (lambda->filter)
{
imgcal = image_create_copyparam(&imgenh, image);
image_convolve_mirror(imgcal, image, lambda->filter);
}
else
{
imgcal = image;
}
image_create_copyparam(&variance, imgcal);
get_variance_mirror(&variance, imgcal, &minvar, &maxvar, lambda->winsize);
bkoef = (R(1.0) - lambda->minlambda)/(minvar - maxvar);
akoef = R(1.0) - (minvar*(R(1.0) - lambda->minlambda))/(minvar - maxvar);
size = lambda->x * lambda->y;
for (i = 0; i < size; i++) {
lambda->lambda[i] = akoef + bkoef*variance.data[i];
}
if (lambda->filter)
{
image_destroy(imgcal);
}
image_destroy(&variance);
return lambda;
}
lambda_t* lambda_calculate_mirror_nl(lambda_t* lambda, image_t* image)
{
image_t imgenh, *imgcal;
image_t variance;
real_t minvar, maxvar;
real_t alpha;
int i, size;
if (lambda->filter)
{
imgcal = image_create_copyparam(&imgenh, image);
image_convolve_mirror(imgcal, image, lambda->filter);
}
else
{
imgcal = image;
}
image_create_copyparam(&variance, imgcal);
get_variance_mirror(&variance, imgcal, &minvar, &maxvar, lambda->winsize);
alpha = (R(1.0)-lambda->minlambda)/(lambda->minlambda*(maxvar-minvar));
size = lambda->x * lambda->y;
for (i = 0; i < size; i++) {
lambda->lambda[i] = R(1.0)/(R(1.0)+alpha*(variance.data[i]-minvar));
}
if (lambda->filter)
{
image_destroy(imgcal);
}
image_destroy(&variance);
return lambda;
}
lambda_t* lambda_calculate(lambda_t* lambda, image_t* image)
{
if (lambda->mirror)
{
if (lambda->nl) return lambda_calculate_mirror_nl(lambda, image);
else return lambda_calculate_mirror(lambda, image);
}
else
{
if (lambda->nl) return lambda_calculate_period_nl(lambda, image);
else return lambda_calculate_period(lambda, image);
}
}
/*
* GET / SET
*/
#if !defined(INLINE) && !defined(INLINE_MACRO)
real_t lambda_get_mirror(lambda_t* lambda, int x, int y)
{
return MACRO_LAMBDA_GET_MIRROR(lambda, x, y);
}
real_t lambda_get_period(lambda_t* lambda, int x, int y)
{
return MACRO_LAMBDA_GET_PERIOD(lambda, x, y);
}
#endif
|