/**************************************************************************** * * MODULE: i.zc * AUTHOR(S): David B. Satnik Central Washington University GIS Laboratory * (original contributor), based on code provided by Bill Hoff * at University of Illinois * Markus Neteler , * Bernhard Reiter , * Brad Douglas , * Glynn Clements , * Jan-Oliver Wagner * PURPOSE: edge detection for imagery using zero crossings method * COPYRIGHT: (C) 1999-2006 by the GRASS Development Team * * This program is free software under the GNU General Public * License (>=v2). Read the file COPYING that comes with GRASS * for details. * *****************************************************************************/ #include #include #include #include #include #include int main(int argc, char *argv[]) { /* Global variable & function declarations */ double Thresh; int NumOrients; int inputfd, zcfd; /* the input and output file descriptors */ struct Cell_head window; CELL *cell_row; float Width; int i, j; /* Loop control variables */ int or, oc; /* Original dimensions of image */ int rows, cols; /* Smallest powers of 2 >= number of rows & columns */ int size; /* the length of one side */ long totsize; /* the Total number of data points */ double *data[2]; /* Data structure containing real & complex values of FFT */ struct GModule *module; struct Option *input_map, *output_map, *width, *threshold, *orientations; G_gisinit(argv[0]); module = G_define_module(); G_add_keyword(_("imagery")); G_add_keyword(_("edges")); module->description = _("Zero-crossing \"edge detection\" raster " "function for image processing."); /* define options */ input_map = G_define_option(); input_map->key = "input"; input_map->type = TYPE_STRING; input_map->required = YES; input_map->multiple = NO; input_map->gisprompt = "old,cell,raster"; input_map->description = _("Name of input raster map"); output_map = G_define_option(); output_map->key = "output"; output_map->type = TYPE_STRING; output_map->required = YES; output_map->multiple = NO; output_map->gisprompt = "new,cell,raster"; output_map->description = _("Zero crossing raster map"); width = G_define_option(); width->key = "width"; width->type = TYPE_INTEGER; width->required = NO; width->multiple = NO; width->description = _("x-y extent of the Gaussian filter"); width->answer = "9"; threshold = G_define_option(); threshold->key = "threshold"; threshold->type = TYPE_DOUBLE; threshold->required = NO; threshold->multiple = NO; threshold->description = _("Sensitivity of Gaussian filter"); threshold->answer = "10"; orientations = G_define_option(); orientations->key = "orientations"; orientations->type = TYPE_INTEGER; orientations->required = NO; orientations->multiple = NO; orientations->description = _("Number of azimuth directions categorized"); orientations->answer = "1"; /* call parser */ if (G_parser(argc, argv)) exit(EXIT_FAILURE); /* open input cell map */ inputfd = Rast_open_old(input_map->answer, ""); sscanf(threshold->answer, "%1lf", &Thresh); if (Thresh <= 0.0) G_fatal_error(_("Threshold less than or equal to zero not allowed")); Thresh /= 100.0; sscanf(width->answer, "%f", &Width); if (Width <= 0.0) G_fatal_error(_("Width less than or equal to zero not allowed")); sscanf(orientations->answer, "%d", &NumOrients); if (NumOrients < 1) G_fatal_error(_("Fewer than 1 orientation classes not allowed")); /* get the current window for later */ G_get_set_window(&window); /* get the rows and columns in the current window */ or = Rast_window_rows(); oc = Rast_window_cols(); rows = G_math_max_pow2((long)or); cols = G_math_max_pow2((long)oc); size = (rows > cols) ? rows : cols; totsize = size * size; G_message(_("Power 2 values : %d rows %d columns"), rows, cols); /* Allocate appropriate memory for the structure containing the real and complex components of the FFT. DATA[0] will contain the real, and DATA[1] the complex component. */ data[0] = (double *)G_malloc(totsize * sizeof(double)); data[1] = (double *)G_malloc(totsize * sizeof(double)); /* Initialize real & complex components to zero */ G_message(_("Initializing data...")); for (i = 0; i < (totsize); i++) { *(data[0] + i) = 0.0; *(data[1] + i) = 0.0; } /* allocate the space for one row of cell map data */ cell_row = Rast_allocate_c_buf(); /* Read in cell map values */ G_message(_("Reading raster map...")); for (i = 0; i < or; i++) { Rast_get_c_row(inputfd, cell_row, i); for (j = 0; j < oc; j++) *(data[0] + (i * size) + j) = (double)cell_row[j]; } /* close input cell map and release the row buffer */ Rast_close(inputfd); G_free(cell_row); /* take the del**2g of image */ del2g(data, size, Width); /* find the zero crossings: Here are several notes - 1) this routine only uses the real values 2) it places the zero crossings in the imaginary array */ G_math_findzc(data[0], size, data[1], Thresh, NumOrients); /* open the output cell maps and allocate cell row buffers */ G_message(_("Writing transformed data to file...")); zcfd = Rast_open_c_new(output_map->answer); cell_row = Rast_allocate_c_buf(); /* Write out result to a new cell map */ for (i = 0; i < or; i++) { for (j = 0; j < oc; j++) { *(cell_row + j) = (CELL) (*(data[1] + i * cols + j)); } Rast_put_row(zcfd, cell_row, CELL_TYPE); } Rast_close(zcfd); G_free(cell_row); /* Release memory resources */ for (i = 0; i < 2; i++) G_free(data[i]); G_done_msg(_("Transform successful")); exit(EXIT_SUCCESS); }