https://github.com/cran/RandomFields
Tip revision: 6b74edb41c4ebbda06af16d7834cba9b7f23e819 authored by Martin Schlather on 27 October 2011, 00:00:00 UTC
version 2.0.52
version 2.0.52
Tip revision: 6b74edb
kleinn.cc
/*
Authors
Martin Schlather, martin.schlather@math.uni-goettingen.de
Definition of correlation functions and derivatives (spectral measures,
tbm operators)
Note:
* Never use the below functions directly, but only by the functions indicated
in RFsimu.h, since there is no error check (e.g. initialization of RANDOM)
* VARIANCE, SCALE are not used here
* definitions for the random coin method can be found in MPPFcts.cc
* definitions for genuinely anisotropic or nonstationary models are in
SophisticatedModel.cc; hyper models also in Hypermodel.cc
Copyright (C) 2001 -- 2003 Martin Schlather
Copyright (C) 2004 -- 2004 Yindeng Jiang & Martin Schlather
Copyright (C) 2005 -- 2011 Martin Schlather
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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
#include <math.h>
#include "RF.h"
#ifdef DoublePrecision
#define boolean int
#define uchar int
#define uword int
#else
#define boolean bool
#define uchar unsigned char
#define uword short unsigned int
#endif
//#define gz unsigned short int
void threshold(int n, res_type *in, double percent, boolean* out) {
int i;
double threshold;
threshold = qnorm(percent, 0.0, 1.0, 1, 0);
for (i=0; i<n; i++) {
out[i] = in[i] < threshold;
}
}
void coarsening(int nr, int nc, boolean *in, int res, uchar *out) {
if (res > 15 && sizeof(uchar)==1) error("max coarse is 15");
int threshold = (res * res) / 2,
ncDres = nc / res,
nrres = nr * res,
nrDres = nr / res;
boolean *pin, *pc, *pinbase, *lastres,
*lastin = in + ncDres * nrres;
uchar *proutende, *pr,
*pout = out,
*lastout = out + nrDres * ncDres;
for(pinbase = in; pinbase<lastin; pout += nrDres) {
pc = pinbase + nrres;
for (; pinbase<pc; pinbase+=nr) {
pin = pinbase;
proutende = pout + nrDres;
for (pr = pout; pr<proutende; pr++) {
for (lastres = pin + res; pin<lastres; pin++) {
*pr += *pin; // *pr = *pr + *pin;
}
}
}
}
for (pout = out; pout<lastout; pout++) {
*pout = *pout > threshold;
}
}
void reference_area(int nr, int nc, uchar *in, int edge, uword *out) {
int i, j, old, laststart,
n = nc * nr,
edgeM1 = edge - 1,
edgeM1nr = nr * edgeM1,
last = n - edge,
edgenr = edge * nr,
nrMedge1 = nr - edge +1;
uword dummy;
out[0] = 0.0;
for(i=0; i<edge; i++) out[0] += in[i];
for(i=1; i<=last; i++) {
out[i] = out[i-1] - in[i-1] + in[i + edgeM1];
}
last = n - 1 - edgeM1nr;
for(i=0; i<nrMedge1; i++) {
old = out[i];
for(j=i + nr, laststart = i + edgenr; j<laststart; j+=nr) {
out[i] += out[j];
}
for(j=i + nr; j<=last; j+=nr) {
dummy = out[j];
out[j] = out[j - nr] - old + out[j + edgeM1nr];
old = dummy;
}
}
}
void analyse_internal(int* keynr,
double *gausspercent, int *ngaussp,
int *coarse, int *ncrs,
int *edges, int *nedg,
double *percent, int *nperc,
res_type *image,
boolean *binary,
uchar * decreased,
int *Nrdecr,
int *Ncdecr,
uword *refarea,
int *Areathreshold,
double *result
) {
int k, Err, nettoedge, igaussp, ires, iedg, iperc,
lastnc, forest, land, j, i, lastnr,
nrdecr=-1,
ncdecr=-1,
areathreshold=-1,
n = KEY[*keynr].loc.totalpoints,
nr = KEY[*keynr].loc.length[0],
nc = KEY[*keynr].loc.length[1],
one = 1,
zero = 0;
double percent_covered;
bool debug = false;
DoSimulateRF(keynr, &one, &zero, image, &Err);
if (Err != 0) {
error("DoSimulateRF does not work");
}
if (debug) PRINTF("simulation done with %dx%d pixels\n", nr, nc);
k = 0;
for (igaussp = 0; igaussp < *ngaussp; igaussp++) {
threshold(n, image, gausspercent[igaussp], binary);
for (ires = 0; ires<*ncrs; ires++) {
nrdecr = nr / coarse[ires];
ncdecr = nc / coarse[ires];
if (debug) PRINTF("thresh=%2.2f coarse=%d new size:%dx%d pixels\n",
gausspercent[igaussp], coarse[ires],
nrdecr, ncdecr);
if (coarse[ires] == 1) {
for (i=0; i<n; i++) decreased[i] = binary[i];
} else {
for (i=0; i<n; i++) decreased[i] = 0;
coarsening(nr, nc, binary, coarse[ires], decreased);
}
for (iedg = 0; iedg<*nedg; iedg++) {
nettoedge = (int) round((double) edges[iedg] /
(double) coarse[ires]);
if (nettoedge > nrdecr) nettoedge = nrdecr;
if (nettoedge > ncdecr) nettoedge = ncdecr;
if (debug) PRINTF(" netto.referenz=%d\n", nettoedge);
for (iperc=0; iperc<*nperc; iperc++) {
areathreshold =
(int) ceil(nettoedge * nettoedge * percent[iperc]);
if (debug) PRINTF(" min=%2.0f%(=%dpx) ",
percent[iperc] * 100, areathreshold
);
reference_area(nrdecr, ncdecr, decreased, nettoedge,
refarea);
lastnc = (ncdecr - nettoedge +1) * nrdecr;
forest = 0;
land = (nrdecr - nettoedge + 1) * (ncdecr - nettoedge + 1);
for (j=0; j<lastnc; j+=nrdecr) {
lastnr = j + nrdecr - nettoedge;
for (i=j; i<=lastnr; i++) {
forest += (refarea[i] >= areathreshold);
}
}
percent_covered = (double) forest / (double) land;
if (debug) PRINTF("-> forest=%d land=(%d-%d+1)^2=%d coverd=%2.2f%\n",
forest, nrdecr, nettoedge, land, percent_covered);
result[k++] = percent_covered;
}
}
}
}
*Nrdecr = nrdecr;
*Ncdecr = ncdecr;
*Areathreshold = areathreshold;
}
void analyseForst(int* keynr,
double *gausspercent, int *ngaussp,
int *coarse, int *ncrs,
int *edges, int *nedg,
double *percent, int *nperc,
double *result) {
int nrdecr, ncdecr, areathreshold,
n = KEY[*keynr].loc.totalpoints;
res_type *image;
boolean *binary;
uchar *decreased;
uword *refarea;
PRINTF("memory needed %e + %e + %e + %e\n",
(double) n * sizeof(res_type),
(double) n * sizeof(boolean),
(double) n * sizeof(uchar),
(double) n * sizeof(uword)
);
if ((image = (res_type*) malloc(n * sizeof(res_type))) == NULL)
error("mem err 1");
if ((binary = (boolean*) malloc(n * sizeof(boolean))) == NULL)
error("mem err 2");
if ((decreased = (uchar*) calloc(n, sizeof(uchar))) == NULL)
error("mem err 3");
if ((refarea = (uword*) malloc(n * sizeof(uword))) == NULL)
error("mem err 4");
analyse_internal(keynr, gausspercent, ngaussp, coarse, ncrs,
edges, nedg, percent, nperc,
image, binary, decreased,
&nrdecr, &ncdecr, refarea,
&areathreshold,
result);
free(refarea);
free(decreased);
free(binary);
free(image);
}
void analyseForstImages(int* keynr,
double *gausspercent, int *ngaussp,
int *coarse, int *ncrs,
int *edges, int *nedg,
double *percent, int *nperc,
//
double *image,
int *binary,
int *decreased,
int *nrdecr,
int *ncdecr,
int *refarea,
//
int *areathreshold,
double *result
) {
if (sizeof(res_type) != sizeof(double)) {
error("images cannot be returned (res_type)\n");
}
if (sizeof(boolean) != sizeof(int)) {
error("images cannot be returned (bool)\n");
}
if (sizeof(uchar) != sizeof(int)) {
error("images cannot be returned (uchar)\n");
}
if (sizeof(uword) != sizeof(int)) {
error("images cannot be returned (uword)\n");
}
analyse_internal(keynr, gausspercent, ngaussp, coarse, ncrs,
edges, nedg, percent, nperc,
(res_type *)image,
(boolean*) binary,
(uchar*) decreased,
nrdecr, ncdecr,
(uword*) refarea, areathreshold,
result);
// assert(*areathreshold < 10);
}