/* 
   disum1.c
   Implement DSF sum with voxel image size = 1
*/

#include "stdio.h"
#include "math.h"

long disum(argc,argv)
int argc;
void *argv[];

{
long   N;
double *dat;
double *i;
double *j;
double *ff;

long k;

long   i0,j0,id,jd;
double di,dj,adi,adj,a0i,a0j;

N   = *(long *)argv[0];   /* Number of pixels in a data array edge */
dat = (double *)argv[1];  /* data array */
i   = (double *)argv[2];  /* i index */
j   = (double *)argv[3];  /* j index */
ff  = (double *)argv[4];  /* multiplier */

for(k=0;k<N*N;k++)
{

  i0  = floor(i[k]+0.5);
  di  = i[k]-i0;
  adi = fabs(di);

  j0  = floor(j[k]+0.5);
  dj  = j[k]-j0;
  adj = fabs(dj);

  id  = i0+((di>0)?1:-1);  /* (di>0)?1:-1 is the sign of di */
  jd  = j0+((dj>0)?1:-1);  /* (dj>0)?1:-1 is the sign of dj */
  a0i = 1-adi;
  a0j = 1-adj;

/* Eliminate spreading
a0i=1;
a0j=1;
adi=0;
adj=0;
*/

  if((i0 >= 0) && (i0 < N))
  { if((j0 >= 0) && (j0 < N))dat[N*j0+i0] += ff[k]*a0i*a0j;
    if((jd >= 0) && (jd < N))dat[N*jd+i0] += ff[k]*a0i*adj;
  }

  if((id >= 0) && (id < N))
  { if((j0 >= 0) && (j0 < N))dat[N*j0+id] += ff[k]*adi*a0j;
    if((jd >= 0) && (jd < N))dat[N*jd+id] += ff[k]*adi*adj;
  }

/*  debug
if((i0==28) && (j0==42))
{ printf("\n1: i0=%d i=%f di=%f adi=%f si=%d id=%d a0i=%f",i0,i[k],di,adi,si,id,a0i);
  printf("\n   j0=%d j=%f dj=%f adj=%f sj=%d jd=%d a0j=%f",j0,j[k],dj,adj,sj,jd,a0j);
  printf("\n   dat[N*j0+i0]=%f dat[N*j0+id]=%f dat[N*jd+i0]=%f dat[N*jd+id]=%fi\n",
  dat[N*j0+i0],dat[N*j0+id],dat[N*jd+i0],dat[N*jd+id]);
  return(N);
}
*/

/*
if((a0i<0)||(a0j<0)||(adi<0)||(adj<0))
printf("\n disum.c error: %f %f %d %d",i[k],j[k],i0,j0);
*/

}

return(10);

}