1
2 /****************************************************************************
3 *
4 * MODULE: i.maxlik
5 * AUTHOR(S): Michael Shapiro (USACERL) & Tao Wen (UIUC)
6 * (original contributors)
7 * Markus Neteler <neteler itc.it>,
8 * Roberto Flor <flor itc.it>,
9 * Bernhard Reiter <bernhard intevation.de>,
10 * Brad Douglas <rez touchofmadness.com>,
11 * Glynn Clements <glynn gclements.plus.com>,
12 * Jan-Oliver Wagner <jan intevation.de>
13 * PURPOSE: maximum likelihood classification of image groups
14 * COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
15 *
16 * This program is free software under the GNU General Public
17 * License (>=v2). Read the file COPYING that comes with GRASS
18 * for details.
19 *
20 *****************************************************************************/
21
22 #include <stdlib.h>
23 #include <grass/gis.h>
24 #include <grass/raster.h>
25 #include <grass/glocale.h>
26 #include "global.h"
27 #include "local_proto.h"
28
29 char *group;
30 char *subgroup;
31 char *sigfile;
32 struct Ref Ref;
33 struct Signature S;
34 DCELL **cell;
35 int *cellfd;
36 CELL *class_cell, *reject_cell;
37 int class_fd, reject_fd;
38 char *reject_name;
39 char class_name[GNAME_MAX];
40 double *B;
41 double *P;
42 CELL cat;
43
main(int argc,char * argv[])44 int main(int argc, char *argv[])
45 {
46 struct Categories cats;
47 struct Colors colr;
48 struct Ref group_ref;
49 int nrows, ncols;
50 int row;
51 int band;
52 int i;
53 struct GModule *module;
54 struct
55 {
56 struct Option *group, *subgroup, *sigfile, *class, *reject;
57 } parm;
58 char xmapset[GMAPSET_MAX];
59
60 G_gisinit(argv[0]);
61
62 module = G_define_module();
63 G_add_keyword(_("imagery"));
64 G_add_keyword(_("classification"));
65 G_add_keyword(_("Maximum Likelihood Classification"));
66 G_add_keyword("MLC");
67 module->label =
68 _("Classifies the cell spectral reflectances in imagery data.");
69 module->description =
70 _("Classification is based on the spectral signature information "
71 "generated by either i.cluster, g.gui.iclass, or i.gensig.");
72
73 parm.group = G_define_standard_option(G_OPT_I_GROUP);
74
75 parm.subgroup = G_define_standard_option(G_OPT_I_SUBGROUP);
76
77 parm.sigfile = G_define_option();
78 parm.sigfile->key = "signaturefile";
79 parm.sigfile->required = YES;
80 parm.sigfile->type = TYPE_STRING;
81 parm.sigfile->gisprompt = "old,sig,sigfile";
82 parm.sigfile->key_desc = "name";
83 parm.sigfile->label = _("Name of input file containing signatures");
84 parm.sigfile->description = _("Generated by either i.cluster, g.gui.iclass, or i.gensig");
85
86 parm.class = G_define_standard_option(G_OPT_R_OUTPUT);
87 parm.class->required = YES;
88 parm.class->description = _("Name for output raster map holding classification results");
89
90 parm.reject = G_define_standard_option(G_OPT_R_OUTPUT);
91 parm.reject->key = "reject";
92 parm.reject->required = NO;
93 parm.reject->description =
94 _("Name for output raster map holding reject threshold results");
95
96 if (G_parser(argc, argv))
97 exit(EXIT_FAILURE);
98
99
100 reject_name = parm.reject->answer;
101 group = parm.group->answer;
102 subgroup = parm.subgroup->answer;
103 sigfile = parm.sigfile->answer;
104
105 if (G_unqualified_name(parm.class->answer, G_mapset(), class_name, xmapset) < 0)
106 G_fatal_error(_("<%s> does not match the current mapset"), xmapset);
107
108 if (G_legal_filename(class_name) < 0)
109 G_fatal_error(_("<%s> is an illegal file name"), class_name);
110
111 open_files();
112
113 nrows = Rast_window_rows();
114 ncols = Rast_window_cols();
115
116 for (row = 0; row < nrows; row++) {
117 G_percent(row, nrows, 2);
118
119 for (band = 0; band < Ref.nfiles; band++)
120 Rast_get_d_row(cellfd[band], cell[band], row);
121
122 classify(class_cell, reject_cell, ncols);
123 Rast_put_row(class_fd, class_cell, CELL_TYPE);
124 if (reject_fd > 0)
125 Rast_put_row(reject_fd, reject_cell, CELL_TYPE);
126 }
127 G_percent(nrows, nrows, 2);
128
129 Rast_close(class_fd);
130 if (reject_fd > 0)
131 Rast_close(reject_fd);
132
133 Rast_init_cats("Maximum Likelihood Classification", &cats);
134 for (i = 0; i < S.nsigs; i++) {
135 if (*S.sig[i].desc) {
136 cat = i + 1;
137 Rast_set_c_cat(&cat, &cat, S.sig[i].desc, &cats);
138 }
139 }
140 Rast_write_cats(class_name, &cats);
141 Rast_free_cats(&cats);
142
143 if (reject_fd > 0) {
144 char title[100];
145 char *label[] = { "no data", "0.1%", "0.5%",
146 "1%", "2%", "5%", "10%",
147 "20%", "30%", "50%", "70%",
148 "80%", "90%", "95%", "98%",
149 "99%", "100%", "bad" };
150 sprintf(title, "Rejection Probability for %s", class_name);
151
152 Rast_init_cats(title, &cats);
153 Rast_set_cats_title(title, &cats);
154 for(i = 0; i < (int) (sizeof(label) / sizeof (char *)); i++) {
155 Rast_set_c_cat(&i, &i, label[i], &cats);
156 }
157 Rast_write_cats(reject_name, &cats);
158 Rast_free_cats(&cats);
159
160 Rast_make_grey_scale_colors(&colr, (CELL) 1, (CELL) 16);
161
162 Rast_set_c_color((CELL) 0, 0, 255, 0, &colr);
163 Rast_set_c_color((CELL) 17, 255, 0, 0, &colr);
164 Rast_write_colors(reject_name, G_mapset(), &colr);
165 Rast_free_colors(&colr);
166 }
167
168 /* associate the output files with the group */
169 I_get_group_ref(group, &group_ref);
170 I_add_file_to_group_ref(class_name, G_mapset(), &group_ref);
171 if (reject_cell)
172 I_add_file_to_group_ref(reject_name, G_mapset(), &group_ref);
173
174 I_put_group_ref(group, &group_ref);
175 make_history(class_name, group, subgroup, sigfile);
176
177 exit(EXIT_SUCCESS);
178 }
179