diff --git a/R/mimetic.R b/R/mimetic.R
index c428df8dfcdc9fe47c9d30d06bd54e35d372587d..ab0c0048dace8b2a49917f8e9597efa3a1d82a45 100755
--- a/R/mimetic.R
+++ b/R/mimetic.R
@@ -2,7 +2,7 @@
#RP 11/06/2013
###################################################
-mimetic<-function(x,upto=NULL,by=NULL,prec=NULL,nsimax=3000,conv=50) {
+mimetic<-function(x,upto=NULL,by=NULL,prec=NULL,nsimax=3000,conv=50,mask=NULL) {
# checking for input parameters
stopifnot(inherits(x,"fads")||inherits(x,"spp"))
if(inherits(x,"fads")) {
@@ -50,13 +50,35 @@ mimetic<-function(x,upto=NULL,by=NULL,prec=NULL,nsimax=3000,conv=50) {
PACKAGE="ads")
}
else {
+ if (!missing(mask) && !is.null(mask)) {
+ stopifnot(inherits(mask, "RasterStack")||inherits(mask, "RasterLayer"))
+ if (inherits(mask, "RasterStack"))
+ mask <- mask[[1]]
+ mask.coord <- coordinates(mask)[!is.na(values(mask)),]
+ mask.n <- nrow(mask.coord)
+ mask.x <- mask.coord[,1]
+ mask.y <- mask.coord[,2]
+ mask.dx <- abs(mask@extent@xmax-mask@extent@xmin)/mask@ncols
+ mask.dy <- abs(mask@extent@ymax-mask@extent@ymin)/mask@nrows
+ mask.spp <- spp(x=mask.coord[,1],y=mask.coord[,2],window=p$window)
+ res<-.C("mimetic_rect",
+ as.integer(p$n),as.double(p$x),as.double(p$y),as.double(surface),
+ as.double(xmin),as.double(xmax),as.double(ymin),as.double(ymax),
+ as.double(mask.spp$n),as.double(mask.spp$x), as.double(mask.spp$y),as.double(mask.dx),as.double(mask.dy),
+ as.double(prec),as.integer(tmax),as.double(by),
+ as.double(lobs),as.integer(nsimax),as.integer(conv),cost=double(nsimax),
+ g=double(tmax),k=double(tmax),xx=double(p$n),yy=double(p$n),mess=as.integer(1),
+ PACKAGE="ads")
+ } else {
res<-.C("mimetic_rect",
as.integer(p$n),as.double(p$x),as.double(p$y),as.double(surface),
as.double(xmin),as.double(xmax),as.double(ymin),as.double(ymax),
+ as.integer(0),as.double(0), as.double(0),as.double(0),as.double(0),
as.double(prec),as.integer(tmax),as.double(by),
as.double(lobs),as.integer(nsimax),as.integer(conv),cost=double(nsimax),
g=double(tmax),k=double(tmax),xx=double(p$n),yy=double(p$n),mess=as.integer(1),
PACKAGE="ads")
+ }
}
}
else if("circle"%in%p$window$type) {
diff --git a/src/Zlibs.c b/src/Zlibs.c
index 362f14b74b01523e759b57696b889f55d18dde3b..e089f9169c5cb4e33b936741c54fc5c638ecee0b 100755
--- a/src/Zlibs.c
+++ b/src/Zlibs.c
@@ -2572,7 +2572,7 @@ int intertype_rect_ic(int *point_nb1, double *x1, double *y1, int *point_nb2, do
r=0;
while (erreur!=0)
{
- erreur=mimetic_rect(point_nb1, x1, y1, surface, xmi, xma, ymi, yma, prec, t2, dt, lt, nsimax, conv, cost, gt, kt, x, y, 0);
+ erreur=mimetic_rect(point_nb1, x1, y1, surface, xmi, xma, ymi, yma, 0, NULL, NULL, 0, 0, prec, t2, dt, lt, nsimax, conv, cost, gt, kt, x, y, 0);
r=r+erreur;
if (r==*rep)
{
@@ -3688,6 +3688,30 @@ void s_alea_rect(int point_nb, double x[], double y[],
PutRNGstate();
}
+void s_alea_rect_mask(int point_nb, double x[], double y[],
+ double nmask, double *xmask, double *ymask, double dx, double dy,
+ double xmi, double xma, double ymi, double yma,
+ double p)
+ {
+ int i, out;
+ size_t ii;
+ double xr, yr;
+ GetRNGstate();
+ for (i=0; i<point_nb; i++)
+ {
+ ii=nmask*unif_rand();
+ out=1;
+ while (out)
+ {
+ x[i]=*(xmask+ii)+((2*unif_rand()-1)*(dx/p))*p;
+ y[i]=*(ymask+ii)+((2*unif_rand()-1)*(dy/p))*p;
+ out=(x[i]<xmi)||(x[i]>xma)||(y[i]<ymi)||(y[i]>yma);
+ }
+ //Rprintf("s_alea_rect_mask, ii %d, xmask %f, ymask %f, x %f, y %f\n", ii, *(xmask+ii), *(ymask+ii), x[i], y[i]);
+ }
+ PutRNGstate();
+ }
+
/*pour un zone circulaire*/
void s_alea_disq(int point_nb, double *x, double *y, double x0, double y0, double r0, double p)
{
@@ -7741,11 +7765,14 @@ int randomdist(int *vec, int nb_type, double *mat, double *matp)
/*Mimetic point process as in Goreaud et al. 2004 */
/******************************************************/
-int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *xmi, double *xma, double *ymi, double *yma,
+int mimetic_rect(int *point_nb, double *x, double *y, double *surface,
+ double *xmi, double *xma, double *ymi, double *yma,
+ double *nmask, double *xmask, double *ymask, double *dx, double *dy,
double *prec, int *t2, double *dt, double *lobs, int *nsimax, int *conv, double *cost,
double *g, double *k, double *xx, double *yy, int *mess)
{
int i, compteur_c=0, r=0, erreur=0;
+ size_t ii;
int compteur=0;
double *l;
double cout, cout_c;
@@ -7756,9 +7783,30 @@ int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *x
double offsetx=((*xmi)<0)?(*xmi):0.;
double offsety=((*ymi)<0)?(*ymi):0.;
decalRect(*point_nb, x, y, xmi, xma, ymi, yma);
+ if ((*nmask)>0)
+ {
+ if (offsetx<0.)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ double avant= *(xmask+ii);
+ *(xmask+ii)-=offsetx;
+ }
+ }
+ if (offsety<0.)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ *(ymask+ii)-=offsety;
+ }
+ }
+ }
//creation of a initial point pattern and cost
- s_alea_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, *prec);
+ if ((*nmask)>0)
+ s_alea_rect_mask(*point_nb, x, y, *nmask, xmask, ymask, *dx, *dy, *xmi, *xma, *ymi, *yma, *prec);
+ else
+ s_alea_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, *prec);
erreur=ripley_rect(point_nb, x, y, xmi, xma, ymi, yma, t2, dt, g, k);
if (erreur!=0) return -1;
cout=0;
@@ -7774,7 +7822,7 @@ int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *x
cout_c=0;
while (compteur<*nsimax)
{
- cout_c=echange_point_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, intensity, *prec, cout, lobs, t2, dt, g, k);
+ cout_c=echange_point_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, *nmask, xmask, ymask, *dx, *dy, intensity, *prec, cout, lobs, t2, dt, g, k);
if (cout==cout_c)
compteur_c++;
else
@@ -7810,17 +7858,40 @@ int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *x
*xma+=offsetx;
*ymi+=offsety;
*yma+=offsety;
+ // Shift mask into original location
+ if ((*nmask)>0)
+ {
+ if (offsetx<0)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ xmask[ii]=xmask[ii]+offsetx;
+ }
+ }
+ if (offsety<0)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ ymask[ii]=ymask[ii]+offsety;
+ }
+ }
+ }
free(l);
return r;
}
-double echange_point_rect(int point_nb, double *x, double *y, double xmi, double xma, double ymi, double yma, double intensity, double p, double cout, double * lobs, int *t2, double *dt, double *g, double *k)
+double echange_point_rect(int point_nb, double *x, double *y,
+ double xmi, double xma, double ymi, double yma,
+ double nmask, double *xmask, double *ymask, double dx, double dy,
+ double intensity, double p, double cout, double * lobs,
+ int *t2, double *dt, double *g, double *k)
{
double xr, yr, xcent[4], ycent[4], n_cout[4], *l, xprec, yprec;
- int erreur, max, i, j, num;
+ int erreur, max, i, j, num, out;
+ size_t ii;
vecalloc(&l, *t2);
GetRNGstate();
- num=unif_rand()* (point_nb); // numero de l'arbre que l'on retire
+ num=unif_rand()*(point_nb); // numero de l'arbre que l'on retire
xprec=x[num];
yprec=y[num];
xr=xma-xmi;
@@ -7835,8 +7906,22 @@ double echange_point_rect(int point_nb, double *x, double *y, double xmi, double
for (j=0; j<4; j++)
{
- xcent[j]=xmi+(unif_rand()*(xr/p))*p; // coordonn�es (x,y) du point tir�
- ycent[j]=ymi+(unif_rand()*(yr/p))*p;
+ if (nmask>0)
+ {
+ ii=nmask*unif_rand();
+ out=1;
+ while (out)
+ {
+ xcent[j]=xmask[ii]+((2*unif_rand()-1)*(dx/p))*p;
+ ycent[j]=ymask[ii]+((2*unif_rand()-1)*(dy/p))*p;
+ out=(xcent[j]<xmi)||(xcent[j]>xma)||(ycent[j]<ymi)||(ycent[j]>yma);
+ }
+ }
+ else
+ {
+ xcent[j]=xmi+(unif_rand()*(xr/p))*p; // coordonn�es (x,y) du point tir�
+ ycent[j]=ymi+(unif_rand()*(yr/p))*p;
+ }
x[num]=xcent[j];
y[num]=ycent[j];
erreur=ripley_rect(&point_nb, x, y, &xmi, &xma, &ymi, &yma, t2, dt, g, k);
diff --git a/src/Zlibs.h b/src/Zlibs.h
index 9130b52accd90d87d9d6aa0063d047ed99d007f1..7b161558a407df1bc2abe5c4c4d91e89e359eabd 100755
--- a/src/Zlibs.h
+++ b/src/Zlibs.h
@@ -78,6 +78,7 @@ int intertypelocal_tr_disq(int *,double *,double *,int *,double *,double *,doubl
double *,double *,double *,double *,double *,double *,int *,double *,double *,double *);
void s_alea_rect(int,double[],double[],double,double,double,double,double);
+void s_alea_rect_mask(int,double[],double[],double,double *,double *,double,double,double,double,double,double,double);
void s_alea_disq(int,double *,double *,double,double,double,double);
void s_alea_tr_rect(int,double *,double *,double,double,double,double,int,double *,double *,double *,double *,
double *,double *,double);
@@ -155,13 +156,13 @@ int rao_tr_disq(int *,double *,double *,double *,double *,double *,int *,double
int rao_tr_disq_ic(int *,double *,double *,double *,double *,double *,int *,double *,double *,double *,double *,double *,double *,int *,double *,int *,int *,double *,
int *,int *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,int *);
-int mimetic_rect(int *,double *,double *, double *,double *,double *,double *,double *,double *, int *,double *,double *,int *,int *,double *,double *,double *,double *,double *,int *);
+int mimetic_rect(int *,double *,double *, double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *, int *,double *,double *,int *,int *,double *,double *,double *,double *,double *,int *);
int mimetic_disq(int *,double *,double *,double *,double *,double *,double *,double *, int *, double *, double *, int *, int *, double *,double *,double *,double *,double *,int *);
int mimetic_tr_rect(int *,double *,double *, double *,double *,double *,double *,double *,int *, double *, double *, double *, double *, double *, double *,
double *, int *, double *, double *, int *, int *, double *,double *, double *,double *,double *,int *);
int mimetic_tr_disq(int *,double *,double *, double *,double *,double *,double *,int *, double *,double *,double *,double *,double *,double *cy,
double *, int *, double *, double *, int *, int *, double *,double *,double *,double *,double *,int *);
-double echange_point_rect(int,double *,double *,double,double,double,double,double,double,double,double *,int *,double *,double *,double *);
+double echange_point_rect(int,double *,double *,double,double,double,double,double,double *,double *,double,double,double,double,double,double *,int *,double *,double *,double *);
double echange_point_disq(int,double *,double *,double,double,double,double,double,double,double *,int *,double *,double *,double *);
double echange_point_tr_rect(int,double *,double *,double,double,double,double,int *,double *,double *,double *,double *,double *,double *,
double,double,double,double *,int *,double *,double *,double *);