1 /* pgmminkowsky.c - read a portable graymap and calculate the Minkowski
2 ** Integrals as a function of the threshold.
3 **
4 ** Copyright (C) 2000 by Luuk van Dijk/Mind over Matter
5 **
6 ** Based on pgmhist.c,
7 ** Copyright (C) 1989 by Jef Poskanzer.
8 **
9 ** Permission to use, copy, modify, and distribute this software and its
10 ** documentation for any purpose and without fee is hereby granted, provided
11 ** that the above copyright notice appear in all copies and that both that
12 ** copyright notice and this permission notice appear in supporting
13 ** documentation. This software is provided "as is" without express or
14 ** implied warranty.
15 */
16
17 #include "pgm.h"
18 #include "mallocvar.h"
19
20
21 #define MAX2(a,b) ( ( (a)>(b) ) ? (a) : (b) )
22 #define MAX4(a,b,c,d) MAX2( MAX2((a),(b)), MAX2((c),(d)) )
23
main(int argc,char ** argv)24 int main( int argc, char** argv ){
25
26 FILE *ifp;
27
28 gray maxval;
29 int cols, rows, format;
30
31 gray* prevrow;
32 gray* thisrow;
33 gray* tmprow;
34
35 int* countTile;
36 int* countEdgeX;
37 int* countEdgeY;
38 int* countVertex;
39
40 int i, col, row;
41
42 int maxtiles, maxedgex, maxedgey, maxvertex;
43 int area, perimeter, eulerchi;
44
45 double l2inv, linv;
46
47 /*
48 * parse arg and initialize
49 */
50
51 pgm_init( &argc, argv );
52
53 if ( argc > 2 ) pm_usage( "[pgmfile]" );
54
55 if ( argc == 2 )
56 ifp = pm_openr( argv[1] );
57 else
58 ifp = stdin;
59
60 /*
61 * initialize
62 */
63
64 pgm_readpgminit( ifp, &cols, &rows, &maxval, &format );
65
66 prevrow = pgm_allocrow( cols );
67 thisrow = pgm_allocrow( cols );
68
69 MALLOCARRAY(countTile , maxval + 1 );
70 MALLOCARRAY(countEdgeX , maxval + 1 );
71 MALLOCARRAY(countEdgeY , maxval + 1 );
72 MALLOCARRAY(countVertex , maxval + 1 );
73
74 if (countTile == NULL || countEdgeX == NULL || countEdgeY == NULL ||
75 countVertex == NULL)
76 pm_error( "out of memory" );
77
78 for ( i = 0; i <= maxval; i++ ) countTile[i] = 0;
79 for ( i = 0; i <= maxval; i++ ) countEdgeX[i] = 0;
80 for ( i = 0; i <= maxval; i++ ) countEdgeY[i] = 0;
81 for ( i = 0; i <= maxval; i++ ) countVertex[i] = 0;
82
83
84
85
86 /* first row */
87
88 pgm_readpgmrow( ifp, thisrow, cols, maxval, format );
89
90 /* tiles */
91
92 for ( col = 0; col < cols; ++col ) ++countTile[thisrow[col]];
93
94 /* y-edges */
95
96 for ( col = 0; col < cols; ++col ) ++countEdgeY[thisrow[col]];
97
98 /* x-edges */
99
100 ++countEdgeX[thisrow[0]];
101
102 for ( col = 0; col < cols-1; ++col )
103 ++countEdgeX[ MAX2(thisrow[col], thisrow[col+1]) ];
104
105 ++countEdgeX[thisrow[cols-1]];
106
107 /* shortcut: for the first row, countVertex == countEdgeX */
108
109 ++countVertex[thisrow[0]];
110
111 for ( col = 0; col < cols-1; ++col )
112 ++countVertex[ MAX2(thisrow[col], thisrow[col+1]) ];
113
114 ++countVertex[thisrow[cols-1]];
115
116
117
118 for ( row = 1; row < rows; ++row ){
119
120 tmprow = prevrow;
121 prevrow = thisrow;
122 thisrow = tmprow;
123
124 pgm_readpgmrow( ifp, thisrow, cols, maxval, format );
125
126 /* tiles */
127
128 for ( col = 0; col < cols; ++col ) ++countTile[thisrow[col]];
129
130 /* y-edges */
131
132 for ( col = 0; col < cols; ++col )
133 ++countEdgeY[ MAX2(thisrow[col], prevrow[col]) ];
134 /* x-edges */
135
136 ++countEdgeX[thisrow[0]];
137
138 for ( col = 0; col < cols-1; ++col )
139 ++countEdgeX[ MAX2(thisrow[col], thisrow[col+1]) ];
140
141 ++countEdgeX[thisrow[cols-1]];
142
143 /* vertices */
144
145 ++countVertex[ MAX2(thisrow[0],prevrow[0]) ];
146
147 for ( col = 0; col < cols-1; ++col )
148 ++countVertex[
149 MAX4(thisrow[col], thisrow[col+1], prevrow[col], prevrow[col+1])
150 ];
151
152 ++countVertex[ MAX2(thisrow[cols-1],prevrow[cols-1]) ];
153
154 } /* for row */
155
156 /* now thisrow contains the top row*/
157
158 /* tiles and x-edges have been counted, now upper
159 y-edges and top vertices remain */
160
161 /* y-edges */
162
163 for ( col = 0; col < cols; ++col ) ++countEdgeY[ thisrow[col] ];
164
165 /* vertices */
166
167 ++countVertex[thisrow[0]];
168
169 for ( col = 0; col < cols-1; ++col )
170 ++countVertex[ MAX2(thisrow[col],thisrow[col+1]) ];
171
172 ++countVertex[ thisrow[cols-1] ];
173
174
175 /* cleanup */
176
177 maxtiles = rows * cols;
178 maxedgex = rows * (cols+1);
179 maxedgey = (rows+1) * cols;
180 maxvertex= (rows+1) * (cols+1);
181
182 l2inv = 1.0/maxtiles;
183 linv = 0.5/(rows+cols);
184
185 /* And print it. */
186 printf( "#threshold\t tiles\tx-edges\ty-edges\tvertices\n" );
187 printf( "#---------\t -----\t-------\t-------\t--------\n" );
188 for ( i = 0; i <= maxval; i++ ){
189
190 if( !(countTile[i] || countEdgeX[i] || countEdgeY[i] || countVertex[i] ) )
191 continue; /* skip empty slots */
192
193 area = maxtiles;
194 perimeter = 2*maxedgex + 2*maxedgey - 4*maxtiles;
195 eulerchi = maxtiles - maxedgex - maxedgey + maxvertex;
196
197 printf( "%f\t%6d\t%7d\t%7d\t%8d\t%g\t%g\t%6d\n", (float) i/(1.0*maxval),
198 maxtiles, maxedgex, maxedgey, maxvertex,
199 area*l2inv, perimeter*linv, eulerchi
200 );
201
202
203 maxtiles -= countTile[i];
204 maxedgex -= countEdgeX[i];
205 maxedgey -= countEdgeY[i];
206 maxvertex-= countVertex[i];
207
208 /* i, countTile[i], countEdgeX[i], countEdgeY[i], countVertex[i] */
209
210 }
211
212 /* these should be zero: */
213 printf( "# check:\t%6d\t%7d\t%7d\t%8d\n",
214 maxtiles, maxedgex, maxedgey, maxvertex );
215
216 pm_close( ifp );
217
218 exit( 0 );
219
220 } /*main*/
221
222
223