Pawel Netzel - software
Repositories
Help
Report an Issue
plClump4p
Code
Commits
Branches
Tags
Search
Tree:
38f25a5
Branches
Tags
master
plClump4p
clump.c
Initial commit
netzel
commited
38f25a5
at 2024-02-21 22:59:53
clump.c
Blame
History
Raw
/**************************************************************************** * * PROGRAM: clump4p * * AUTHOR(S): Pawel Netzel - SIL UC * Jarek Jasiewicz - SIL UC * (the clumping algorithm for a region) * * PURPOSE: Recategorizes data in a raster map layer by grouping cells * that form physically discrete areas into unique categories. * * COPYRIGHT: (C) 2012 by Pawel Netzel * * This program is a free software under the GNU General Public * License (>=v2). * ***************************************************************************/ #include <stdio.h> #include <string.h> #include <time.h> #include <gdal.h> #include <ezgdal.h> #include <omp.h> #include "main.h" #include "data.h" static int nextr[] = {-1, 0, 1, 0 ,-1, 1,-1, 1}; static int nextc[] = { 0,-1, 0, 1 ,-1, 1, 1,-1}; #define NR(x) 1 + r + nextr[(x)] #define NC(x) 1 + c + nextc[(x)] #define INDEX(r,c) ((r)-1)*ncols+(c)-1 void clumpRegion(int regionNum, int nDirections) { int** map; int** map_clump; int nrows; int ncols; int colMax,colMin; unsigned current_clump, index, i; unsigned row, col, r, c; int last, first, next_r, next_c; //double sum=0; unsigned* queue; /* The size of the region */ colMin=regions[regionNum].colMin; colMax=regions[regionNum].colMax; ncols=colMax-colMin+1; nrows=regions[regionNum].rowMax-regions[regionNum].rowMin+1; /* Allocating memory and copying input data */ queue=(unsigned*)malloc((nrows+2)*(ncols+2)*sizeof(unsigned)); map=malloc((nrows+2)*sizeof(int*)); map_clump=malloc((nrows+2)*sizeof(int*)); map_clump[0]=calloc(ncols+2,sizeof(int)); map_clump[nrows+1]=calloc(ncols+2,sizeof(int)); map[0]=calloc(ncols+2,sizeof(int)); map[nrows+1]=calloc(ncols+2,sizeof(int)); for(r=1; r<=nrows; r++) { map_clump[r]=calloc(ncols+2,sizeof(int)); map[r]=calloc(ncols+2,sizeof(int)); // Data copying memcpy(map[r]+1,dataBuffer[r-1]+colMin,ncols*sizeof(int)); } /* start clumping */ current_clump=0; for(row=1;row<=nrows;row++) for(col=1;col<=ncols;col++) { if(map_clump[row][col] || !map[row][col]) /* NULL == 0 */ continue; current_clump++; map_clump[row][col]=current_clump; last=1; first=0; queue[first]=INDEX(row,col); do { index=queue[first++]; r=index/ncols; c=index%ncols; for(i=0; i<nDirections; i++) { next_r=NR(i); next_c=NC(i); if(map[r+1][c+1]==map[next_r][next_c]) { if(!map_clump[next_r][next_c]) { map_clump[next_r][next_c]=current_clump; queue[last++]=INDEX(next_r,next_c); } } } } while(first!=last); } /* Creating labels list */ if(current_clump>INCR) { free(regions[regionNum].labels); regions[regionNum].labels=(int *)malloc((current_clump+2)*sizeof(int)); } for(i=1; i<=current_clump; i++) regions[regionNum].labels[i]=i+1; regions[regionNum].nLabels=current_clump; /* Copying resultant map of clumps */ for(r=1;r<=nrows;r++) memcpy(resultBuffer[r-1]+colMin,map_clump[r]+1,ncols*sizeof(int)); /* Releasing memmory */ for(r=0;r<=nrows+1;r++) { free(map[r]); free(map_clump[r]); } free(map); free(map_clump); free(queue); } void shiftLablesForRegion(int regionNum) { int i,nLabels; int lbl,maxLbl, *labels; labels=regions[regionNum].labels; nLabels=regions[regionNum].nLabels; maxLbl=0; for(i=1; i<=nLabels; i++) { lbl=labels[i]; labels[i]+=totalIndex; if(lbl>maxLbl) maxLbl=lbl; } totalIndex+=maxLbl; } void resolveEquivalence(int labelP, int regionP, int labelQ, int regionQ) { int idxP,idxQ,maxLab, minLab,n;//,*kLabels; int j,k; Region* reg; int* lab; idxP=regions[regionP].labels[labelP]; idxQ=regions[regionQ].labels[labelQ]; if(idxP!=idxQ) { maxLab=(idxP>idxQ)?idxP:idxQ; minLab=(idxP>idxQ)?idxQ:idxP; #pragma omp parallel for private(k,j,n,lab,reg) shared(maxLab,minLab,regions,regionP) for(k=0; k<=regionP; k++) { reg=regions+k; n=(*reg).nLabels; lab=(*reg).labels+1; if(k<regionP) reg++; for(j=1; j<=n; j++) { if(*lab>maxLab) (*lab)--; else if(*lab==maxLab) *lab=minLab; lab++; } } totalIndex--; } } void merge4(int regionNum) { int colMin,colMax,i,j,n; colMin=regions[regionNum].colMin; colMax=regions[regionNum].colMax; /* Left vertical border */ if(colMin>0) { n=regions[regionNum].rowMax-regions[regionNum].rowMin+1; j=colMin-1; for(i=0; i<n; i++) if(resultBuffer[i][colMin]>0) if(resultBuffer[i][j]>0 && dataBuffer[i][j]==dataBuffer[i][colMin]) resolveEquivalence(resultBuffer[i][colMin],regionNum, resultBuffer[i][j],regionNum-1); } /* Top horizontal border */ if(regions[regionNum].rowMin>0) for(i=colMin; i<=colMax; i++) { if(resultBuffer[0][i]>0) if(borderResultBuffer[i]>0 && dataBuffer[0][i]==borderDataBuffer[i]) resolveEquivalence(resultBuffer[0][i],regionNum, borderResultBuffer[i], regions[regionNum].regionUp); } } void merge8(int regionNum) { int colMin,colMax,i,j,n; colMin=regions[regionNum].colMin; colMax=regions[regionNum].colMax; /* Left vertical border */ if(colMin>0) { n=regions[regionNum].rowMax-regions[regionNum].rowMin; j=colMin-1; for(i=1; i<n; i++) if(resultBuffer[i][colMin]>0) { if(resultBuffer[i-1][j]>0 && dataBuffer[i-1][j]==dataBuffer[i][colMin]) resolveEquivalence(resultBuffer[i][colMin],regionNum, resultBuffer[i-1][j],regionNum-1); if(resultBuffer[i][j]>0 && dataBuffer[i][j]==dataBuffer[i][colMin]) resolveEquivalence(resultBuffer[i][colMin],regionNum, resultBuffer[i][j],regionNum-1); if(resultBuffer[i+1][j]>0 && dataBuffer[i+1][j]==dataBuffer[i][colMin]) resolveEquivalence(resultBuffer[i][colMin],regionNum, resultBuffer[i+1][j],regionNum-1); } if(resultBuffer[n][colMin]>0) { if(resultBuffer[n-1][j]>0 && dataBuffer[n-1][j]==dataBuffer[n][colMin]) resolveEquivalence(resultBuffer[n][colMin],regionNum, resultBuffer[n-1][j],regionNum-1); if(resultBuffer[n][j]>0 && dataBuffer[n][j]==dataBuffer[n][colMin]) resolveEquivalence(resultBuffer[n][colMin],regionNum, resultBuffer[n][j],regionNum-1); } } /* Left, top corner */ if(resultBuffer[0][colMin]>0) { /* Left region*/ if(colMin>0) { if(resultBuffer[0][colMin-1]>0 && dataBuffer[0][colMin]==dataBuffer[0][colMin-1]) resolveEquivalence(resultBuffer[0][colMin],regionNum, resultBuffer[0][colMin-1],regionNum-1); if(resultBuffer[1][colMin-1]>0 && dataBuffer[0][colMin]==dataBuffer[1][colMin-1]) resolveEquivalence(resultBuffer[0][colMin],regionNum, resultBuffer[1][colMin-1],regionNum-1); } /* Left-Top region*/ if((regions[regionNum].rowMin>0)&&(colMin>0)) if(borderResultBuffer[colMin-1]>0 && dataBuffer[0][colMin]==borderDataBuffer[colMin-1]) resolveEquivalence(resultBuffer[0][colMin],regionNum, borderResultBuffer[colMin-1], regions[regionNum].regionUp-1); /* Top region */ if(regions[regionNum].rowMin>0) { if(borderResultBuffer[colMin]>0 && dataBuffer[0][colMin]==borderDataBuffer[colMin]) resolveEquivalence(resultBuffer[0][colMin],regionNum, borderResultBuffer[colMin], regions[regionNum].regionUp); if(borderResultBuffer[colMin+1]>0 && dataBuffer[0][colMin]==borderDataBuffer[colMin+1]) resolveEquivalence(resultBuffer[0][colMin],regionNum, borderResultBuffer[colMin+1], regions[regionNum].regionUp); } } /* Top horizontal border */ if(regions[regionNum].rowMin>0) { for(i=colMin+1; i<colMax; i++) { if(resultBuffer[0][i]>0) { if(borderResultBuffer[i-1]>0 && dataBuffer[0][i]==borderDataBuffer[i-1]) resolveEquivalence(resultBuffer[0][i],regionNum, borderResultBuffer[i-1], regions[regionNum].regionUp);\ if(borderResultBuffer[i]>0 && dataBuffer[0][i]==borderDataBuffer[i]) resolveEquivalence(resultBuffer[0][i],regionNum, borderResultBuffer[i], regions[regionNum].regionUp);\ if(borderResultBuffer[i+1]>0 && dataBuffer[0][i]==borderDataBuffer[i+1]) resolveEquivalence(resultBuffer[0][i],regionNum, borderResultBuffer[i+1], regions[regionNum].regionUp);\ } } if(resultBuffer[0][i]>0) { if(borderResultBuffer[colMax-1]>0 && dataBuffer[0][colMax]==borderDataBuffer[colMax-1]) resolveEquivalence(resultBuffer[0][colMax],regionNum, borderResultBuffer[colMax-1], regions[regionNum].regionUp); if(borderResultBuffer[colMax]>0 && dataBuffer[0][colMax]==borderDataBuffer[colMax]) resolveEquivalence(resultBuffer[0][colMax],regionNum, borderResultBuffer[colMax-1], regions[regionNum].regionUp); } } } void merge(int n, int nDirections) { if(nDirections==4) merge4(n); else merge8(n); } void clump(EZGDAL_LAYER *in_fd, EZGDAL_LAYER *out_fd, int nDirections) { int rowMin,rowMax,i,j,n,reg; totalIndex=0; for(i=0; i<nRegionRows; i++) { reg=i*nRegionCols; ezgdal_show_progress(stdout,reg,nRegions); n=nDataCols*sizeof(int); if(i==0) { memset(borderDataBuffer,0,n); memset(borderResultBuffer,0,n); } else { memcpy(borderDataBuffer,dataBuffer[rowMax-rowMin],n); memcpy(borderResultBuffer,resultBuffer[rowMax-rowMin],n); } rowMin=regions[reg].rowMin; rowMax=regions[reg].rowMax; loadDataBuffer(in_fd,rowMin,rowMax); #pragma omp parallel for private(j) shared(reg,nRegionCols,regions, dataBuffer, resultBuffer, borderDataBuffer,borderResultBuffer) for(j=0; j<nRegionCols; j++) clumpRegion(reg+j, nDirections); for(j=0; j<nRegionCols; j++) { n=reg+j; shiftLablesForRegion(n); merge(n,nDirections); } saveResultBuffer(out_fd,rowMin,rowMax); } ezgdal_show_progress(stdout,nRegions,nRegions); } void labelsOrdering() { int *index,lbl; int i,j,x; index=(int *)calloc(totalIndex+1,sizeof(int)); x=1; for(i=0; i<nRegions; i++) for(j=1; j<=regions[i].nLabels; j++) { lbl=regions[i].labels[j]; if(index[lbl]==0) index[lbl]=x++; } #pragma omp parallel for private(i,j) for(i=0; i<nRegions; i++) for(j=1; j<=regions[i].nLabels; j++) regions[i].labels[j]=index[regions[i].labels[j]]; free(index); } void replaceLabels(EZGDAL_LAYER *in_fd, EZGDAL_LAYER *out_fd) { int rowMin,rowMax,i,colMin,colMax; int row,rows,col,cols,colAbs; ezgdal_show_progress(stdout,0,nRegions); labelsOrdering(); rowMin=regions[0].rowMin; rowMax=regions[0].rowMax; loadDataBuffer(in_fd,rowMin,rowMax); for(i=0; i<nRegions; i++) { ezgdal_show_progress(stdout,i,nRegions); if(rowMin<regions[i].rowMin) { saveResultBuffer(out_fd,rowMin,rowMax); rowMin=regions[i].rowMin; rowMax=regions[i].rowMax; loadDataBuffer(in_fd,rowMin,rowMax); } colMin=regions[i].colMin; colMax=regions[i].colMax; cols=colMax-colMin+1; rows=rowMax-rowMin+1; #pragma omp parallel for private(row,col,colAbs) shared(rows,cols,dataBuffer,resultBuffer,regions) for(row=0; row<rows; row++) for(col=0; col<cols; col++) { colAbs=col+colMin; if(!ezgdal_is_null(in_fd,(double)dataBuffer[row][colAbs])) resultBuffer[row][colAbs]=regions[i].labels[dataBuffer[row][colAbs]]; else resultBuffer[row][colAbs]=0; } } saveResultBuffer(out_fd,rowMin,rowMax); ezgdal_show_progress(stdout,i,nRegions); }