1 /* ----------------------------------------------------------------------
2 *
3 * Convert a single-image stereogram to a red/cyan anaglyphic image
4 *
5 * By Scott Pakin <scott+pbm@pakin.org>
6 *
7 * ----------------------------------------------------------------------
8 *
9 * Copyright (C) 2009 Scott Pakin <scott+pbm@pakin.org>
10 *
11 * This program is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or (at
14 * your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful, but
17 * WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program. If not, see http://www.gnu.org/licenses/.
23 *
24 * ----------------------------------------------------------------------
25 */
26
27 #include <stdio.h>
28 #include <string.h>
29 #include <math.h>
30
31 #include "mallocvar.h"
32 #include "nstring.h"
33 #include "shhopt.h"
34 #include "pam.h"
35
36
37 struct cmdlineInfo {
38 /* This structure represents all of the information the user
39 supplied in the command line but in a form easy for the program
40 to use.
41 */
42 int separation;
43 /* Exact separation in pixels between the left and right eye,
44 or -1
45 */
46 int minSeparation;
47 /* Minimum separation in pixels between the left and right eye */
48 gray maxGrayVal;
49 /* Maximum grayscale value to which to scale the image */
50 int swapEyes;
51 /* 0=left red, right cyan; 1=left cyan, right red */
52 const char *inputFilename; /* '-' if stdin */
53 };
54
55
56
57 static void
parseCommandLine(int argc,const char ** const argv,struct cmdlineInfo * const cmdlineP)58 parseCommandLine( int argc, const char ** const argv,
59 struct cmdlineInfo * const cmdlineP ) {
60 /*--------------------------------------------------------------------
61 Parse the command line into a structure.
62 ----------------------------------------------------------------------*/
63
64 optEntry * option_def;
65 /* Instructions to OptParseOptions3 on how to parse our options */
66 optStruct3 opt;
67 unsigned int option_def_index;
68 int maxgrayval;
69
70 maxgrayval = 63; /* default */
71
72 MALLOCARRAY(option_def, 100);
73 option_def_index = 0; /* Incremented by OPTENTRY */
74 MEMSZERO(cmdlineP);
75 cmdlineP->separation = -1;
76
77 OPTENT3('s', "sep", OPT_INT, &cmdlineP->separation, NULL, 0);
78 OPTENT3('g', "gray", OPT_INT, &maxgrayval, NULL, 0);
79 OPTENT3('i', "invert", OPT_FLAG, &cmdlineP->swapEyes, NULL, 0);
80 OPTENT3('m', "minsep", OPT_INT, &cmdlineP->minSeparation, NULL, 0);
81
82 opt.opt_table = option_def;
83 opt.short_allowed = 1;
84 opt.allowNegNum = 0;
85
86 pm_optParseOptions3( &argc, (char **)argv, opt, sizeof(opt), 0 );
87
88 if (argc-1 < 1)
89 cmdlineP->inputFilename = "-";
90 else {
91 cmdlineP->inputFilename = argv[1];
92 if (argc-1 > 1)
93 pm_error("Too many arguments: %u. The only argument is the "
94 "optional input file name", argc-1);
95 }
96 cmdlineP->maxGrayVal = (gray) maxgrayval;
97 }
98
99
100
101 static gray **
readAsGray(const char * const fileName,gray const maxGrayVal,struct pam * const pamP)102 readAsGray( const char * const fileName,
103 gray const maxGrayVal,
104 struct pam * const pamP) {
105 /*--------------------------------------------------------------------
106 Read the input image and convert it to grayscale to reduce the
107 number of "almost but not quite" equal pixels. Return the
108 grayscale array and the initialized PAM structure.
109 ----------------------------------------------------------------------*/
110
111 FILE * fileP;
112 tuple * tuplerow;
113 gray ** grayArray;
114 unsigned int row;
115
116 fileP = pm_openr( fileName );
117
118 pnm_readpaminit( fileP, pamP, PAM_STRUCT_SIZE(tuple_type) );
119
120 tuplerow = pnm_allocpamrow( pamP );
121
122 grayArray = pgm_allocarray( pamP->width, pamP->height );
123
124 for (row = 0; row < pamP->height; ++row) {
125 unsigned int col;
126 pnm_readpamrow( pamP, tuplerow );
127 for (col = 0; col < pamP->width; ++col) {
128 double YP, CbP, CrP;
129
130 pnm_YCbCrtuple( tuplerow[col], &YP, &CbP, &CrP );
131 grayArray[row][col] = (gray)
132 (YP * maxGrayVal / (double)pamP->maxval);
133 }
134 }
135 pnm_freepamrow( tuplerow );
136 pm_close( fileP );
137 return grayArray;
138 }
139
140
141
142 static int
bestEyeSepWeEncountered(int const bestSeparation[3],int const altBestSeparation)143 bestEyeSepWeEncountered(int const bestSeparation[3],
144 int const altBestSeparation) {
145
146 int i;
147
148 for (i = 2; i >= 0; --i) {
149 if (bestSeparation[i] != 0)
150 return bestSeparation[i];
151 }
152 return altBestSeparation;
153 }
154
155
156
157 static int
findRegionEyeSeparation(gray ** const grayArray,int const width,int const height)158 findRegionEyeSeparation( gray ** const grayArray,
159 int const width,
160 int const height ) {
161 /*----------------------------------------------------------------------
162 Determine the number of pixels that corresponds to the separation
163 between the viewer's left eye and right eye. We do this by counting
164 the number of pixels that match N pixels ahead in the image for all
165 N in [1, W/2]. The first big spike in the number of matched pixels
166 determines the N to use for the eye separation. More specifically,
167 if a spike that exceeds 3*stdev+mean is found, the corresponding
168 value of N is taken as the eye separation; otherwise, a spike
169 exceeding 2*stdev+mean is used, then 1*stdev+mean, and finally, the
170 eye separation that produces the minimum average distance between
171 matched pixels. A return value of zero indicates that no eye
172 separation could be determined.
173 ------------------------------------------------------------------------*/
174 int bestSeparation[3];
175 /* Eye separation corresponding to spikes of N+1 standard deviations */
176 int hShift;
177 /* Current horizontal shift */
178 double sumMatches;
179 /* Sum of all matches seen so far */
180 double sumSqMatches;
181 /* Sum of the squares of all matches seen so far */
182 double meanMatches;
183 /* Mean of all matches seen so far */
184 double stdMatches;
185 /* Standard deviation of all matches seen so far */
186 double minAvgDist;
187 /* Min. average distance between matches */
188 int altBestSeparation;
189 /* Shift distance corresponding to the above */
190 unsigned int i;
191
192 /* Try in turn each horizontal shift value from 1 to width/2. A
193 shift of 0 is defined to be a perfect match. A shift of more
194 than width/2 implies that the right-eye image is truncated, which
195 is an unnatural way to construct a crosseyed stereogram.
196 */
197 for (i = 0; i < 3; ++i)
198 bestSeparation[i] = 0;
199
200 altBestSeparation = 0;
201 sumMatches = sumSqMatches = 0.0;
202 meanMatches = stdMatches = minAvgDist = width * height;
203
204 for (hShift = 1; hShift <= width/2; ++hShift) {
205 unsigned int row;
206 unsigned long numMatches; /* Number of matched pixels */
207 double avgDist; /* Average distance between matches */
208
209 numMatches = 0; /* initial value */
210
211 /* Tally the number of matches for this shift distance. */
212 for (row = 0; row < height; ++row) {
213 unsigned int col;
214 for (col = 0; col < width - hShift; ++col)
215 if (grayArray[row][col] == grayArray[row][col + hShift])
216 ++numMatches;
217
218 /* See if the number of matches exceeds the running mean plus N
219 standard deviations. Also, keep track of the shortest average
220 distance between matches seen so far.
221 */
222 if (hShift > 1) {
223 int i;
224 for (i = 2; i >= 0; --i)
225 if (bestSeparation[i] == 0 &&
226 numMatches > meanMatches + (i+1)*stdMatches) {
227 bestSeparation[i] = hShift;
228 break;
229 }
230 }
231 avgDist = (height * (width-hShift)) / (double)numMatches;
232 if (minAvgDist > avgDist) {
233 minAvgDist = avgDist;
234 altBestSeparation = hShift;
235 }
236
237 /* Compute the new mean and standard deviation. */
238 sumMatches += (double)numMatches;
239 sumSqMatches += (double)numMatches * (double)numMatches;
240 meanMatches = sumMatches / (double)hShift;
241 stdMatches = sqrt(sumSqMatches/hShift - meanMatches*meanMatches);
242 }
243 }
244
245 return bestEyeSepWeEncountered(bestSeparation, altBestSeparation);
246 }
247
248
249
250 #ifndef LITERAL_FN_DEF_MATCH
251 static qsort_comparison_fn compareInts;
252 #endif
253
254 static int
compareInts(const void * const a,const void * const b)255 compareInts(const void * const a,
256 const void * const b) {
257
258 const int * const firstP = a;
259 const int * const secondP = b;
260
261 int const first = *firstP;
262 int const second = *secondP;
263
264 int retval;
265
266 if (first < second)
267 retval = -1;
268 else if (first > second)
269 retval = +1;
270 else
271 retval = 0;
272
273 return retval;
274 }
275
276
277
278 static int
findEyeSeparation(struct pam * const pamP,gray ** const grayArray,int const minSeparation)279 findEyeSeparation( struct pam * const pamP,
280 gray ** const grayArray,
281 int const minSeparation ) {
282 /*----------------------------------------------------------------------
283 Compute the eye separation for each row of the grayscale image.
284 Ignore rows for which the eye separation could not be determined and
285 return the median of the remaining rows, aborting with an error
286 message if there are no remaining rows. Out of laziness we use
287 qsort() to help find the median; if this turns out to be a
288 performance problem, it should be replaced with a linear-time median
289 finder.
290 ------------------------------------------------------------------------*/
291 int bestSeparation; /* Best eye separation found */
292
293 /* First attempt: Find the best eye separation across the image as a
294 whole. This works well when the image consists of relatively
295 small foreground objects in front of a comparatively large
296 background plane.
297 */
298 bestSeparation =
299 findRegionEyeSeparation( grayArray, pamP->width, pamP->height );
300
301 /* Second attempt: Compute the best eye separation for each row
302 independently and return the median of the best eye
303 separations.
304 */
305 if (bestSeparation < minSeparation) {
306 int * rowSeparation; /* Per-row best separation distance */
307 unsigned int numValidRows;
308 /* Number of entries in the above (<= #rows) */
309 unsigned int row;
310
311 numValidRows = 0; /* initial value */
312
313 MALLOCARRAY_NOFAIL( rowSeparation, pamP->height );
314 for (row = 0; row < pamP->height; ++row) {
315 int const sep =
316 findRegionEyeSeparation( &grayArray[row], pamP->width, 1);
317 if (sep >= minSeparation)
318 rowSeparation[numValidRows++] = sep;
319 }
320 if (numValidRows > 0) {
321 qsort(rowSeparation, numValidRows, sizeof(int), compareInts);
322 bestSeparation = rowSeparation[numValidRows/2];
323 }
324 free( rowSeparation );
325 }
326
327 if (bestSeparation < minSeparation)
328 pm_error("Failed to determine the separation between "
329 "the left and right views");
330
331 return bestSeparation;
332 }
333
334
335
336 static void
writeAnaglyph(FILE * const ofP,gray ** const grayArray,gray const maxGrayVal,int const eyeSep,int const swapEyes,struct pam * const pamP)337 writeAnaglyph( FILE * const ofP,
338 gray ** const grayArray,
339 gray const maxGrayVal,
340 int const eyeSep,
341 int const swapEyes,
342 struct pam * const pamP) {
343 /*----------------------------------------------------------------------
344 Output an anaglyphic stereogram from the given grayscale array and
345 eye-separation value.
346 ------------------------------------------------------------------------*/
347 struct pam outPam;
348 tuple * tuplerow;
349
350 outPam.size = sizeof(struct pam);
351 outPam.len = PAM_STRUCT_SIZE(tuple_type);
352 outPam.file = ofP;
353 outPam.format = PAM_FORMAT;
354 outPam.plainformat = 0;
355 outPam.height = pamP->height;
356 outPam.width = pamP->width - eyeSep;
357 /* Avoid color bands on the left/right edges. */
358 outPam.depth = 3;
359 outPam.maxval = (sample) maxGrayVal;
360 strcpy(outPam.tuple_type, PAM_PPM_TUPLETYPE);
361
362 pnm_writepaminit( &outPam );
363
364 tuplerow = pnm_allocpamrow( &outPam );
365
366 if (swapEyes) {
367 unsigned int row;
368
369 for (row = 0; row < outPam.height; ++row) {
370 unsigned int col;
371 for (col = 0; col < outPam.width; ++col) {
372 tuplerow[col][PAM_RED_PLANE] = grayArray[row][col+eyeSep];
373 tuplerow[col][PAM_GRN_PLANE] = grayArray[row][col];
374 tuplerow[col][PAM_BLU_PLANE] = grayArray[row][col];
375 }
376 pnm_writepamrow( &outPam, tuplerow );
377 }
378 } else {
379 unsigned int row;
380 for (row = 0; row < outPam.height; ++row) {
381 unsigned int col;
382 for (col = 0; col < outPam.width; ++col) {
383 tuplerow[col][PAM_RED_PLANE] = grayArray[row][col];
384 tuplerow[col][PAM_GRN_PLANE] = grayArray[row][col+eyeSep];
385 tuplerow[col][PAM_BLU_PLANE] = grayArray[row][col+eyeSep];
386 }
387 pnm_writepamrow( &outPam, tuplerow );
388 }
389 }
390 pnm_freepamrow( tuplerow );
391 }
392
393
394
395 int
main(int argc,const char * argv[])396 main(int argc, const char *argv[]) {
397 struct pam inPam;
398 gray ** inImage;
399 int eyeSep;
400 struct cmdlineInfo cmdline;
401
402 pm_proginit( &argc, argv );
403 parseCommandLine( argc, argv, &cmdline );
404
405 inImage = readAsGray( cmdline.inputFilename, cmdline.maxGrayVal, &inPam );
406
407 if (cmdline.separation >= 0)
408 eyeSep = cmdline.separation;
409 else {
410 int const minSeparation =
411 cmdline.minSeparation > 0
412 ? cmdline.minSeparation : inPam.width / 10;
413
414 /* Minimum separation in pixels between eyes.
415 Heuristic: Eye separation must be at least 10% of image width.
416 */
417 eyeSep = findEyeSeparation ( &inPam, inImage, minSeparation );
418 }
419
420 pm_message( "Separation between left/right views = %d pixels", eyeSep );
421
422 writeAnaglyph ( stdout, inImage, cmdline.maxGrayVal,
423 eyeSep, cmdline.swapEyes,
424 &inPam );
425
426 return 0;
427 }
428
429