Pawel Netzel - software
Repositories
Help
Report an Issue
plGeoAdaptels
Code
Commits
Branches
Tags
Search
Tree:
e8ae190
Branches
Tags
master
plGeoAdaptels
generate.c
new distances
netzel
commited
e8ae190
at 2024-03-25 09:42:52
generate.c
Blame
History
Raw
#include <stdlib.h> #include <math.h> #include <stdio.h> #include "minlist.h" #include "seeds.h" #include "labels.h" int* create_shiftIdx(int cols) { int* t = (int*)calloc(8, sizeof(int)); t[0] = -1; t[1] = -cols; t[2] = 1; t[3] = cols; t[4] = -1-cols; t[5] = 1-cols; t[6]= 1+cols; t[7] = -1+cols; return t; } void add_cumul_sum(double* cumul_sum, long index, int nodata, double** layers, int n_layers) { int l; if(nodata) return; for(l = 0; l < n_layers; l++) cumul_sum[l] += layers[l][index]; } void init_cumul_sum(double* cumul_sum, long index, int nodata, double** layers, int n_layers) { int l; if(nodata) for(l = 0; l < n_layers; l++) cumul_sum[l] = 0; else for(l = 0; l < n_layers; l++) cumul_sum[l] = layers[l][index]; } double calc_distance(double** layers, int n_layers, double* cumul_sum, long index, long superpixel_size, int distance, double minkowski_p) { double diff, dist = 0.0, A = 1.0, B = 1.0, AB = 1.0; int i; switch(distance) { case 0: // minkowski for(i=0; i<n_layers; i++) { diff = fabs(cumul_sum[i] - layers[i][index]*superpixel_size); dist += pow(diff, minkowski_p); } dist = pow(dist, 1.0/minkowski_p)/superpixel_size; break; case 1: // cosine for(i=0; i<n_layers; i++) { diff = cumul_sum[i]/superpixel_size; A += pow(diff, 2.0); B += pow(layers[i][index],2.0); AB += layers[i][index]*diff; } dist = AB/sqrt(A*B); break; case 2: // angular for(i=0; i<n_layers; i++) { diff = cumul_sum[i]/superpixel_size; A += pow(diff, 2.0); B += pow(layers[i][index],2.0); AB += layers[i][index]*diff; } dist = acos(AB/sqrt(A*B))/M_PI; break; default: dist = 0.0; } return dist; } void get_cell_neighbour(CELL* c, int n, int* x, int* y) { const int dx[] = {-1, 0, 1, 0, -1, 1, 1, -1}; const int dy[] = { 0, -1, 0, 1, -1, -1, 1, 1}; *x = c->x + dx[n]; *y = c->y + dy[n]; } void grow_adaptel(int i, LABELS* labels, MINLIST* minlist, SEEDS* seeds, double* distances, double** layers, int n_layers, unsigned char* mask, int cols, int rows, long size, double threshold, int connectivity, int* dIdx, int distance, double minkowski_p) { long idx1, idx2, index, superpixel_size; double* cumul_sum = (double*)malloc(sizeof(double)*n_layers); double dist; int conn, nodata, neighbour_x, neighbour_y; CELL cell; // get seed index = get_seed_index(seeds, i); if(labels->labels[index] < 0) { insert_minlist(minlist, get_seed_x(seeds, i), get_seed_y(seeds, i), index, 0); distances[index] = 0.0; init_label(labels, index); nodata = (mask[index]==1); superpixel_size = 1; init_cumul_sum(cumul_sum, index, nodata, layers, n_layers); while(minlist->N) { getmin_minlist(minlist, &cell); idx1 = cell.index; for(conn = 0; conn < connectivity; conn++) { get_cell_neighbour(&cell, conn, &neighbour_x, &neighbour_y); if(!(neighbour_x < 0 || neighbour_x >= cols || neighbour_y < 0 || neighbour_y >= rows)) { idx2 = cell.index + dIdx[conn]; if(distances[idx1] < threshold || mask[idx1]==1) { if(labels->labels[idx2] != labels->labels[idx1]) { nodata = (mask[idx2]==1); if(nodata) dist = 0.0; else dist = distances[idx1] + calc_distance(layers, n_layers, cumul_sum, idx2, superpixel_size, distance, minkowski_p); if(dist < distances[idx2] || labels->labels[idx2] < 0) { distances[idx2] = dist; labels->labels[idx2] = labels->labels[idx1]; add_cumul_sum(cumul_sum, idx2, nodata, layers, n_layers); if(nodata) insert_minlist(minlist, neighbour_x, neighbour_y, idx2, 0.0); else insert_minlist(minlist, neighbour_x, neighbour_y, idx2, distances[idx2]); superpixel_size++; } } } else if(labels->labels[idx2] < 0) // insert new seed insert_seed(seeds, i, neighbour_x, neighbour_y, idx2); } } } labels->max_label++; } } void remove_nodata_labels(LABELS* labels, unsigned char* mask, long size) { int* labels_list = (int*)calloc(labels->max_label, sizeof(int)); long i; int j, lab; int id = 0; for(i=0; i<size; i++) if(mask[i]==1) labels->labels[i] = -9999; for(j=0; j<labels->max_label; j++) labels_list[j] = -1; for(i=0; i<size; i++) { if(mask[i]==0) { lab = labels->labels[i]; if(labels_list[lab]==-1) labels_list[lab] = id++; } } for(i=0; i<size; i++) if(mask[i]==0) labels->labels[i] = labels_list[labels->labels[i]]; labels->max_label = id; free(labels_list); } LABELS* create_adaptels( double** layers, const int n_layers, unsigned char* mask, int cols, int rows, const double threshold, const int connectivity, const int distance, const double minkowski_p) { const long size = (long)cols*rows; long i; LABELS* labels = create_labels(size); int* dIdx = create_shiftIdx(cols); // output array // layer with distances double* distances = (double*)calloc(sizeof(double), size); MINLIST* minlist = create_minlist(1000000); SEEDS* seeds = create_seeds(size); // insert initial seed insert_seed(seeds, 0, cols/2, rows/2, (long)((rows/2)*cols+cols/2)); // create adaptels for(i = 0; i < seeds->N; i++) grow_adaptel(i, labels, minlist, seeds, distances, layers, n_layers, mask, cols, rows, size, threshold, connectivity, dIdx, distance, minkowski_p); free_minlist(minlist); free_seeds(seeds); free(dIdx); free(distances); // remove labels for nodata adaptels remove_nodata_labels(labels, mask, size); return labels; }