[GRASS-SVN] r33136 - in grass-addons: . vector
vector/v.adehabitat.clusthr vector/v.adehabitat.kernelUD
vector/v.adehabitat.mcp
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 28 03:02:08 EDT 2008
Author: neteler
Date: 2008-08-28 03:02:08 -0400 (Thu, 28 Aug 2008)
New Revision: 33136
Added:
grass-addons/vector/v.adehabitat.clusthr/
grass-addons/vector/v.adehabitat.clusthr/Makefile
grass-addons/vector/v.adehabitat.clusthr/description.html
grass-addons/vector/v.adehabitat.clusthr/main.c
grass-addons/vector/v.adehabitat.kernelUD/
grass-addons/vector/v.adehabitat.kernelUD/Makefile
grass-addons/vector/v.adehabitat.kernelUD/description.html
grass-addons/vector/v.adehabitat.kernelUD/function.c
grass-addons/vector/v.adehabitat.kernelUD/global.h
grass-addons/vector/v.adehabitat.kernelUD/main.c
grass-addons/vector/v.adehabitat.mcp/
grass-addons/vector/v.adehabitat.mcp/Makefile
grass-addons/vector/v.adehabitat.mcp/description.html
grass-addons/vector/v.adehabitat.mcp/main.c
Modified:
grass-addons/contributors.csv
Log:
Cl?\195?\169ment Calenge <calenge biomserv univ-lyon1.fr>: ported some adehabitat R package functions to GRASS
Modified: grass-addons/contributors.csv
===================================================================
--- grass-addons/contributors.csv 2008-08-28 01:06:15 UTC (rev 33135)
+++ grass-addons/contributors.csv 2008-08-28 07:02:08 UTC (rev 33136)
@@ -7,6 +7,7 @@
-,Tomas Cebecauer,<tomas.cebecauer jrc.it>,tomce,yes
-,Dylan Beaudette,<dylan.beaudette gmail.com>,dylan,yes
-,William Richard,<willster3021 gmail.com>,willr,yes
+-,Clément Calenge,<calenge biomserv univ-lyon1.fr>-,-
bundala,Daniel Bundala,<bundala gmail.com>,-,yes
cannata,Massimiliano Cannata,<massimiliano.cannata supsi.ch>,maxi,yes
chemin,Yann Chemin,<ychemin gmail com>,ychemin,yes
Added: grass-addons/vector/v.adehabitat.clusthr/Makefile
===================================================================
--- grass-addons/vector/v.adehabitat.clusthr/Makefile (rev 0)
+++ grass-addons/vector/v.adehabitat.clusthr/Makefile 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.adehabitat.clusthr
+
+LIBES = $(SITESLIB) $(VECTLIB) $(GISLIB) $(CURSES)
+DEPENDENCIES= $(VECTDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/vector/v.adehabitat.clusthr/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/vector/v.adehabitat.clusthr/description.html
===================================================================
--- grass-addons/vector/v.adehabitat.clusthr/description.html (rev 0)
+++ grass-addons/vector/v.adehabitat.clusthr/description.html 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,92 @@
+<h2>DESCRIPTION</h2>
+
+<p><em>v.adehabitat.clusthr</em> computes the home range of one animal using
+a clustering algorithm (Kenward et al. 2001). This method has been
+"translated" for GRASS from the R function "clusthr" in the
+package adehabitat (Calenge, 2006). This method estimates
+home range using a modification of single-linkage cluster analysis
+developped by Kenward et al. (2001). The clustering process is
+described hereafter: the three locations with the minimum mean
+of nearest-neighbour joining distances (NNJD) form the first cluster.
+At each step, two distances are computed: (i) the minimum mean
+NNJD between three locations (which corresponds to the next
+potential cluster) and (ii) the minimum of the NNJD between a
+cluster "c" and the closest location. If (i) is smaller that (ii),
+another cluster is defined with these three locations. If (ii)
+is smaller than (i), the cluster "c" gains a new location.
+If this new location belong to another cluster, the two cluster fuses.
+The process stop when all relocations are assigned to the same cluster.</p>
+
+<p>At each step of the clustering process, the proportion of all
+relocations which are assigned to a cluster is computed (so that
+the home range can be defined to enclose a given proportion of the
+relocations at hand, i.e. to an uncomplete process). At a given
+step, the home range is defined as the set of minimum convex
+polygon enclosing the relocations in the clusters.</p>
+
+<p>The user may either: (i) compute the home range at a given
+percentage: this corresponds to the set of clusters enclosing a
+given percentage of the relocations (default), (ii) compute the
+home range corresponding to a given step of the algorithm (flag -s, a
+first run of the module indicates the maximum number of steps
+supported by the data), or (iii) have a visual display of the
+whole clustering algorithm (flag -a). Note that for the first
+case (i), we may have different home range configuration for a given
+value of 'percent'. For example, if percent = 100, it is
+possible that all relocations have been assigned to one out of
+two clusters, just before the last step of the algorithm. So that the
+last step of the algorithm just consists into the "fusing" of the
+two clusters. The two choices (two clusters or one cluster)
+correspond to configurations where 100% of the relocations
+are enclosed in the home range. When such
+cases happen, the module returns the configuration with the
+largest number of clusters, i.e., the minimum home range size
+(in the example, two clusters of relocations).</p>
+
+
+
+<H2>REFERENCES</H2>
+
+<p><EM>Kenward R.E., Clarke R.T., Hodder K.H. and Walls S.S. (2001)
+Density and linkage estimators of homre range: nearest neighbor
+clustering defines multinuclear cores. Ecology, 82, 1905-1920.</EM></p>
+
+<p><EM>Calenge, C. (2006) The package adehabitat for the R software: a tool
+for the analysis of space and habitat use by animals. Ecological
+Modelling, 197, 516-519.</EM></p>
+
+
+<h2>EXAMPLE</h2>
+
+Estimation of the 95% home range:<br>
+
+<div class="code"><pre>
+v.adehabitat.clusthr input=localisations output=homerange percent=95
+</pre></div>
+
+Estimation of the home range corresponding to the 25th step of
+the clustering algorithm:<br>
+
+<div class="code"><pre>
+v.adehabitat.clusthr input=localisations output=homerange step=25 -s
+</pre></div>
+
+Summary of the clustering algorithm (show all the steps):<br>
+
+<div class="code"><pre>
+v.adehabitat.clusthr input=localisations output=homerange -a
+</pre></div>
+
+
+<h2>SEE ALSO</h2>
+
+<em><a HREF="v.hull.html">v.hull</a></em>,
+<em><a HREF="v.adehabitat.mcp.html">v.adehabitat.mcp</a></em>,
+<em><a HREF="v.adehabitat.kernelUD.html">v.adehabitat.kernelUD</a></em>
+
+
+<h2>AUTHOR</h2>
+
+Clement Calenge
+
+<p><i>Last changed: $Date: 2007/08/23 17:12:30 $</i>
Added: grass-addons/vector/v.adehabitat.clusthr/main.c
===================================================================
--- grass-addons/vector/v.adehabitat.clusthr/main.c (rev 0)
+++ grass-addons/vector/v.adehabitat.clusthr/main.c 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,1198 @@
+/******************************************************************************
+*
+* Computes the Cluster home range from a file of relocations
+*
+* 5 September 2007
+*
+* Clement Calenge
+* This file is part of GRASS GIS. It 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.
+******************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/site.h>
+#include <grass/glocale.h>
+
+
+/* *************************************************************
+
+* The sources from ADE-4
+
+**************************************************************** */
+
+
+void vecalloc (double **vec, int n)
+/*--------------------------------------------------
+ * Memory Allocation for a vector of length n
+ --------------------------------------------------*/
+{
+ if ( (*vec = (double *) calloc(n+1, sizeof(double))) != 0) {
+ **vec = n;
+ return;
+ } else {
+ return;
+ }
+}
+
+/*****************/
+void vecintalloc (int **vec, int n)
+/*--------------------------------------------------
+ * Memory allocation for an integer vector of length n
+ --------------------------------------------------*/
+{
+ if ( (*vec = (int *) calloc(n+1, sizeof(int))) != NULL) {
+ **vec = n;
+ return;
+ } else {
+ return;
+ }
+}
+
+
+void taballoc (double ***tab, int l1, int c1)
+/*--------------------------------------------------
+ * Dynamic Memory Allocation for a table (l1, c1)
+ --------------------------------------------------*/
+{
+ int i, j;
+
+ if ( (*tab = (double **) calloc(l1+1, sizeof(double *))) != 0) {
+ for (i=0;i<=l1;i++) {
+ if ( (*(*tab+i)=(double *) calloc(c1+1, sizeof(double))) == 0 ) {
+ return;
+ for (j=0;j<i;j++) {
+ free(*(*tab+j));
+ }
+ }
+ }
+ }
+
+ **(*tab) = l1;
+ **(*tab+1) = c1;
+}
+
+
+
+
+
+/***********************************************************************/
+void freevec (double *vec)
+/*--------------------------------------------------
+ * Free memory for a vector
+ --------------------------------------------------*/
+{
+ free((char *) vec);
+}
+
+/***********************************************************************/
+void freeintvec (int *vec)
+/*--------------------------------------------------
+* Free memory for an integer vector
+--------------------------------------------------*/
+{
+
+ free((char *) vec);
+
+}
+
+
+
+
+void freetab (double **tab)
+/*--------------------------------------------------
+ * Free memory for a table
+ --------------------------------------------------*/
+{
+ int i, n;
+
+ n = *(*(tab));
+ for (i=0;i<=n;i++) {
+ free((char *) *(tab+i) );
+ }
+ free((char *) tab);
+}
+
+
+
+/* *************************************************************
+
+* Code from v.hull
+
+**************************************************************** */
+
+
+int read_points( struct Map_info *In, double ***coordinate)
+{
+ int line, nlines, npoints, ltype, i = 0;
+ double **xySites;
+ static struct line_pnts *Points = NULL;
+
+ if (!Points)
+ Points = Vect_new_line_struct ();
+
+ /* Allocate array of pointers */
+ npoints = Vect_get_num_primitives(In,GV_POINT);
+ xySites = (double **) G_calloc ( npoints, sizeof(double*) );
+
+ nlines = Vect_get_num_lines(In);
+
+ for ( line = 1; line <= nlines; line++){
+ ltype = Vect_read_line (In, Points, NULL, line);
+ if ( !(ltype & GV_POINT ) ) continue;
+
+ xySites[i] = (double *)G_calloc((size_t)2, sizeof(double));
+
+ xySites[i][0] = Points->x[0];
+ xySites[i][1] = Points->y[0];
+ i++;
+ }
+
+ *coordinate = xySites;
+
+ return (npoints);
+}
+
+
+
+
+
+struct Point {
+ double x;
+ double y;
+};
+
+
+int rightTurn(struct Point *P, int i, int j, int k) {
+ double a, b, c, d;
+ a = P[i].x - P[j].x;
+ b = P[i].y - P[j].y;
+ c = P[k].x - P[j].x;
+ d = P[k].y - P[j].y;
+ return a*d - b*c < 0;
+}
+
+
+int cmpPoints(const void* v1, const void* v2) {
+ struct Point *p1, *p2;
+ p1 = (struct Point*) v1;
+ p2 = (struct Point*) v2;
+ if( p1->x > p2->x )
+ return 1;
+ else if( p1->x < p2->x )
+ return -1;
+ else
+ return 0;
+}
+
+
+int convexHull(struct Point* P, const int numPoints, int **hull) {
+ int pointIdx, upPoints, loPoints;
+ int *upHull, *loHull;
+
+ /* sort points in ascending x order*/
+ qsort(P, numPoints, sizeof(struct Point), cmpPoints);
+
+ *hull = (int*) G_malloc(numPoints * 2 * sizeof(int));
+
+ /* compute upper hull */
+ upHull = *hull;
+ upHull[0] = 0;
+ upHull[1] = 1;
+ upPoints = 1;
+ for(pointIdx = 2; pointIdx < numPoints; pointIdx++) {
+ upPoints++;
+ upHull[upPoints] = pointIdx;
+ while( upPoints > 1 &&
+ !rightTurn(P, upHull[upPoints], upHull[upPoints-1],
+ upHull[upPoints-2])
+ ) {
+ upHull[upPoints-1] = upHull[upPoints];
+ upPoints--;
+ }
+ }
+
+ /*
+ printf("upPoints: %d\n", upPoints);
+ for(pointIdx = 0; pointIdx <= upPoints; pointIdx ++)
+ printf("%d ", upHull[pointIdx]);
+ printf("\n");
+ */
+
+ /* compute lower hull, overwrite last point of upper hull */
+ loHull = &(upHull[upPoints]);
+ loHull[0] = numPoints - 1;
+ loHull[1] = numPoints - 2;
+ loPoints = 1;
+ for(pointIdx = numPoints - 3; pointIdx >= 0; pointIdx--) {
+ loPoints++;
+ loHull[loPoints] = pointIdx;
+ while( loPoints > 1 &&
+ !rightTurn(P, loHull[loPoints], loHull[loPoints-1],
+ loHull[loPoints-2])
+ ) {
+ loHull[loPoints-1] = loHull[loPoints];
+ loPoints--;
+ }
+ }
+
+ G_debug(3, "numPoints:%d loPoints:%d upPoints:%d",
+ numPoints, loPoints, upPoints);
+ /*
+ printf("loPoints: %d\n", loPoints);
+ for(pointIdx = 0; pointIdx <= loPoints; pointIdx ++)
+ printf("%d ", loHull[pointIdx]);
+ printf("\n");
+ */
+
+ /* reclaim uneeded memory */
+ *hull = (int *) G_realloc(*hull, (loPoints + upPoints) * sizeof(int));
+ return loPoints + upPoints;
+}
+
+
+
+
+
+/* *********************************************************************
+ * *
+ * Home range by Clustering (Kenward et al. 2001) *
+ * Sources from adehabitat *
+ * *
+ ***********************************************************************/
+
+
+/* finds the cluster with the minimum average distance between the 3 points
+ not assigned to a cluster */
+
+void trouveclustmin(double **xy, int *clust, int *lo1, int *lo2,
+ int *lo3, double *dist)
+{
+ /* Declaration */
+ int i, j, k, m, npas, nr, *indice;
+ double **xy2, di1, di2, di3, ditmp;
+
+ /* Memory allocation */
+ nr = (int) xy[0][0];
+ npas = 0;
+ di1 = 0;
+ di2 = 0;
+ di3 = 0;
+ ditmp = 0;
+
+ /* Number of non assigned points */
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] == 0) {
+ npas++;
+ }
+ }
+ taballoc(&xy2, npas, 2);
+ vecintalloc(&indice, npas);
+
+ /* The non assigned points are stored in xy2 */
+ k = 1;
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] == 0) {
+ xy2[k][1] = xy[i][1];
+ xy2[k][2] = xy[i][2];
+ indice[k] = i;
+ k++;
+ }
+ }
+
+ /* Computes the distane between the relocations */
+ *dist = 0;
+ m=0;
+ for (i = 1; i <= (npas-2); i++) {
+ for (j = (i+1); j <= (npas-1); j++) {
+ for (k = (j+1); k <= npas; k++) {
+ di1 = sqrt((xy2[i][1] - xy2[j][1]) * (xy2[i][1] - xy2[j][1]) +
+ (xy2[i][2] - xy2[j][2]) * (xy2[i][2] - xy2[j][2]) );
+ di2 = sqrt((xy2[i][1] - xy2[k][1]) * (xy2[i][1] - xy2[k][1]) +
+ (xy2[i][2] - xy2[k][2]) * (xy2[i][2] - xy2[k][2]));
+ di3 = sqrt((xy2[k][1] - xy2[j][1]) * (xy2[k][1] - xy2[j][1]) +
+ (xy2[k][2] - xy2[j][2]) * (xy2[k][2] - xy2[j][2]));
+ /* average distance */
+ ditmp = (di1 + di2 + di3) / 3;
+ /* minimum distance */
+ if ((m == 0) || (ditmp < *dist)) {
+ *dist = ditmp;
+ *lo1 = indice[i];
+ *lo2 = indice[j];
+ *lo3 = indice[k];
+ }
+ m = 1;
+ }
+ }
+ }
+ /* free memory */
+ freeintvec(indice);
+ freetab(xy2);
+}
+
+
+
+/* Finds the distance between a cluster of points and the nearest point*/
+void nndistclust(double **xy, double *xyp, double *dist)
+{
+ /* Declaration */
+ int n, i, m;
+ double di;
+
+ m = 0;
+ di =0;
+ n = (int) xy[0][0];
+ *dist = 0;
+
+ /* finds the distance and the corresponding point */
+ for (i = 1; i <= n; i++) {
+ di = sqrt( (xy[i][1] - xyp[1]) * (xy[i][1] - xyp[1]) +
+ (xy[i][2] - xyp[2]) * (xy[i][2] - xyp[2]) );
+ if ( (di < *dist) || (m == 0) ) {
+ *dist = di;
+ }
+ m = 1;
+ }
+}
+
+
+/* The function nndistclust is applied for all available clusters */
+void parclust(double **xy, int *clust, int *noclust,
+ int *noloc, double *dist)
+{
+ /* Declaration */
+ int i, k, m, nr2, nr, nocl;
+ double **xy2, *xyp, di, di2;
+
+ /* Memory allocation */
+ nocl = *noclust;
+ nr = xy[0][0];
+ nr2 = 0;
+
+ /* The number of available clusters */
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] == nocl) {
+ nr2++;
+ }
+ }
+
+ taballoc(&xy2, nr2, 2);
+ vecalloc(&xyp, 2);
+
+ /* stores the non assigned points in xy2 */
+ k = 1;
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] == nocl) {
+ xy2[k][1] = xy[i][1];
+ xy2[k][2] = xy[i][2];
+ k++;
+ }
+ }
+
+ /* Finds the minimum distance between a point and a cluster,
+ performed for all clusters */
+ di = 0;
+ di2 = 0;
+ m = 0;
+ *dist = 0;
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] != nocl) {
+ xyp[1] = xy[i][1];
+ xyp[2] = xy[i][2];
+ nndistclust(xy2, xyp, &di);
+ if ( (di < *dist) || (m == 0) ) {
+ *dist = di;
+ *noloc = i;
+ }
+ m = 1;
+ }
+ }
+
+ /* Free memory */
+ freetab(xy2);
+ freevec(xyp);
+}
+
+
+/* The function trouveminclust identifies the cluster for which the nearest
+ point is the closest */
+void trouveminclust(double **xy, int *liclust, int *clust,
+ int *noclust, int *noloc, double *dist)
+{
+ /* Declaration */
+ int i, nr, nc, m, labclust, nolo;
+ double di;
+
+ nr = (int) xy[0][0];
+ nc = 0;
+ di = 0;
+ labclust = 0;
+ nolo = 0;
+
+ /* Assigned clusters */
+ for (i = 1; i <= nr; i++) {
+ if (liclust[i] > 0) {
+ nc++;
+ }
+ }
+
+ /* finds the minimum distance between a cluster and its nearest point
+ (the cluster name and the point ID are searched) */
+ m = 0;
+ *dist = 0;
+ for (i = 1; i <= nc; i++) {
+ labclust = liclust[i];
+ parclust(xy, clust, &labclust, &nolo, &di);
+ if ( (m == 0) || (di < *dist) ) {
+ *dist = di;
+ *noloc = nolo;
+ *noclust = labclust;
+ }
+ m = 1;
+ }
+}
+
+
+/* What should be done: create a new cluster or add a relocation
+ to an existing one ? */
+
+void choisnvclust(double **xy, int *liclust, int *clust, int *ordre)
+{
+ /* Declaration */
+ int i, k, nr, noloat, cluat, nolo1, nolo2, nolo3, maxclust;
+ int maxindiceclust, clu1, *liclub, nz;
+ double dmoyclust, dminloc;
+
+ /* Memory allocation */
+ nz = 0;
+ i = 0;
+ k = 0;
+ nr = (int) xy[0][0];
+ maxclust = 0;
+ maxindiceclust = 0;
+ nolo1 = 0;
+ nolo2 = 0;
+ nolo3 = 0;
+ noloat = 0;
+ cluat = 0;
+ clu1 = 0;
+ vecintalloc(&liclub, nr);
+
+ /* finds the max label for the cluster */
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] != 0) {
+ if (clust[i] > maxclust) {
+ maxclust = clust[i];
+ }
+ if (liclust[i] != 0) {
+ maxindiceclust = i;
+ }
+ }
+ }
+
+ /* Finds the min distance between 3 relocations */
+ trouveminclust(xy, liclust, clust, &cluat, &noloat, &dminloc);
+
+ /* Computes the average distance between the locs of the smaller cluster */
+ /* First, one verifies that there is at least Three non assigned locs */
+ dmoyclust = dminloc +1;
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] == 0) {
+ nz++;
+ }
+ }
+ if (nz > 3) {
+ dmoyclust = 0;
+ trouveclustmin(xy, clust, &nolo1, &nolo2, &nolo3, &dmoyclust);
+ }
+
+ /* First case: A new cluster independent from the others */
+ if (dmoyclust < dminloc) {
+ ordre[nolo1] = 1;
+ ordre[nolo2] = 1;
+ ordre[nolo3] = 1;
+
+ clust[nolo1] = maxclust + 1;
+ clust[nolo2] = maxclust + 1;
+ clust[nolo3] = maxclust + 1;
+
+ liclust[maxindiceclust+1] = maxclust + 1;
+
+ } else {
+ /* Second case: one loc is added to a cluster */
+
+ /* Case 2.1: the loc does not belong to one cluster */
+ if (clust[noloat] == 0) {
+ ordre[noloat] = 1;
+ clust[noloat] = cluat;
+ } else {
+
+ /* Case 2.2: the loc belong to one cluster: fusion */
+ clu1 = clust[noloat];
+ for (i = 1; i <= nr; i++) {
+ if (clust[i] == clu1) {
+ clust[i] = cluat;
+ ordre[i] = 1;
+ }
+ if (liclust[i] == clu1) {
+ liclust[i] = 0;
+ }
+ }
+ /* and cleaning of liclust */
+ k = 1;
+ for (i = 1; i <= nr; i++) {
+ if (liclust[i] != 0) {
+ liclub[k] = liclust[i];
+ k++;
+ }
+ }
+ for (i = 1; i <= nr; i++) {
+ liclust[i] = liclub[i];
+ }
+ }
+ }
+ freeintvec(liclub);
+}
+
+
+/* The main function for home range computation */
+
+void clusterhr(double **xy, int *facso, int *nolocso, int *cluso)
+{
+ /* Declaration */
+ int i, nr, lo1, lo2, lo3, *clust, len, con, *ordre, *liclust, courant;
+ double di;
+
+ /* Memory allocation */
+ courant = 1;
+ nr = (int) xy[0][0];
+ vecintalloc(&clust, nr);
+ vecintalloc(&ordre, nr);
+ vecintalloc(&liclust, nr);
+ lo1 = 0;
+ lo2 = 0;
+ lo3 = 0;
+ di = 0;
+ con = 1;
+ len = 0;
+
+ /* Begin: Search for the first cluster */
+ trouveclustmin(xy, clust, &lo1, &lo2,
+ &lo3, &di);
+
+ clust[lo1] = 1;
+ clust[lo2] = 1;
+ clust[lo3] = 1;
+ liclust[1] = 1;
+ len = 3;
+
+ /* We store it in the output */
+ cluso[1] = 1;
+ cluso[2] = 1;
+ cluso[3] = 1;
+ nolocso[1] = lo1;
+ nolocso[2] = lo2;
+ nolocso[3] = lo3;
+ facso[1] = 1;
+ facso[2] = 1;
+ facso[3] = 1;
+
+ /* Then repeat until all relocations belong to the same cluster */
+ while (con == 1) {
+ courant++;
+
+ for (i = 1; i <= nr; i++) {
+ ordre[i] = 0;
+ }
+
+ choisnvclust(xy, liclust, clust, ordre);
+
+ for (i = 1; i <= nr; i++) {
+ if (ordre[i] != 0) {
+ len++;
+ cluso[len] = clust[i];
+ nolocso[len] = i;
+ facso[len] = courant;
+ }
+ }
+
+ con = 0;
+ for (i = 2; i <= nr; i++) {
+ if (clust[i] != clust[1])
+ con = 1;
+ }
+ if (con == 0) {
+ con = 0;
+ }
+ }
+
+ /* Free memory */
+ freeintvec(clust);
+ freeintvec(ordre);
+ freeintvec(liclust);
+}
+
+
+
+/* Finds the length of the output for the table containing the home range */
+
+void longfacclust(double **xy, int *len2)
+{
+ /* Declaration */
+ int i, nr, lo1, lo2, lo3, *clust, len, con, *ordre, *liclust, courant;
+ double di;
+
+ /* Memory allocation */
+ courant = 1;
+ nr = (int) xy[0][0];
+ vecintalloc(&clust, nr);
+ vecintalloc(&ordre, nr);
+ vecintalloc(&liclust, nr);
+ lo1 = 0;
+ lo2 = 0;
+ lo3 = 0;
+ di = 0;
+ con = 1;
+ len = 0;
+
+ /* Begin: search for the first cluster */
+ trouveclustmin(xy, clust, &lo1, &lo2,
+ &lo3, &di);
+ clust[lo1] = 1;
+ clust[lo2] = 1;
+ clust[lo3] = 1;
+ liclust[1] = 1;
+ len = 3;
+
+ /* Counts the number of rows needed for the table, which will contain the results */
+ while (con == 1) {
+ courant++;
+
+ for (i = 1; i <= nr; i++) {
+ ordre[i] = 0;
+ }
+
+ choisnvclust(xy, liclust, clust, ordre);
+
+ for (i = 1; i <= nr; i++) {
+ if (ordre[i] != 0) {
+ len++;
+ }
+ }
+ con = 0;
+ for (i = 2; i <= nr; i++) {
+ if (clust[i] != clust[1])
+ con = 1;
+ }
+ if (con == 0) {
+ con = 0;
+ }
+ }
+
+ *len2 = len;
+
+ /* Free memory */
+ freeintvec(clust);
+ freeintvec(ordre);
+ freeintvec(liclust);
+}
+
+
+
+
+
+
+
+
+/* *********************************************************************
+ * *
+ * Main code *
+ * *
+ ***********************************************************************/
+
+
+int main(int argc, char **argv) {
+
+ struct GModule *module;
+ struct Option *input, *output, *percent, *step;
+ struct Flag *all;
+ struct Cell_head window;
+ struct Map_info Map, In;
+ struct Point *points;
+ struct line_pnts *Pointsb;
+ struct line_cats *Cats;
+ struct Flag *flag_s, *flag_a;
+
+ char *mapset;
+
+ double **coordinate, **coor, *tmpx, *tmpy;
+ int numSitePoints, numHullPoints;
+ int ntab, *facso, *nolocso, *cluso, *assigne, *numclust;
+ int *nomclust, *nbptsclu, *hull;
+ double perc, tmp, xc, yc;
+ int stepn, nolig, maxstep, i, k, j, m, ok, kk, pointIdx, maxid;
+
+
+ G_gisinit (argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("vector, geometry");
+ module->description = "Uses a GRASS vector points map to compute the cluster home range";
+
+ input = G_define_option ();
+ input->key = "input";
+ input->type = TYPE_STRING;
+ input->required = YES;
+ input->description = "Name of map containing the relocations to be input";
+ input->gisprompt = "old,vector,vector,input";
+
+ output = G_define_option ();
+ output->key = "output";
+ output->type = TYPE_STRING;
+ output->required = YES;
+ output->description = "Name of a vector area map to be output";
+ output->gisprompt = "new,dig,binary file,output";
+
+ percent = G_define_option ();
+ percent->key = "percent";
+ percent->type = TYPE_STRING;
+ percent->required = NO;
+ percent->description = "Percentage of the relocations to keep";
+ percent->answer = "100";
+
+ step = G_define_option ();
+ step->key = "step";
+ step->type = TYPE_STRING;
+ step->required = NO;
+ step->description = "Step number";
+ step->answer = "1";
+
+ flag_s = G_define_flag();
+ flag_s->key = 's';
+ flag_s->description = _("Use step number");
+
+ flag_a = G_define_flag();
+ flag_a->key = 'a';
+ flag_a->description = _("Return all steps");
+
+ /* Arguments OK ? */
+ if (G_parser (argc, argv)) exit (1);
+
+ /* Output name OK? */
+ Vect_check_input_output_name ( input->answer, output->answer,
+ GV_FATAL_EXIT );
+
+ /* récupère la région */
+ G_get_window (&window);
+
+
+ /* open relocations file */
+ if ((mapset = G_find_vector2 (input->answer, "")) == NULL)
+ G_fatal_error (_("Could not find input map '%s'."), input->answer);
+
+ Vect_set_open_level (2); /* Essential before opening a map!
+ defines the level */
+ Vect_open_old (&In, input->answer, mapset); /* Opens the vector map */
+
+
+ /* checks for legal file name */
+ if(G_legal_filename( output->answer ) < 0)
+ G_fatal_error(_("illegal file name [%s]."), output->answer);
+
+
+ /* load coordinates */
+ numSitePoints = read_points ( &In, &coordinate );
+
+ /* verifier l'allocation de memoire pour coordinate
+ mais ça a l'air d'être bon */
+
+ if(numSitePoints < 3 )
+ G_fatal_error ("Cluster Home Range calculation requires at least three points");
+
+ /* load the percent or the step */
+ if (!(flag_s->answer)) {
+ perc = atof(percent->answer);
+ } else {
+ stepn = ((int) atof(step->answer));
+ }
+
+
+ /* create vector file */
+ if (0 > Vect_open_new (&Map, output->answer, 0) )
+ G_fatal_error ("Unable to open vector file <%s>\n", output->answer);
+ Vect_hist_command ( &Map );
+
+
+ /* ***************************************************
+
+ * Ici, on passe aux instructions pour ADE-4 *
+
+ ****************************************************** */
+
+ /* Les coordonees */
+ taballoc(&coor, numSitePoints, 2);
+
+ for (i = 1; i <= numSitePoints; i++) {
+ coor[i][1] = coordinate[i-1][0];
+ coor[i][2] = coordinate[i-1][1];
+ }
+
+
+ /* size of the output */
+ longfacclust(coor, &ntab);
+
+ /* Memory allocation */
+ vecintalloc(&facso, ntab);
+ vecintalloc(&nolocso, ntab);
+ vecintalloc(&cluso, ntab);
+
+ /* Clustering */
+ clusterhr(coor, facso, nolocso, cluso);
+
+
+ /* Le nombre de lignes à facso */
+ nolig = facso[0];
+
+ /* Le nombre d'étapes à l'algorithme */
+ maxstep = facso[nolig];
+ G_message(_("Maximum step = %d"), maxstep);
+ if (flag_s->answer) {
+ if (stepn > maxstep) {
+ G_fatal_error(_("The step number should be lower than the max = %d."), maxstep);
+ }
+ if (stepn < 1) {
+ G_fatal_error(_("The step number should be greater than 1"));
+ }
+
+ }
+
+ /* On détermine la limite du numero de step correspondant
+ au pourcent donne et le numero de ligne max correspondant */
+
+ if (!(flag_a->answer)) {
+ tmp=0.;
+ vecintalloc(&assigne, numSitePoints);
+ vecintalloc(&numclust, numSitePoints);
+
+ if (!(flag_s->answer)) {
+ for (i = 1; i <= maxstep; i++) {
+ if (tmp < perc) {
+ for (j = 1; j <= nolig; j++) {
+ if (facso[j] == i) {
+ if (assigne[nolocso[j]] == 0) {
+ tmp += 100. * (1. / ((double) numSitePoints));
+ assigne[nolocso[j]] = 1;
+ }
+ /* dans tous les cas */
+ maxid = j;
+ }
+ }
+ }
+ }
+ } else {
+ for (j = 1; j <= nolig; j++) {
+ if (facso[j] == stepn) {
+ maxid = j;
+ }
+ }
+
+ }
+
+ /* on pose assigne = 0 pour tout j */
+ for (j = 1; j <= numSitePoints; j++) {
+ assigne[j] = 0;
+ }
+
+ /* Puis, on determine les clusters encore pas "resolved" pour
+ la limite donnée */
+ vecintalloc(&numclust, numSitePoints);
+ for (j = 1; j <= numSitePoints; j++) {
+ numclust[j] = 0;
+ }
+
+ for (i = maxid; i>=1; i--) {
+ j = nolocso[i];
+ if (assigne[j] == 0) {
+ assigne[j] = 1;
+ numclust[j] = cluso[i];
+ }
+ }
+
+
+
+ /* et leur nombre */
+ k = 1;
+ vecintalloc(&nomclust, numSitePoints);
+ for (i = 1; i <= numSitePoints; i++) {
+ ok = 0;
+ for (j = 1; j <= k; j++) {
+ if (nomclust[j] == numclust[i]) {
+ if (nomclust[j] != 0) {
+ ok = 1;
+ }
+ }
+ }
+ if (ok == 0) {
+ if (numclust[i] != 0) {
+ nomclust[k] = numclust[i];
+ k++;
+ }
+ }
+ }
+ k--;
+
+
+ /* et le nombre de locs pour chacun */
+ vecintalloc(&nbptsclu, k);
+ for (i = 1; i <= k; i++) {
+ for (j = 1; j <= numSitePoints; j++) {
+ if (numclust[j] == nomclust[i])
+ nbptsclu[i]++;
+ }
+ }
+
+
+ /* *****************************************************
+ * *
+ * Enfin, on repasse à GRASS: on calcule un MCP *
+ * par cluster de localisations, que le sort dans *
+ * la carte de sortie *
+ * *
+ * ***************************************************** */
+
+ /* Pour chaque cluster */
+
+ for (i = 1; i <= k; i++) {
+
+ /* memory allocation for the points */
+ points = (struct Point *) G_malloc((nbptsclu[i]) * sizeof(struct Point));
+ /* gets the relocations corresponding to the current cluster */
+ kk=0;
+ for (j = 1; j <= numSitePoints; j++) {
+ if (numclust[j] == nomclust[i]) {
+ points[kk].x = coor[j][1];
+ points[kk].y = coor[j][2];
+ kk++;
+ }
+ }
+
+ /* computes the convex hull */
+ numHullPoints = convexHull(points, nbptsclu[i], &hull);
+
+ /* Memory allocation for tmpx and tmpy */
+ tmpx = (double *) G_malloc((numHullPoints+1) * sizeof(double));
+ tmpy = (double *) G_malloc((numHullPoints+1) * sizeof(double));
+
+ /* gets the coordinates of the vertices of the hull */
+ xc = 0;
+ yc = 0;
+ for (j = 0; j < numHullPoints; j++) {
+ pointIdx = hull[j];
+ tmpx[j] = points[pointIdx].x;
+ tmpy[j] = points[pointIdx].y;
+ xc += tmpx[j] / ((double) numHullPoints);
+ yc += tmpy[j] / ((double) numHullPoints);
+ }
+
+ /* adds the first point of the hull at the end, to have
+ a closed boundary */
+ tmpx[numHullPoints] = points[hull[0]].x;
+ tmpy[numHullPoints] = points[hull[0]].y;
+
+ /* Creates a new structure for lines and categories */
+ Pointsb = Vect_new_line_struct ();
+ Cats = Vect_new_cats_struct ();
+
+ /* Now copies the vertices in the line structure */
+ Vect_copy_xyz_to_pnts(Pointsb, tmpx, tmpy, 0, numHullPoints+1);
+
+ /* And write the line in the map */
+ Vect_write_line (&Map, GV_BOUNDARY, Pointsb, Cats);
+
+ /* */
+ Vect_reset_line (Pointsb);
+ Vect_append_point (Pointsb,xc, yc, 0.0);
+ Vect_cat_set ( Cats, 1, i );
+ Vect_write_line (&Map, GV_CENTROID, Pointsb, Cats);
+
+ /* remove everything that has been created in this loop */
+
+ Vect_destroy_line_struct(Pointsb);
+ Vect_destroy_cats_struct(Cats);
+ G_free (points);
+ G_free (tmpx);
+ G_free (tmpy);
+ G_free (hull);
+ }
+ freeintvec(nbptsclu);
+
+
+
+ } else {
+
+ vecintalloc(&assigne, numSitePoints);
+ vecintalloc(&numclust, numSitePoints);
+ vecintalloc(&numclust, numSitePoints);
+ vecintalloc(&nomclust, numSitePoints);
+
+ for (m = 1; m <= maxstep; m++) {
+ tmp=0.;
+
+ for (j = 1; j <= nolig; j++) {
+ if (facso[j] == m) {
+ maxid = j;
+ }
+ }
+
+
+ /* on pose assigne = 0 pour tout j */
+ for (j = 1; j <= numSitePoints; j++) {
+ assigne[j] = 0;
+ }
+
+ /* Puis, on determine les clusters encore pas "resolved" pour
+ la limite donnée */
+ for (j = 1; j <= numSitePoints; j++) {
+ numclust[j] = 0;
+ }
+ for (j = 1; j <= numSitePoints; j++) {
+ nomclust[j] = 0;
+ }
+
+ for (i = maxid; i>=1; i--) {
+ j = nolocso[i];
+ if (assigne[j] == 0) {
+ assigne[j] = 1;
+ numclust[j] = cluso[i];
+ }
+ }
+
+
+ /* et leur nombre */
+ k = 1;
+ for (i = 1; i <= numSitePoints; i++) {
+ ok = 0;
+ for (j = 1; j <= k; j++) {
+ if (nomclust[j] == numclust[i]) {
+ if (nomclust[j] != 0) {
+ ok = 1;
+ }
+ }
+ }
+ if (ok == 0) {
+ if (numclust[i] != 0) {
+ nomclust[k] = numclust[i];
+ k++;
+ }
+ }
+ }
+ k--;
+
+ if (m == 64) {
+ for (i = 1; i <= numSitePoints; i++)
+ G_message(_("Maximum step = %d"), nomclust[i]);
+ }
+
+
+ /* et le nombre de locs pour chacun */
+ vecintalloc(&nbptsclu, k);
+ for (i = 1; i <= k; i++) {
+ for (j = 1; j <= numSitePoints; j++) {
+ if (numclust[j] == nomclust[i])
+ nbptsclu[i]++;
+ }
+ }
+
+
+ /* *****************************************************
+ * *
+ * Enfin, on repasse à GRASS: on calcule un MCP *
+ * par cluster de localisations, que le sort dans *
+ * la carte de sortie *
+ * *
+ * ***************************************************** */
+
+ /* Pour chaque cluster */
+
+ for (i = 1; i <= k; i++) {
+
+ /* memory allocation for the points */
+ points = (struct Point *) G_malloc((nbptsclu[i]) * sizeof(struct Point));
+ /* gets the relocations corresponding to the current cluster */
+ kk=0;
+ for (j = 1; j <= numSitePoints; j++) {
+ if (numclust[j] == nomclust[i]) {
+ points[kk].x = coor[j][1];
+ points[kk].y = coor[j][2];
+ kk++;
+ }
+ }
+
+ /* computes the convex hull */
+ numHullPoints = convexHull(points, nbptsclu[i], &hull);
+
+ /* Memory allocation for tmpx and tmpy */
+ tmpx = (double *) G_malloc((numHullPoints+1) * sizeof(double));
+ tmpy = (double *) G_malloc((numHullPoints+1) * sizeof(double));
+
+ /* gets the coordinates of the vertices of the hull */
+ for (j = 0; j < numHullPoints; j++) {
+ pointIdx = hull[j];
+ tmpx[j] = points[pointIdx].x;
+ tmpy[j] = points[pointIdx].y;
+ }
+
+ /* adds the first point of the hull at the end, to have
+ a closed line */
+ tmpx[numHullPoints] = points[hull[0]].x;
+ tmpy[numHullPoints] = points[hull[0]].y;
+
+ /* Creates a new structure for lines and categories */
+ Pointsb = Vect_new_line_struct ();
+ Cats = Vect_new_cats_struct ();
+
+ /* Now copies the vertices in the line structure */
+ Vect_copy_xyz_to_pnts(Pointsb, tmpx, tmpy, 0, numHullPoints+1);
+
+ /* And write the line in the map */
+ Vect_cat_set ( Cats, 1, m );
+ Vect_write_line (&Map, GV_LINE, Pointsb, Cats);
+
+ /* remove everything that has been created in this loop */
+ Vect_destroy_line_struct(Pointsb);
+ Vect_destroy_cats_struct(Cats);
+ G_free (points);
+ G_free (tmpx);
+ G_free (tmpy);
+ G_free (hull);
+ }
+ }
+ freeintvec(nbptsclu);
+ }
+
+
+ /* clean up and bye bye */
+ Vect_build (&Map, stdout);
+ Vect_close (&Map);
+ freetab(coor);
+ freeintvec(facso);
+ freeintvec(nolocso);
+ freeintvec(cluso);
+ freeintvec(assigne);
+ freeintvec(numclust);
+ freeintvec(nomclust);
+
+}
+
+
Property changes on: grass-addons/vector/v.adehabitat.clusthr/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/vector/v.adehabitat.kernelUD/Makefile
===================================================================
--- grass-addons/vector/v.adehabitat.kernelUD/Makefile (rev 0)
+++ grass-addons/vector/v.adehabitat.kernelUD/Makefile 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,15 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../../
+
+PGM = v.adehabitat.kernelUD
+
+LIBES = $(VECTLIB) $(GISLIB)
+DEPENDENCIES= $(VECTDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/vector/v.adehabitat.kernelUD/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/vector/v.adehabitat.kernelUD/description.html
===================================================================
--- grass-addons/vector/v.adehabitat.kernelUD/description.html (rev 0)
+++ grass-addons/vector/v.adehabitat.kernelUD/description.html 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,100 @@
+<H2>DESCRIPTION</H2>
+
+<p><em>v.adehabitat.kernelUD</em> is used to estimate the utilization distribution (UD)
+of animals monitored by radio-tracking, with the classical kernel
+method.</p>
+
+<p>The Utilization Distribution (UD) is the bivariate function giving
+the probability density that an animal is found at a point
+according to its geographical coordinates. Using this model, one
+can define the home range as the minimum area in which an animal
+has some specified probability of being located. The module v.adehabitat.kernelUD
+correspond to the approach described in Worton (1995).</p>
+
+<p>The kernel method has been recommended by many authors for the
+estimation of the utilization distribution (e.g. Worton, 1989,
+1995). The default method for the estimation of the smoothing
+parameter is the <em>ad hoc</em> method, i.e. for a bivariate normal
+kernel: h = Sigma*n^(-1/6), where Sigma = 0.5*(sd(x)+sd(y)), which supposes
+that the UD is bivariate normal.</p>
+
+<p> Alternatively, the smoothing parameter h may be computed by
+Least Square Cross Validation (LSCV). The estimated value then
+minimizes the Mean Integrated Square Error (MISE), i.e. the difference
+in volume between the true UD and the estimated UD. Note
+that the cross-validation criterion cannot be minimized in
+some cases (indicated by "convergence: no" by the module).
+According to Seaman and Powell (1998) <em>"This is a difficult problem
+that has not been worked out by statistical theoreticians, so no
+definitive response is available at this time"</em>
+(see Seaman and Powell, 1998 for further details and tricky solutions).
+
+<p>The flag -r returns the UD so that the contour of the equal value of the
+UD displayed by the module r.contour corresponds to the different
+percentage levels of home-range estimation.</p>
+
+<p>The results returned by this module are roughly identical to
+the results of the function kernelUD of the adehabitat package for the
+R software. The only difference is related to the LSCV method:
+because the minimization of the MISE requires numerical methods,
+the smoothing parameter found by v.adehabitat.kernelUD differ slightly from
+those returned by the function kernelUD in adehabitat (the correlation
+between the two is often greater than 0.99, pers.obs.).</p>
+
+<p>This module heavily relies on the code written by Stefano Menegon
+for his module v.kernel. It has just been extended to allow the LSCV
+and ad hoc estimation of the smoothing parameter, and the computation
+of the volume under the UD for home range estimation.</p>
+
+
+<h2>EXAMPLE</h2>
+
+<p>Estimation of the home range with the LSCV smoothing parameter:<br>
+
+<div class="code"><pre>
+v.adehabitat.kernelUD input=localisations output=ud -o
+</pre></div>
+</p>
+
+
+<p>
+Estimation of the UD with the ad hoc smoothing parameter:<br>
+
+<div class="code"><pre>
+v.adehabitat.kernelUD input=localisations output=ud
+</pre></div>
+</p>
+
+
+<p>
+Estimation of the volume under the UD with a LSCV smmoothing parameter:<br>
+
+<div class="code"><pre>
+v.adehabitat.kernelUD input=localisations output=ud -o -r
+</pre></div>
+</p>
+
+
+
+<H2>REFERENCES</H2>
+
+<p><EM>Worton, B.J. (1995) Using Monte Carlo simulation to evaluate
+kernel-based home range estimators. Journal of Wildlife Management,
+59, 794-800.</EM></p>
+
+<p><EM>Calenge, C. (2006) The package adehabitat for the R software: a tool
+for the analysis of space and habitat use by animals. Ecological
+Modelling, 197, 516-519.</EM></p>
+
+
+<H2>SEE ALSO</H2>
+<EM>
+<A HREF="v.kernel.html">v.kernel</A></EM>, <A HREF="v.adehabitat.clusthr.html">v.adehabitat.clusthr</A></EM>, <A HREF="v.adehabitat.mcp.html">v.adehabitat.mcp</A></EM>.
+
+
+
+<H2>AUTHOR</H2>
+Clement Calenge, Universite Lyon 1, France<br>
+Original code by Stefano Menegon, ITC-irst, Trento, Italy.
+
+<p><i>Last changed: $Date: 2007/08/23 17:13:12 $</i></p>
Added: grass-addons/vector/v.adehabitat.kernelUD/function.c
===================================================================
--- grass-addons/vector/v.adehabitat.kernelUD/function.c (rev 0)
+++ grass-addons/vector/v.adehabitat.kernelUD/function.c 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,42 @@
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include "global.h"
+
+
+
+/* probability for gaussian distribution */
+double gaussianKernel(double x, double term)
+{
+ return(term * exp(-(x*x)/2.));
+}
+
+
+
+
+double segno(double x)
+{
+ double y;
+
+ y = (x > 0 ? 1. : 0.) + (x < 0 ? -1. : 0.);
+
+ return y;
+}
+
+
+
+
+/* euclidean distance between vectors x and y of length n */
+double euclidean_distance(double *x, double *y, int n)
+{
+ int j;
+ double out = 0.0;
+ double tmp;
+
+ for(j=0;j<n;j++){
+ tmp = x[j] - y[j];
+ out += tmp * tmp;
+ }
+
+ return sqrt(out);
+}
Property changes on: grass-addons/vector/v.adehabitat.kernelUD/function.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/vector/v.adehabitat.kernelUD/global.h
===================================================================
--- grass-addons/vector/v.adehabitat.kernelUD/global.h (rev 0)
+++ grass-addons/vector/v.adehabitat.kernelUD/global.h 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,10 @@
+#include <grass/gis.h>
+
+double euclidean_distance(double *x, double *y, int n);
+double gaussianKernel(double x, double term);
+
+int read_points(struct Map_info *In, double ***coordinate);
+double compute_all_distances(double **coordinate, double **dists, int n, double dmax);
+void compute_distance( double N, double E, struct Map_info *In,
+ double sigma, double term, double *gaussian, double dmax);
+
Property changes on: grass-addons/vector/v.adehabitat.kernelUD/global.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/vector/v.adehabitat.kernelUD/main.c
===================================================================
--- grass-addons/vector/v.adehabitat.kernelUD/main.c (rev 0)
+++ grass-addons/vector/v.adehabitat.kernelUD/main.c 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,577 @@
+/****************************************************************************
+*
+* MODULE: v.kernelUD
+*
+* AUTHOR(S): Clement Calenge (univ. Lyon 1), original code by
+* Stefano Menegon, ITC-irst, Trento, Italy
+* PURPOSE: Generates a raster density map from vector points data using
+* a moving 2D isotropic Gaussian kernel
+* COPYRIGHT: (C) 2004 by the GRASS Development Team
+*
+* This program is free software under the GNU General Public
+* License (>=v2). Read the file COPYING that comes with GRASS
+* for details.
+*
+*****************************************************************************/
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/gmath.h>
+#include <grass/Vect.h>
+#include "global.h"
+
+
+static int ndists; /* number of distances in dists */
+static double *dists; /* array of all distances < dmax */
+static int npoints;
+static int verbose = 1 ;
+static double dimension = 2.;
+
+
+/* -----------------------------------------------------
+ Sources for memory allocation from ade4 package for R
+ -----------------------------------------------------*/
+
+void taballoc (double ***tab, int l1, int c1)
+
+{
+ int i, j;
+
+ if ( (*tab = (double **) calloc(l1+1, sizeof(double *))) != 0) {
+ for (i=0;i<=l1;i++) {
+ if ( (*(*tab+i)=(double *) calloc(c1+1, sizeof(double))) == 0 ) {
+ return;
+ for (j=0;j<i;j++) {
+ free(*(*tab+j));
+ }
+ }
+ }
+ }
+
+}
+
+void freetab (double **tab)
+{
+ int i, n;
+
+ n = *(*(tab));
+ for (i=0;i<=n;i++) {
+ free((char *) *(tab+i) );
+ }
+ free((char *) tab);
+}
+
+
+
+void vecintalloc (int **vec, int n)
+{
+ if ( (*vec = (int *) calloc(n+1, sizeof(int))) != NULL) {
+ **vec = n;
+ return;
+ } else {
+ return;
+ }
+}
+
+void freeintvec (int *vec)
+{
+
+ free((char *) vec);
+
+}
+
+
+/*--------------------------------------------------
+* End of sources of the ade4 package
+ --------------------------------------------------*/
+
+
+
+
+/* The score function for LSCV minimization: choice of the
+ best smoothing parameter */
+
+double L(double smooth)
+{
+ int ii;
+ double resL,n;
+
+ n = (double) npoints;
+ resL = 0.;
+
+ for(ii=0; ii < ndists; ii++){
+ resL+= (exp(-pow(dists[ii],2)/(4. * pow(smooth,2)))) - (4. * (exp(-pow(dists[ii],2)/(2. * pow(smooth,2.)))));
+ }
+
+ resL = 1./(M_PI * pow(smooth,2.) * n) + (2*resL -3*n)/(M_PI * 4. * pow(smooth,2.) * pow(n, 2.));
+
+ G_debug(3, "smooth = %e resL = %e", smooth, resL);
+ if(verbose)
+ G_message(_("\tScore Value=%f\tsmoothing parameter (standard deviation)=%f"),resL, smooth);
+
+ return(resL);
+}
+
+
+
+
+/* Main function for the kernelUD */
+
+int main(int argc, char **argv)
+{
+ /* structures of type "option" allow to specify the options of
+ the function (defined in gis.h) */
+ struct Option *in_opt, *out_opt;
+ struct Option *stddev_opt;
+
+ /* structures of type "flag" specify the flags of the function */
+ struct Flag *flag_o, *flag_v, *flag_q, *flag_h, *flag_r;
+
+ char *mapset;
+
+ /* structure "Map_info" contains the info concerning the
+ map formats (defined in dig_structs.i) */
+ struct Map_info In, Out;
+ int fdout = 0;
+ int row,col,ii,i;
+ struct Cell_head window; /* Info concerning the map header */
+ double gaussian;
+ double moyx = 0;
+ double moyy = 0;
+ double errtx = 0;
+ double errty = 0;
+ double cellsize;
+ double hvec;
+ double itt,tmple,tmpsco;
+ double N,E;
+ DCELL *output_cell = NULL;
+ double sigma, dmax;
+
+ double **coordinate;
+ double **tutu, **attribue;
+ int *colid, *rowid;
+ int rowb = 0;
+ int colb = 0;
+ double sigmaOptimal, cumsu;
+ struct GModule *module; /* infos related to the module
+ (name, description) */
+ double term;
+ double gausmax = 0;
+
+ /* Initialize the GIS calls */
+ G_gisinit(argv[0]); /* this function initialize the environment:
+ it verifies that the database and the
+ mapset are correct */
+
+ module = G_define_module(); /* define the info related to the module */
+ module->keywords = _("vector, kernel UD");
+ module->description =
+ _("Estimates the Utilization distribution on a raster map from vector points data using a moving 2D isotropic Gaussian kernel");
+
+
+/* Options of the function: */
+ in_opt = G_define_option();
+ in_opt->key = "input";
+ in_opt->type = TYPE_STRING;
+ in_opt->required = YES;
+ in_opt->description = _("input relocations");
+ in_opt->gisprompt = "old,vector,vector,input";
+
+ out_opt = G_define_option();
+ out_opt->key = "output";
+ out_opt->type = TYPE_STRING;
+ out_opt->required = YES;
+ out_opt->description = _("output raster map");
+
+ stddev_opt = G_define_option() ;
+ stddev_opt->key = "stddeviation";
+ stddev_opt->type = TYPE_DOUBLE;
+ stddev_opt->required = NO;
+ stddev_opt->description = _("Suggested smoothing parameter");
+
+/* Les flags */
+ flag_o = G_define_flag();
+ flag_o->key = 'o';
+ flag_o->description = _("LSCV smoothing parameter");
+
+ flag_q = G_define_flag();
+ flag_q->key = 'q';
+ flag_q->description = _("Only calculate LSCV smoothing parameter and exit (no map is written)");
+
+ flag_h = G_define_flag();
+ flag_h->key = 'h';
+ flag_h->description = _("Given smoothing parameter");
+
+ flag_v = G_define_flag();
+ flag_v->key = 'v';
+ flag_v->description = _("Run verbosely");
+
+ flag_r = G_define_flag();
+ flag_r->key = 'r';
+ flag_r->description = _("Compute home-range volume");
+
+
+ if (G_parser(argc, argv)) /* Any problem with the arguments? */
+ exit(EXIT_FAILURE);
+
+ /*read options*/
+ if ( flag_h->answer )
+ sigma = atof(stddev_opt->answer);
+ verbose = flag_v->answer;
+
+ if( flag_q->answer ) {
+ flag_o->answer=1;
+ }
+
+
+ /* Get the region */
+ G_get_window(&window);
+
+ G_message("RES: %f\tROWS: %d\tCOLS: %d",
+ window.ew_res, window.rows, window.cols);
+ cellsize= pow(window.ew_res, 2);
+
+ /* Open input vector */
+ /* G_find_vector2 finds a vector map (look but dont touch)*/
+ if ((mapset = G_find_vector2 (in_opt->answer, "")) == NULL)
+ G_fatal_error (_("Could not find input map '%s'."), in_opt->answer);
+
+ Vect_set_open_level (2); /* Essential before opening a map!
+ defines the level */
+ Vect_open_old (&In, in_opt->answer, mapset); /* Opens the vector map */
+
+ /* check and open the name of output map */
+ if( !flag_q->answer ) {
+ if(G_legal_filename( out_opt->answer ) < 0)
+ G_fatal_error(_("illegal file name [%s]."), out_opt->answer);
+
+ G_set_fp_type (DCELL_TYPE); /* storage type for floating-point maps. */
+
+ /* Open the output file, and return a descriptor
+ of the file in fdout */
+ if((fdout = G_open_raster_new(out_opt->answer,DCELL_TYPE)) < 0)
+ G_fatal_error(_("error opening raster map [%s]."), out_opt->answer);
+
+
+ /* allocate output raster */
+ output_cell=G_allocate_raster_buf(DCELL_TYPE);
+ }
+
+ /* Read points */
+ npoints = read_points ( &In, &coordinate );
+
+
+ if (!(flag_h->answer)) {
+
+ /* average of x and y */
+ for (ii = 0; ii < npoints; ii++) {
+ moyx += (coordinate[ii][0] / ((double) npoints));
+ moyy += (coordinate[ii][1] / ((double) npoints));
+ }
+
+ /* variance of x and y */
+ for (ii = 0; ii < npoints; ii++) {
+ errtx += ((coordinate[ii][0] - moyx) *
+ (coordinate[ii][0] - moyx)) / ((double) (npoints-1));
+ errty += ((coordinate[ii][1] - moyy) *
+ (coordinate[ii][1] - moyy)) / ((double) (npoints-1));
+ }
+ /* Compute href in sigma */
+ sigma = ((sqrt((errtx+errty)/2)) * pow( ((double) npoints), -1./6.));
+
+ /* maximum distance 4*sigma (3.9*sigma ~ 1.0000),
+ * keep it small, otherwise it takes
+ * too much points and calculation on network becomes slow */
+ }
+
+ dmax = 400*sigma; /* used as maximum value */
+
+ G_message(_("Using maximum distance between points: %f"), dmax);
+
+ ndists = compute_all_distances(coordinate,&dists,npoints,dmax);
+ /* here, one computes all the distances LOWER THAN dmax */
+
+ G_message(_("Number of input points: %d."), npoints);
+ G_message(_("%d distances read from the map."), ndists);
+
+ if (ndists == 0)
+ G_fatal_error(_("distances between all points are beyond %e (4 * "
+ "standard deviation) cannot calculate optimal value."), dmax);
+
+
+ /* Computes the optimal smoothing parameter by LSCV */
+ if( flag_o->answer ) {
+
+ tmple=(( (1.5 * sigma) - (0.1 * sigma) )/100.);
+ tmpsco = L((0.1 * sigma));
+ sigmaOptimal=(0.1 * sigma);
+
+ for (itt=(0.1 * sigma); itt<(1.5 * sigma); itt+=tmple) {
+ hvec=L(itt);
+ if (hvec < tmpsco) {
+ sigmaOptimal=itt;
+ tmpsco=hvec;
+ }
+ }
+ /* Reset sigma to calculated optimal value */
+ tmpsco=(0.1 * sigma)+tmple;
+ if (sigmaOptimal<tmpsco)
+ G_message(_("Convergence: no"));
+ if (sigmaOptimal>=tmpsco)
+ G_message(_("Convergence: yes"));
+ sigma=sigmaOptimal;
+ }
+ G_message(_("Smoothing parameter: %f"), sigma);
+
+
+ /* in case just LSCV h is desired */
+ if( flag_q->answer ) {
+ Vect_close (&In);
+ exit (0);
+ }
+
+
+ /* Now, compute the UD */
+
+ term=1./(pow(sigma,2.) * 2. * M_PI * npoints);
+ dmax= sigma*4.;
+
+
+ /* case where the volume is not desired */
+ if (!(flag_r->answer)) {
+ G_message(_("\nComputing UD using smoothing parameter=%f."), sigma);
+ for(row=0; row<window.rows; row++){
+ G_percent(row,window.rows,2);
+
+ for(col=0; col<window.cols; col++) {
+
+ N = G_row_to_northing(row+0.5,&window);
+ E = G_col_to_easting(col+0.5,&window);
+
+ compute_distance ( N, E, &In, sigma, term, &gaussian,dmax );
+ output_cell[col] = gaussian;
+ if ( gaussian > gausmax ) gausmax = gaussian;
+ }
+ G_put_raster_row(fdout,output_cell,DCELL_TYPE);
+ }
+ G_close_cell(fdout);
+ }
+
+ /* case the volume is desired */
+ taballoc(&tutu, window.rows, window.cols);
+ taballoc(&attribue, window.rows, window.cols);
+ vecintalloc(&rowid, (window.rows * window.cols));
+ vecintalloc(&colid, (window.rows * window.cols));
+
+ if (flag_r->answer) {
+ G_message(_("\nComputing UD using smoothing parameter=%f."), sigma);
+
+ /* computes the UD */
+ for(row=0; row<window.rows; row++){
+ G_percent(row,window.rows,2);
+
+ for(col=0; col<window.cols; col++) {
+
+ N = G_row_to_northing(row+0.5,&window);
+ E = G_col_to_easting(col+0.5,&window);
+
+ compute_distance ( N, E, &In, sigma, term, &gaussian,dmax );
+ ii=0;
+ tutu[row+1][col+1]= gaussian;
+ }
+ }
+
+ /* compute the volume */
+ cumsu=0;
+ for(row=0; row<window.rows; row++){
+ for(col=0; col<window.cols; col++) {
+ cumsu+=tutu[row+1][col+1];
+ }
+ }
+
+ G_message(_("\nVolume under the UD=%f."), cumsu*cellsize);
+
+ /* standardization of the volume */
+ for(row=0; row<window.rows; row++){
+ for(col=0; col<window.cols; col++) {
+ tutu[row+1][col+1]=(tutu[row+1][col+1] / cumsu);
+ }
+ }
+
+ /* Sort the values */
+ G_message(_("\nComputing home ranges..."));
+ ii=0;
+ for(row=0; row<window.rows; row++){
+ G_percent(row,window.rows,2);
+ for(col=0; col<window.cols; col++) {
+
+ tmple=-1;
+ for (rowb=0; rowb<window.rows; rowb++){
+ for (colb=0; colb<window.cols; colb++) {
+ if (tutu[rowb+1][colb+1] >= tmple) {
+ if (abs(attribue[rowb+1][colb+1])<0.5) {
+ tmple=tutu[rowb+1][colb+1];
+ rowid[ii+1]=rowb;
+ colid[ii+1]=colb;
+ }
+ }
+ }
+ }
+ attribue[rowid[ii+1]+1][colid[ii+1]+1]=1;
+ ii++;
+ }
+ }
+
+ /* set the levels of the corresponding home range */
+ G_message(_("\nOutput"));
+ cumsu=0;
+ for (i=0; i < ii; i++) {
+ if (cumsu<0.9999999) {
+ cumsu+=tutu[rowid[i+1]+1][colid[i+1]+1];
+ } else {
+ cumsu=1;
+ }
+ tutu[rowid[i+1]+1][colid[i+1]+1]=cumsu;
+
+ }
+
+
+ /* output */
+ for(row=0; row<window.rows; row++){
+ for(col=0; col<window.cols; col++) {
+ output_cell[col] = tutu[row+1][col+1];
+ }
+ G_put_raster_row(fdout,output_cell,DCELL_TYPE);
+ }
+
+ G_close_cell(fdout);
+
+ }
+
+
+ Vect_close (&In);
+ freetab(tutu);
+ freetab(attribue);
+ freeintvec(rowid);
+ freeintvec(colid);
+ exit(0);
+}
+
+
+/* Read points to array return number of points */
+int read_points( struct Map_info *In, double ***coordinate)
+{
+ int line, nlines, npoints, ltype, i = 0;
+ double **xySites;
+ static struct line_pnts *Points = NULL;
+
+ if (!Points)
+ Points = Vect_new_line_struct ();
+
+ /* Allocate array of pointers */
+ npoints = Vect_get_num_primitives(In,GV_POINT);
+ xySites = (double **) G_calloc ( npoints, sizeof(double*) );
+
+ nlines = Vect_get_num_lines(In);
+
+ for ( line = 1; line <= nlines; line++){
+ ltype = Vect_read_line (In, Points, NULL, line);
+ if ( !(ltype & GV_POINT ) ) continue;
+
+ xySites[i] = (double *)G_calloc((size_t)2, sizeof(double));
+
+ xySites[i][0] = Points->x[0];
+ xySites[i][1] = Points->y[0];
+ i++;
+ }
+
+ *coordinate = xySites;
+
+ return (npoints);
+}
+
+
+/* Calculate distances < dmax between all sites in coordinate
+ * Return: number of distances in dists */
+double compute_all_distances(double **coordinate, double **dists, int n, double dmax)
+{
+ int ii,jj,kk;
+ size_t nn;
+
+ nn = n*(n-1)/2;
+ *dists = (double *) G_calloc(nn,sizeof(double));
+ kk=0;
+
+ for(ii=0; ii < n-1; ii++){
+ for(jj=ii+1; jj<n; jj++){
+ double dist;
+
+ dist = euclidean_distance(coordinate[ii],coordinate[jj],2);
+ G_debug (3, "dist = %f", dist);
+
+ if ( dist <= dmax ) {
+ (*dists)[kk] = dist;
+ kk++;
+ }
+ }
+ }
+
+ return (kk);
+}
+
+
+
+void compute_distance( double N, double E, struct Map_info *In,
+ double sigma, double term, double *gaussian, double dmax)
+{
+ int line, nlines;
+ double a[2],b[2];
+ double dist;
+
+ /* spatial index handling, borrowed from lib/vector/Vlib/find.c */
+ BOUND_BOX box;
+ static struct ilist *NList = NULL;
+ static struct line_pnts *Points = NULL;
+
+ a[0] = E;
+ a[1] = N;
+
+ if (!NList)
+ NList = Vect_new_list ();
+ if (!Points)
+ Points = Vect_new_line_struct ();
+
+ /* create bounding box 2x2*dmax size from the current cell center */
+ box.N = N + dmax;
+ box.S = N - dmax;
+ box.E = E + dmax;
+ box.W = E - dmax;
+ box.T = HUGE_VAL;
+ box.B = -HUGE_VAL;
+
+ /* number of lines within dmax box */
+ nlines = Vect_select_lines_by_box (In, &box, GV_POINT, NList);
+
+ *gaussian=.0;
+
+ for ( line = 0; line < nlines; line++ ) {
+
+ Vect_read_line (In, Points, NULL, NList->value[line]);
+
+ b[0] = Points->x[0];
+ b[1] = Points->y[0];
+
+ dist = euclidean_distance(a,b,2);
+
+ if(dist<=dmax)
+ *gaussian += gaussianKernel(dist/sigma,term);
+
+ }
+}
+
+
+
+
Property changes on: grass-addons/vector/v.adehabitat.kernelUD/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/vector/v.adehabitat.mcp/Makefile
===================================================================
--- grass-addons/vector/v.adehabitat.mcp/Makefile (rev 0)
+++ grass-addons/vector/v.adehabitat.mcp/Makefile 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.adehabitat.mcp
+
+LIBES = $(SITESLIB) $(VECTLIB) $(GISLIB) $(CURSES)
+DEPENDENCIES= $(VECTDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/vector/v.adehabitat.mcp/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/vector/v.adehabitat.mcp/description.html
===================================================================
--- grass-addons/vector/v.adehabitat.mcp/description.html (rev 0)
+++ grass-addons/vector/v.adehabitat.mcp/description.html 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,40 @@
+<h2>DESCRIPTION</h2>
+
+<p><em>v.adehabitat.mcp</em> computes the home range of one animal using
+the Minimum Convex Polygon estimator. This function computes the
+Minimum Convex Polygon estimation after the removal of '(100 minus percent)'
+percent of the relocations the farthest away from the centroid of the
+home range (computed by the arithmetic mean of the coordinates of
+the relocations for each animal). This function is just a slight
+modification of the module v.hull of Andrea Aime (aaime at libero.it).</p>
+
+
+<H2>REFERENCES</H2>
+
+<p><EM> Mohr, C.O. (1947) Table of equivalent populations of north
+american small mammals. The American Midland Naturalist, 37,
+223-249.</EM></p>
+
+
+
+<h2>EXAMPLE</h2>
+
+Estimation of the 95% home range:<br>
+
+<div class="code"><pre>
+v.adehabitat.mcp input=localisations output=homerange percent=95
+</pre></div>
+
+
+<h2>SEE ALSO</h2>
+
+<em><a HREF="v.hull.html">v.hull</a></em>,
+<em><a HREF="v.adehabitat.clusthr.html">v.adehabitat.clusthr</a></em>,
+<em><a HREF="v.adehabitat.kernelUD.html">v.adehabitat.kernelUD</a></em>
+
+
+<h2>AUTHOR</h2>
+
+Clement Calenge
+
+<p><i>Last changed: $Date: 2007/08/23 17:12:30 $</i>
Added: grass-addons/vector/v.adehabitat.mcp/main.c
===================================================================
--- grass-addons/vector/v.adehabitat.mcp/main.c (rev 0)
+++ grass-addons/vector/v.adehabitat.mcp/main.c 2008-08-28 07:02:08 UTC (rev 33136)
@@ -0,0 +1,414 @@
+/******************************************************************************
+*
+* Computes the minimum convex polygon from a file of relocations
+*
+* 29 August 2006
+*
+* Clement Calenge, Original code from Andrea Aime
+* This file is part of GRASS GIS. It 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.
+******************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/site.h>
+#include <grass/glocale.h>
+
+
+struct Point {
+ double x;
+ double y;
+};
+
+void vecalloc (double **vec, int n)
+/*--------------------------------------------------
+ * Memory Allocation for a vector of length n
+ --------------------------------------------------*/
+{
+ if ( (*vec = (double *) calloc(n+1, sizeof(double))) != 0) {
+ **vec = n;
+ return;
+ } else {
+ return;
+ }
+}
+
+/*****************/
+void vecintalloc (int **vec, int n)
+/*--------------------------------------------------
+ * Memory allocation for an integer vector of length n
+ --------------------------------------------------*/
+{
+ if ( (*vec = (int *) calloc(n+1, sizeof(int))) != NULL) {
+ **vec = n;
+ return;
+ } else {
+ return;
+ }
+}
+
+
+/***********************************************************************/
+void freevec (double *vec)
+/*--------------------------------------------------
+ * Free memory for a vector
+ --------------------------------------------------*/
+{
+ free((char *) vec);
+}
+
+/***********************************************************************/
+void freeintvec (int *vec)
+/*--------------------------------------------------
+* Free memory for an integer vector
+--------------------------------------------------*/
+{
+
+ free((char *) vec);
+
+}
+
+
+
+int read_points( struct Map_info *In, double ***coordinate)
+{
+ int line, nlines, npoints, ltype, i = 0;
+ double **xySites;
+ static struct line_pnts *Points = NULL;
+
+ if (!Points)
+ Points = Vect_new_line_struct ();
+
+ /* Allocate array of pointers */
+ npoints = Vect_get_num_primitives(In,GV_POINT);
+ xySites = (double **) G_calloc ( npoints, sizeof(double*) );
+
+ nlines = Vect_get_num_lines(In);
+
+ for ( line = 1; line <= nlines; line++){
+ ltype = Vect_read_line (In, Points, NULL, line);
+ if ( !(ltype & GV_POINT ) ) continue;
+
+ xySites[i] = (double *)G_calloc((size_t)2, sizeof(double));
+
+ xySites[i][0] = Points->x[0];
+ xySites[i][1] = Points->y[0];
+ i++;
+ }
+
+ *coordinate = xySites;
+
+ return (npoints);
+}
+
+
+
+int rightTurn(struct Point *P, int i, int j, int k) {
+ double a, b, c, d;
+ a = P[i].x - P[j].x;
+ b = P[i].y - P[j].y;
+ c = P[k].x - P[j].x;
+ d = P[k].y - P[j].y;
+ return a*d - b*c < 0;
+}
+
+
+int cmpPoints(const void* v1, const void* v2) {
+ struct Point *p1, *p2;
+ p1 = (struct Point*) v1;
+ p2 = (struct Point*) v2;
+ if( p1->x > p2->x )
+ return 1;
+ else if( p1->x < p2->x )
+ return -1;
+ else
+ return 0;
+}
+
+
+
+
+int convexHull(struct Point* P, const int numPoints, int **hull) {
+ int pointIdx, upPoints, loPoints;
+ int *upHull, *loHull;
+
+ /* sort points in ascending x order*/
+ qsort(P, numPoints, sizeof(struct Point), cmpPoints);
+
+ *hull = (int*) G_malloc(numPoints * 2 * sizeof(int));
+
+ /* compute upper hull */
+ upHull = *hull;
+ upHull[0] = 0;
+ upHull[1] = 1;
+ upPoints = 1;
+ for(pointIdx = 2; pointIdx < numPoints; pointIdx++) {
+ upPoints++;
+ upHull[upPoints] = pointIdx;
+ while( upPoints > 1 &&
+ !rightTurn(P, upHull[upPoints], upHull[upPoints-1],
+ upHull[upPoints-2])
+ ) {
+ upHull[upPoints-1] = upHull[upPoints];
+ upPoints--;
+ }
+ }
+
+ /*
+ printf("upPoints: %d\n", upPoints);
+ for(pointIdx = 0; pointIdx <= upPoints; pointIdx ++)
+ printf("%d ", upHull[pointIdx]);
+ printf("\n");
+ */
+
+ /* compute lower hull, overwrite last point of upper hull */
+ loHull = &(upHull[upPoints]);
+ loHull[0] = numPoints - 1;
+ loHull[1] = numPoints - 2;
+ loPoints = 1;
+ for(pointIdx = numPoints - 3; pointIdx >= 0; pointIdx--) {
+ loPoints++;
+ loHull[loPoints] = pointIdx;
+ while( loPoints > 1 &&
+ !rightTurn(P, loHull[loPoints], loHull[loPoints-1],
+ loHull[loPoints-2])
+ ) {
+ loHull[loPoints-1] = loHull[loPoints];
+ loPoints--;
+ }
+ }
+
+ G_debug(3, "numPoints:%d loPoints:%d upPoints:%d",
+ numPoints, loPoints, upPoints);
+ /*
+ printf("loPoints: %d\n", loPoints);
+ for(pointIdx = 0; pointIdx <= loPoints; pointIdx ++)
+ printf("%d ", loHull[pointIdx]);
+ printf("\n");
+ */
+
+ /* reclaim uneeded memory */
+ *hull = (int *) G_realloc(*hull, (loPoints + upPoints) * sizeof(int));
+ return (loPoints + upPoints);
+}
+
+
+
+
+/*
+ * Outputs the points that comprises the convex hull as a single closed line
+ * and the hull baricenter as the label points (as it is a linear combination
+ * of points on the hull is guaranteed to be inside the hull, follow from the
+ * definition of convex polygon)
+ */
+
+int outputHull(struct Map_info *Map, struct Point* P, int *hull,
+ int numPoints) {
+ struct line_pnts *Points;
+ struct line_cats *Cats;
+ double *tmpx, *tmpy;
+ int i, pointIdx;
+ double xc, yc;
+
+ tmpx = (double *) G_malloc((numPoints + 1) * sizeof(double));
+ tmpy = (double *) G_malloc((numPoints + 1) * sizeof(double));
+
+ xc = yc = 0;
+ for(i = 0; i < numPoints; i++) {
+ pointIdx = hull[i];
+ tmpx[i] = P[pointIdx].x;
+ tmpy[i] = P[pointIdx].y;
+ /* average coordinates calculation... may introduce a little
+ numerical error but guaratees that no overflow will occurr */
+ xc = xc + tmpx[i] / numPoints;
+ yc = yc + tmpy[i] / numPoints;
+ }
+ tmpx[numPoints] = P[hull[0]].x;
+ tmpy[numPoints] = P[hull[0]].y;
+
+ Points = Vect_new_line_struct ();
+ Cats = Vect_new_cats_struct ();
+ Vect_copy_xyz_to_pnts(Points, tmpx, tmpy, 0, numPoints+1);
+ G_free(tmpx);
+ G_free(tmpy);
+
+ /* write out convex hull */
+ Vect_write_line (Map, GV_BOUNDARY, Points, Cats);
+
+ /* find and add centroid */
+ Vect_reset_line (Points);
+ Vect_append_point (Points,xc, yc, 0.0);
+ Vect_cat_set ( Cats, 1, 1 );
+ Vect_write_line (Map, GV_CENTROID, Points, Cats);
+ Vect_destroy_line_struct (Points);
+
+
+ return 0;
+}
+
+
+
+int main(int argc, char **argv) {
+ struct GModule *module;
+ struct Option *input, *output, *percent;
+ struct Flag *all;
+ struct Cell_head window;
+
+ char *mapset;
+ FILE* fdsite;
+ struct Map_info Map, In;
+ struct Point *points; /* point loaded from site file */
+ struct Point *pointsb; /* point after removing outliers */
+ int *hull; /* index of points located on the convex hull */
+ int numSitePoints, numHullPoints;
+ int *cons;
+ int nnn, pascons, ran, kk, i, j;
+ double *dist, yc, xc, perc, tmpd;
+ double **coordinate;
+
+ G_gisinit (argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("vector, geometry");
+ module->description = "Uses a GRASS vector points map to compute the Minimum convex polygon home range";
+
+ input = G_define_option ();
+ input->key = "input";
+ input->type = TYPE_STRING;
+ input->required = YES;
+ input->description = "Name of map containing the relocations to be input";
+ input->gisprompt = "old,vector,vector,input";
+
+ output = G_define_option ();
+ output->key = "output";
+ output->type = TYPE_STRING;
+ output->required = YES;
+ output->description = "Name of a vector area map to be output";
+ output->gisprompt = "new,dig,binary file,output";
+
+ percent = G_define_option ();
+ percent->key = "percent";
+ percent->type = TYPE_STRING;
+ percent->required = YES;
+ percent->description = "Percentage of the relocations to keep";
+ percent->answer = "95";
+
+ if (G_parser (argc, argv)) exit (1);
+
+ Vect_check_input_output_name ( input->answer, output->answer, GV_FATAL_EXIT );
+
+ /* récupère la région */
+ G_get_window (&window);
+
+
+ /* open relocations file */
+ if ((mapset = G_find_vector2 (input->answer, "")) == NULL)
+ G_fatal_error (_("Could not find input map '%s'."), input->answer);
+
+ Vect_set_open_level (2); /* Essential before opening a map!
+ defines the level */
+ Vect_open_old (&In, input->answer, mapset); /* Opens the vector map */
+
+
+ /* checks for legal file name */
+ if(G_legal_filename( output->answer ) < 0)
+ G_fatal_error(_("illegal file name [%s]."), output->answer);
+
+
+ /* load coordinates */
+ numSitePoints = read_points ( &In, &coordinate );
+ if(numSitePoints < 3 )
+ G_fatal_error ("MCP calculation requires at least three points");
+
+ /* creates the points data structure */
+ points = (struct Point *) G_malloc(numSitePoints * sizeof(struct Point));
+ for (i = 0; i < numSitePoints; i++) {
+ points[i].x = coordinate[i][0];
+ points[i].y = coordinate[i][1];
+ }
+
+
+ /* load the percent */
+ perc = atof(percent->answer);
+
+ /* create vector file */
+ if (0 > Vect_open_new (&Map, output->answer, 0) )
+ G_fatal_error ("Unable to open vector file <%s>\n", output->answer);
+
+ Vect_hist_command ( &Map );
+
+ /* Compute the barycenter */
+
+ vecalloc(&dist, numSitePoints);
+ xc=0;
+ yc=0;
+ for(i = 0; i < numSitePoints; i++) {
+ xc = xc + (points[i].x) / ((double) numSitePoints);
+ yc = yc + (points[i].y) / ((double) numSitePoints);
+ }
+
+ for(i = 0; i < numSitePoints; i++) {
+ dist[i+1] = sqrt(pow((points[i].x -xc),2) + pow((points[i].y - yc),2));
+ }
+
+ nnn = ((int) ((double) numSitePoints - (perc * ((double) numSitePoints) / 100)));
+ vecintalloc(&cons, nnn);
+
+ kk = 1;
+ for (i = 1; i <= numSitePoints; i++) {
+ tmpd = dist[i];
+ ran = 0;
+ for (j = 1; j <= numSitePoints; j++) {
+ if (dist[j] > tmpd)
+ ran++;
+ }
+ if (ran < nnn) {
+ cons[kk] = i;
+ kk++;
+ }
+ }
+
+ pointsb = (struct Point *) G_malloc((numSitePoints - nnn) * sizeof(struct Point));
+ kk = 0;
+ for (i = 1; i <= numSitePoints; i++) {
+ pascons = 0;
+ for (j = 1; j <= nnn; j++) {
+ if (i==cons[j])
+ pascons = 1;
+ }
+ if (pascons == 0) {
+ pointsb[kk].x = points[i-1].x;
+ pointsb[kk].y = points[i-1].y;
+ kk++;
+ }
+ }
+ numSitePoints = numSitePoints - nnn;
+ if(numSitePoints < 3 )
+ G_fatal_error ("MCP calculation requires at least three points");
+
+
+ /* compute convex hull */
+ numHullPoints = convexHull(pointsb, numSitePoints, &hull);
+
+ /* output vector file */
+ outputHull(&Map, pointsb, hull, numHullPoints);
+
+ /* clean up and bye bye */
+ Vect_build (&Map, stdout);
+ Vect_close (&Map);
+
+ freevec(dist);
+ freeintvec(cons);
+ return 0;
+}
Property changes on: grass-addons/vector/v.adehabitat.mcp/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
More information about the grass-commit
mailing list