Revision 13a20a05ba53faf0cbd40e9aefdcf87f9ca6d450 authored by David J. Marchette on 11 November 2010, 00:00:00 UTC, committed by Gabor Csardi on 11 November 2010, 00:00:00 UTC
1 parent c8cf498
cccd.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <float.h>
#define MIN(x,y) (((x)<(y))?(x):(y))
#ifndef MAXFLOAT
#define MAXFLOAT 10E20
#endif
/* sorts smallest to biggest */
extern int cmp1(const void *x, const void *y);
/*
* DX is nx by nx
* DY is nx by ny
*
* Returns
* radii in R,
* covered points in C.
*/
void cccd(double *DX, /* distances between x and x */
double *DY, /* distances between x and y */
int *NX, /* number of x points */
int *NY, /* number of y points */
double *R, /* radiuses (returned) */
int *C, /* points covered (returned) */
int *K, /* measure of y covered, AKA Beta */
int *SHRINK, /* shrink the radii back to first covered x */
int *VERBOSE)
{
int i,j,k;
int nx=*NX;
int ny=*NY;
double *dy;
int *covered;
double *p;
int verbose=*VERBOSE;
double maxxd;
double temp;
int shrink=*SHRINK;
k = *K;
if(verbose) fprintf(stderr,"running cccd\n");
dy = (double *)calloc(ny,sizeof(double));
p = DY;
for(i=0;i<nx;i++){
covered = C+i*nx;
memset(covered,0,sizeof(int)*nx);
memcpy(dy,p,ny*sizeof(double));
qsort(dy,ny,sizeof(double),cmp1);
if(k>=ny) {
R[i] = MAXFLOAT;
fprintf(stderr,"\n%g",R[i]);
}
else{
R[i] = dy[k];
}
maxxd = dy[0];
/* find the farthest x covered by the ball */
for(j=0;j<nx;j++){
temp = *(DX+i*nx+j);
if(temp<R[i]){
maxxd = (maxxd>temp)?maxxd:temp;
covered[j] = 1;
}
}
/* shrink the ball down so that it covers only those x's */
if(shrink){
for(j=0;j<k+1;j++){
if(dy[j]>=maxxd) {
maxxd = dy[j];
break;
}
}
/* shrink the ball so it covers only those x's */
R[i] = maxxd;
}
p += ny;
}
if(verbose) fprintf(stderr,"done\n");
free(dy);
}
Computing file changes ...