/* Last edited on 2021-07-13 21:14:14 by jstolfi */

/* INDEX ROLL OPERATION - CORRECT BUT TOO CONFUSING TO BE USEFUL */

void ppv_roll ( ppv_array_t *A, ppv_axis_t i, ppv_axis_t j, int n );
  /* Modifies the descriptor {A} so that it describes the same samples
    as before, with indices {i..j} cyclically shifted by {n} places
    forward. Thus, element {A[c,h,v,d,t,u]} before the call
    {ppv_roll(A, 2,4, 1)} is equivalent to element {A[c,h,t,v,d,u]}
    after the call. 
    
    If the shift {n} is negative, or greater than {M=j-i+1}, it is
    reduced modulo {M} to the range {0..j-i}. */

void ppv_roll ( ppv_array_t *A, ppv_axis_t i, ppv_axis_t j, int n )
  { ix_roll(N, A->size, &(A->base), A->step, i, j, n); }

/* ERROR MESSAGES */

void ppv_progerror 
  ( const char *msg, 
    const char *file, 
    const unsigned int line, 
    const char* proc 
  );
  /* Prints {file ":" line ": (" *proc ")" *msg} to {stderr} and stops. */

#define ppv_error(msg) \
  ppv_progerror((msg), __FILE__, __LINE__, __FUNCTION__)
  /* Prints {*msg} and the source location.  */

  
void ppv_progerror 
  ( const char *msg, 
    const char *file, 
    const unsigned int line, 
    const char* proc 
  )
  { fprintf (stderr, "%s:%u: *** (%s) %s\n", file, line, proc, msg);
    exit(1);
  }

   for (ppv_index_t r0 = 0;  r0 < sz; r0++)
      { ppv_index_t d0 = r0 - b->ix0min;
        /* Create slice number {r0} of brush: voxels with {ix[0]=d0}. */ 
        ppv_sample_count_t nvr; /* Number of voxels in slice. */
        collect_slice_voxels(d0, &(b->ix1[r0]), &nvr); 
        if (debug) { fprintf(stderr, " %lu", nvr); }

        b->step[r0] = notnull(malloc(nvr*sizeof(ppv_step_t)), "no mem");
        for (ppv_sample_count_t k = 0; k < nvr; k++) { b->step[r0][k] = 0; }

        b->nv[r0] = nvr;
        nvtot += nvr;
      }

/* The vector {b.step} is an auxiliary area used to speed up local
  operations on a specific {ppv_array_desc_t} {A}. Element
  {b.step[r0][k]} is the relative position (bit address, {ppv_pos_t})
  in {A} of the neighborhood voxel that coresponds to voxel {V[k]} of
  slice {r0} of the brush. If {ixc[0..d-1]} is the index vector of the
  center of the neighborhood in {A}, the increment {b.step[r0][k]}
  will be relative to the position of the neighborhood voxel with
  indices {ixv[0..d-1]} where {ixv[0] = ixc[0]+r0-b.ix0min} and
  {ixv[ax] = ixc[ax]} for {ax} in {1..d-1}. Note that the displacement
  {b.step[r0][k]} is within the slice of {A} that contains that
  voxel. */

    auto void collect_voxels(ppv_index_t d0, ppv_index_t **ix1rP, ppv_sample_count_t *nvrP);
      /* Finds all voxels of the ball that have index {ix[0]} equal to {d0}.
        Returns in {*ix1rP} a newly allocated vector that has ix1ices {1..d-1}
        of all those voxels. Returns in {*nvrP} the count of such voxels.
        Requres {iabs(d0) <= rad}. */
    
    void collect_slice_voxels(ppv_index_t d0, ppv_index_t **ix1rP, ppv_sample_count_t *nvrP)
      { 
        assert(labs(d0) <= rad);
        double radr = sqrt(rad*rad - (double)(d0*d0)) + 1.0e-8; /* Radius of ball on this slice. */
        ppv_size_t iradr = (ppv_size_t)floor(fabs(radr)); /* Max radius of digital ball on this slice. */
        ppv_size_t szr = 2*iradr + 1; /* Max count of voxels on this slice along axes {1..d-1}. */
        ppv_sample_count_t nvmaxr = ipow(szr, d-1); /* Max total count of voxels on this slice. */
        ppv_index_t *ix1r = notnull(malloc(nvmaxr*(d-1)*sizeof(ppv_index_t)), "no mem");
        ppv_sample_count_t nvr = 0; /* Number of voxels in slice. */
        /* Now we enumerate all index tuples with {ix[0]} fixed at {r0}
          and the other indices ranging in {0..szr-1}. */
        ppv_index_t ix[d]; /* Indices of a brush voxel in slice, shifted. */
        ppv_size_t sz[d]; /* Max voxel count along each axis in slice. */
        ix[0] = 0; sz[0] = 1;
        for (ppv_axis_t ax = 1; ax < d; ax++) { ix[ax] = 0; sz[ax] = szr; }
        bool_t done = FALSE; /* Exhausted all possible voxel indices in slice. */
        for (ppv_sample_count_t kv = 0; kv < nvmaxr; kv++)
          { assert(! done);
            double dist2 = (double)(d0*d0); /* Squared distance from center voxel. */
            for (ppv_axis_t ax = 1; ax < d; ax++) 
              { ppv_index_t da = ix[ax]-iradr; dist2 += (double)(da*da); }
            if (dist2 <= rad*rad) 
              { assert(nvr < nvmaxr);
                for (ppv_index_t ax = 1; ax < d; ax++)
                  { ix1r[nvr*(d-1) + ax] = ix[ax] - iradr;
                nvr++;
              }
             }
             done = ix_next(d, ix, sz, ix_order_L, NULL, NULL, NULL, NULL, NULL, NULL);
           }
        assert(done);
        assert(nvr > 0);
        ix1r = realloc(ix1r, nvr*(d-1)*sizeof(ppv_index_t));
        (*ix1rP) = ix1r;
        (*nvrP) = nvr;
       }

    /* Its position in {T} will be

        {T.base + ((ix[0]+dx[0])%sz0T)*T.step[0] + (ix[1]+dx[1])*T.step[1] + ...(ix[d-1]+dx[d-1])*T.step[d-1] }
        
      which is the position in {T} of the voxel with indices 
      
        {{(ix[0]+dx[0])%sz0T, ix[1],.. ix[d-1]}}
        
      plus the displacement {dx[1]*T.step[1] + ... + dx[d-1]*T.step[d-1]}. 
      
      This displacement is obtained as {b->step[iabs(dx[0])][iv]} where {iv} ranges
      in {0..b->nv[iabs(dx[0])]-1}.  The corresponding voxel index increment 
      pairs {(dx,dx[1])} for each of those pixels are stored in {indb[iabs(dx[0])][iv]}. 
    */

    for (ppv_index_t r0 = 0; r0 <= sz0b; r0++)
      { ppv_index_t ix0b = r0 + b->ix0min;  /* Brush slice index in {b} (rel center) */
        ppv_index_t ix0A = ixA[0] + ix0b;   /* Brush slice index in {A}. */
        if (debug) { fprintf(stderr, "  checking slice %ld of original A\n", ix0A); }
        if ((ix0A >= 0) && (ix0A < sz0A))
          { /* Scan ball neighborhood elements of original {A} in slice {ix0A}. */
            ppv_sample_count_t nvz = b->nv[r0];
            ppv_index_t *ix1r = b->ix1[r0];
            ppv_step_t *stepr = b->step[r0];
            /* Compute address of voxel {{ix0,ixA[1..d-1]}} on this slice of {T}. */ 
            ixT[0] = ix0 % sz0T;
            if (debug) { fprintf(stderr, "  assuming it is in slice %ld of T\n", ixT[0]); }
            for(ppv_axis_t ax = 1; ax < d; ax++) { ixT[ax] = ixA[ax]; }
            ppv_pos_t baser = ppv_sample_pos(T, ixT); 
            for (ppv_sample_count_t k = 0; k < nvr; k++)
              { /* Coompute the indices {1..d-1] of this voxel in {T}: */
                int32_t k 
                ppv_index_t rx = kx + indk.c[0]; 
                ppv_index_t ry = ky + indk.c[1]; 
                if ((rx >= 0) && (rx < nxA) && (ry >= 0) && (ry < nyA))
                  { /* Get the ball voxel from {T}: */
                    ppv_pos_t pos = baser + stepr[k];
                    ppv_sample_t vsmp = ppv_get_sample_at_pos(T->el,T->bps,T->bpw, pos);
                    if (debug) { fprintf(stderr, "    T[%ld,%ld,%ld] = %d\n", ix0,ry,rx, vsmp); }
                    if (vsmp == vexp) { return vexp; }
                  }
              }
          }
      }
    return 1-vexp;

/* Pseudo-pseudorandom sample value generator: */

    ppv_sample_count_t mu = 37; /* Multiplicative coeff. */
    ppv_sample_count_t so = 417; /* Additive coeff. */
    ppv_sample_count_t maxstate = (ppv_MAX_SAMPLES - so)/mu;
    ppv_sample_count_t state = (ppv_sample_count_t)((seed | 1) & maxval);
    if (bps > 0) { while (state < (maxstate >> bps)) { state = (state << bps) | seed; } }
    
    for(...)
      { 
        ppv_sample_t val = (ppv_sample_t)(state & maxval);
        ...
        ppv_sample_count_t temp = mu*state + so;
        state = (temp % maxstate);
      }

    ppv_sample_count_t nv = 1;
    for (ppv_axis_t ax = 0; ax < A->d; ax++)
      { ppv_step_t stepk = A->step[ax];
        ppv_size_t sizek = A->size[ax];
        if ((sizek >= 2) && (stepk == 0) && (! reptoo)) { sizek = 1; }
        if (sizek == 0) 
          { nv = 0; break; }
        else
          { assert(nv <= ppv_MAX_SIZE/sizek);
            nv = nv*sizek;
          }
      }