dupl.cpp

Go to the documentation of this file.

#include "dupl.h"

#include "geometry.h"

#include "mathlib.h"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


namespace midaspace {


void DuplicatePoints :: add_point_to_cell (long id, long pos[])
{
  int *ccellp = ccellpoints[pos[0]][pos[1]][pos[2]];
  long *cellp =  cellpoints[pos[0]][pos[1]][pos[2]];
  
  if (ccellp == NULL) {
    ccellp = new int;
    *ccellp = 0;
  }
  allocate_another_val (*ccellp, cellp);
  
  cellp[*ccellp] = id;
  (*ccellp)++;
  
  ccellpoints[pos[0]][pos[1]][pos[2]] = ccellp;
  cellpoints [pos[0]][pos[1]][pos[2]] =  cellp;
}

bool DuplicatePoints :: point_is_in_cell (long id, long pos[], double offset)
{
  if (pos[0]==-1    || pos[1]==-1    || pos[2]==-1   ) return false;
  if (pos[0]==cc[0] || pos[1]==cc[1] || pos[2]==cc[2]) return false;
  
  if ( points[id]->x > (mincoo[0] +  pos[0]*dcoo[0] - offset)  &&  points[id]->x < (mincoo[0] + (pos[0]+1)*dcoo[0] + offset)  &&
       points[id]->y > (mincoo[1] +  pos[1]*dcoo[1] - offset)  &&  points[id]->y < (mincoo[1] + (pos[1]+1)*dcoo[1] + offset)  &&
       points[id]->z > (mincoo[2] +  pos[2]*dcoo[2] - offset)  &&  points[id]->z < (mincoo[2] + (pos[2]+1)*dcoo[2] + offset)     )
    return true;
  
  return false;
}

void DuplicatePoints :: assign_cellpoints (void)
{
  long i,j, k, l, m;
  
  // set boundary box
  for (i=0; i<3; i++) {
    mincoo[i] = (*points[0])[i];
    maxcoo[i] = (*points[0])[i];
  }
  
  for (i=1; i<np; i++)
    for (j=0; j<3; j++)
      if      (mincoo[j] > (*points[i])[j])  mincoo[j] = (*points[i])[j];
      else if (maxcoo[j] < (*points[i])[j])  maxcoo[j] = (*points[i])[j];
  
  for (i=0; i<3; i++) {
    mincoo[i] -= 1.2*Zero_dist;   // az to bude chodit, tak tady dej 1.2*Zero_dist misto Zero
    maxcoo[i] += 1.2*Zero_dist;   // !!!!!!!!!! 
  }
  
  // set count of cells in direction of axes
  for (i=0; i<3; i++) {
    cc[i] = 1 + (maxcoo[i] - mincoo[i]) / averdist;
    dcoo[i] = (maxcoo[i] - mincoo[i]) / cc[i];
  }
  
  // alloc cellpoints
  allocate (ccellpoints, cc[0], cc[1], cc[2]);
  allocate ( cellpoints, cc[0], cc[1], cc[2]);
  fill_all_by (ccellpoints, cc[0], cc[1], cc[2],  (int*)NULL);
  fill_all_by ( cellpoints, cc[0], cc[1], cc[2], (long*)NULL);
  
  // allocate points to cells
  long pos[3],pos2[3]; // position
  for (i=0; i<np; i++) {
    for (j=0; j<3; j++)
      pos[j] = ((*points[i])[j] - mincoo[j]) / dcoo[j] ;
    
    add_point_to_cell (i, pos);
    
    if (! point_is_in_cell (i, pos, -1.1*Zero_dist))
      for (k=-1; k<2; k++) {
    pos2[0] = pos[0] + k;
    for (l=-1; l<2; l++) {
      pos2[1] = pos[1] + l;
      for (m=-1; m<2; m++) {
        pos2[2] = pos[2] + m;
        
        if (k==0 && l==0 && m==0) continue;
        
        if (point_is_in_cell (i, pos2, 1.1*Zero_dist))
          add_point_to_cell (i, pos2);
      }
    }
      }
  }
  
  // print 
  if (PRINT) {
    long aux;
    
    fprintf (stdout,"\n Class DuplicatePoints - control print START \n\n");
    
    fprintf (stdout," mincoo = [%lf][%lf][%lf]\n",mincoo[0],mincoo[1],mincoo[2]);
    fprintf (stdout," maxcoo = [%lf][%lf][%lf]\n",maxcoo[0],maxcoo[1],maxcoo[2]);
    fprintf (stdout," cc     = [%ld][%ld][%ld]\n",cc[0],cc[1],cc[2]);
    for (k=0; k<cc[0]; k++)
      for (l=0; l<cc[1]; l++)
    for (m=0; m<cc[2]; m++) {
      if (ccellpoints[k][l][m] == NULL)  aux = 0;
      else                               aux = ccellpoints[k][l][m][0];
      
      fprintf (stdout," cc[%ld][%ld][%ld] :: ccellps %ld  cellps ", k,l,m, aux);
      
      for (i=0; i<aux; i++)
        fprintf (stdout," %ld", cellpoints[k][l][m][i]+1);
      
      fprintf (stdout,"\n");
      
    }
    
    fprintf (stdout,"\n Class DuplicatePoints - control print END \n");
    
  }
}

void DuplicatePoints :: find_duplicitys (long &nd, long *Cduplicity, long **duplicity)
{
  // finds duplicity nodes in each cell
  long i,j, k, l, m, aux;
  long  cP,  *P, **D;
  double d2, d2max;
  
  d2max = Zero_dist * Zero_dist;
  nd=0;
  for (k=0; k<cc[0]; k++) {
    for (l=0; l<cc[1]; l++) {
      for (m=0; m<cc[2]; m++) {
    // JEDNA BUNKA
    if (ccellpoints[k][l][m] == NULL)  continue; // no point in cell
    
    cP = ccellpoints[k][l][m][0]; // count      of points in one cell
    P  =  cellpoints[k][l][m];    // id numbers of points in one cell
    
    allocate (D, cP, cP);
    for (i=0; i<cP; i++)  D[i][0] = 0;
    
    // pokud se najde duplicita mezi points dvou bodu
    // tak do radku 1. bodu se ulozi lid(local id = cistlo radku) 2. bodu a naopak
    // do D[i][0] se uklada pocet pocet bodu duplicitnich s i-th bodem
    // od D[i][1] dal, se ukladaji lid duplicitnich bodu
    for (i=0; i<cP-1; i++)
      for (j=i+1; j<cP; j++) {
        d2 = points[P[i]]->dist2_to (points[P[j]]);
        if (d2 < d2max) {
          D[i][++D[i][0]] = j;
          D[j][++D[j][0]] = i;
        }
      }
    
    // predpokad:
    // * v bunce muze byt nekolik duplicitnich skupinek
    // * kazda duplicitni skupinka o 2 a vice bodech musi byt 'uzavrena'
    //   tj. vsechny body musi byt duplicitni kazdy s kazdym
    //   predpoklada se totiz ze retez duplicit neprekracuje d2max
    
    // najdou se dupl. skupinky a ulozi se do answer poli
    for (i=0; i<cP-1; i++) {
      if (D[i][0] == 0)  continue; // nalezena dupl. skupina
      
      // z matice D se zrusi radky se stejnou dupl. skupinou
      for (j=i+1; j<cP; j++)
        if (D[j][0] && is_member_of_array(i, D[j][0], D[j]+1)) {
          if (D[i][0] != D[j][0]) _errorr("error");
          D[j][0] = 0;
        }
      
      aux = D[i][0]+1;
      D[i][0] = i;
      
      // hledam, zda dupl.skup. uz nebyla objevena v jine cell
      for (j=0; j<nd; j++)
        if (duplicity[j][0] == P[i]) break;
      
      // byla objevena
      if (j!=nd) {
        // zkontroluju ze jsou opravdu stejne
        if (Cduplicity[j] != aux) _errorr("error");
        aux = j;
        for (j=0; j<Cduplicity[aux]; j++)
          if (duplicity[aux][j] != P[D[i][j]]) _errorr("error");
      }
      // nebyla objevena
      else { 
        // duplicitni skupinka se lokalizuje do answer poli
        Cduplicity[nd] = aux;
        duplicity[nd] = new long[Cduplicity[nd]];
        for (j=0; j<Cduplicity[nd]; j++)
          duplicity[nd][j] = P[D[i][j]];
        
        nd++;
      }
    }
    
    deallocateCheck ( D, cP);
    
    delete    ccellpoints[k][l][m];
    delete []  cellpoints[k][l][m];
    // JEDNA BUNKA
      }
      delete [] ccellpoints[k][l];
      delete []  cellpoints[k][l];
    }
    delete [] ccellpoints[k];
    delete []  cellpoints[k];
  }
  delete [] ccellpoints;
  delete []  cellpoints;
  
  // kontrola, ve vsech dupl skupinkach dohromady by nemel byt zadny stejny point
  long ndp = 0;
  for (k=0; k<nd; k++)  ndp += Cduplicity[k];
  
  long *dp = new long[ndp];
  
  ndp = 0;
  for (k=0; k<nd; k++)
    for (l=0; l<Cduplicity[k]; l++)
      dp[ndp++] = duplicity[k][l];
  
  if (! members_are_unique (ndp, dp)) _errorr ("error is not unique");
  
  delete [] dp;
  
  
  // print
  if (PRINT) {
    fprintf (stdout,"\n DUPLICITY: START \n");
    
    fprintf (stdout,"\n Number of duplicities %ld\n",nd);
    for (i=0; i<nd; i++) {
      fprintf (stdout," duplicity[%ld] :: ", i+1);
      for (j=0; j<Cduplicity[i]; j++)
    fprintf (stdout,", %ld", duplicity[i][j]+1);
      //if ((i+1)%10 == 0)  fprintf (stdout,"\n");
      fprintf (stdout,"\n");
    }
    
    fprintf (stdout,"\n DUPLICITY: END \n");
  }
}

} // namespace midaspace