diff options
Diffstat (limited to 'analyzer')
| -rw-r--r-- | analyzer/Makefile | 48 | ||||
| -rw-r--r-- | analyzer/pamfile.c | 243 | ||||
| -rw-r--r-- | analyzer/pamsharpmap.c | 188 | ||||
| -rw-r--r-- | analyzer/pamsharpness.c | 153 | ||||
| -rw-r--r-- | analyzer/pamslice.c | 199 | ||||
| -rw-r--r-- | analyzer/pamsumm.c | 231 | ||||
| -rw-r--r-- | analyzer/pamtilt.c | 437 | ||||
| -rw-r--r-- | analyzer/pbmminkowski.c | 168 | ||||
| -rw-r--r-- | analyzer/pgmhist.c | 87 | ||||
| -rw-r--r-- | analyzer/pgmminkowski.c | 222 | ||||
| -rw-r--r-- | analyzer/pgmtexture.c | 1047 | ||||
| -rw-r--r-- | analyzer/pnmhistmap.c | 483 | ||||
| -rw-r--r-- | analyzer/pnmpsnr.c | 209 | ||||
| -rw-r--r-- | analyzer/ppmhist.c | 290 |
14 files changed, 4005 insertions, 0 deletions
diff --git a/analyzer/Makefile b/analyzer/Makefile new file mode 100644 index 00000000..6a447ea7 --- /dev/null +++ b/analyzer/Makefile @@ -0,0 +1,48 @@ +ifeq ($(SRCDIR)x,x) + SRCDIR = $(CURDIR)/.. + BUILDDIR = $(SRCDIR) +endif +SUBDIR = analyzer +VPATH=.:$(SRCDIR)/$(SUBDIR) + +include $(BUILDDIR)/Makefile.config + +# We tend to separate out the build targets so that we don't have +# any more dependencies for a given target than it really needs. +# That way, if there is a problem with a dependency, we can still +# successfully build all the stuff that doesn't depend upon it. +# This package is so big, it's useful even when some parts won't +# build. + +PORTBINARIES = pamfile pamslice pamsumm pamtilt \ + pgmhist pnmhistmap ppmhist pgmminkowski +MATHBINARIES = pamsharpmap pamsharpness pgmtexture pnmpsnr + +BINARIES = $(PORTBINARIES) $(MATHBINARIES) +SCRIPT = + +OBJECTS = $(BINARIES:%=%.o) + +# We don't include programs that have special library dependencies in the +# merge scheme, because we don't want those dependencies to prevent us +# from building all the other programs. + +MERGEBINARIES = $(BINARIES) +MERGE_OBJECTS = $(MERGEBINARIES:%=%.o2) + +.PHONY: all +all: $(BINARIES) + +include $(SRCDIR)/Makefile.common + +install.bin: install.bin.local +.PHONY: install.bin.local +install.bin.local: $(PKGDIR)/bin +# Remember that $(SYMLINK) might just be a copy command. + cd $(PKGDIR)/bin ; \ + $(SYMLINK) pamslice$(EXE) pgmslice + cd $(PKGDIR)/bin ; \ + $(SYMLINK) pamfile$(EXE) pnmfile + +FORCE: + diff --git a/analyzer/pamfile.c b/analyzer/pamfile.c new file mode 100644 index 00000000..efb2cbee --- /dev/null +++ b/analyzer/pamfile.c @@ -0,0 +1,243 @@ +/* pamfile.c - describe a portable anymap +** +** Copyright (C) 1991 by Jef Poskanzer. +** +** Permission to use, copy, modify, and distribute this software and its +** documentation for any purpose and without fee is hereby granted, provided +** that the above copyright notice appear in all copies and that both that +** copyright notice and this permission notice appear in supporting +** documentation. This software is provided "as is" without express or +** implied warranty. +*/ + +#include "mallocvar.h" +#include "nstring.h" +#include "shhopt.h" +#include "pam.h" + +struct cmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + int inputFileCount; /* Number of input files */ + const char ** inputFilespec; /* Filespecs of input files */ + unsigned int allimages; /* -allimages or -count */ + unsigned int count; /* -count */ + unsigned int comments; /* -comments */ +}; + + + +static void +parseCommandLine(int argc, char ** argv, + struct cmdlineInfo * const cmdlineP) { +/*---------------------------------------------------------------------------- + Note that the file spec array we return is stored in the storage that + was passed to as as the argv array. +-----------------------------------------------------------------------------*/ + optEntry * option_def; + /* Instructions to optParseOptions3 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + + MALLOCARRAY_NOFAIL(option_def, 100); + + option_def_index = 0; /* incremented by OPTENT3 */ + OPTENT3(0, "allimages", OPT_FLAG, NULL, &cmdlineP->allimages, 0); + OPTENT3(0, "count", OPT_FLAG, NULL, &cmdlineP->count, 0); + OPTENT3(0, "comments", OPT_FLAG, NULL, &cmdlineP->comments, 0); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ + opt.allowNegNum = FALSE; /* We have no parms that are negative numbers */ + + optParseOptions3(&argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others */ + + cmdlineP->inputFilespec = (const char **)&argv[1]; + cmdlineP->inputFileCount = argc - 1; +} + + + +static void +dumpHeader(struct pam const pam) { + + switch (pam.format) { + case PAM_FORMAT: + printf("PAM, %d by %d by %d maxval %ld\n", + pam.width, pam.height, pam.depth, pam.maxval); + printf(" Tuple type: %s\n", pam.tuple_type); + break; + + case PBM_FORMAT: + printf("PBM plain, %d by %d\n", pam.width, pam.height ); + break; + + case RPBM_FORMAT: + printf("PBM raw, %d by %d\n", pam.width, pam.height); + break; + + case PGM_FORMAT: + printf("PGM plain, %d by %d maxval %ld\n", + pam.width, pam.height, pam.maxval); + break; + + case RPGM_FORMAT: + printf("PGM raw, %d by %d maxval %ld\n", + pam.width, pam.height, pam.maxval); + break; + + case PPM_FORMAT: + printf("PPM plain, %d by %d maxval %ld\n", + pam.width, pam.height, pam.maxval); + break; + + case RPPM_FORMAT: + printf("PPM raw, %d by %d maxval %ld\n", + pam.width, pam.height, pam.maxval); + break; + } +} + + + +static void +dumpComments(const char * const comments) { + + const char * p; + bool startOfLine; + + printf("Comments:\n"); + + for (p = &comments[0], startOfLine = TRUE; *p; ++p) { + if (startOfLine) + printf(" #"); + + fputc(*p, stdout); + + if (*p == '\n') + startOfLine = TRUE; + else + startOfLine = FALSE; + } + if (!startOfLine) + fputc('\n', stdout); +} + + + +static void +doOneImage(const char * const name, + unsigned int const imageDoneCount, + FILE * const fileP, + bool const allimages, + bool const justCount, + bool const wantComments, + bool * const eofP) { + + struct pam pam; + const char * comments; + enum pm_check_code checkRetval; + + pam.comment_p = &comments; + + pnm_readpaminit(fileP, &pam, PAM_STRUCT_SIZE(comment_p)); + + if (!justCount) { + if (allimages) + printf("%s:\tImage %d:\t", name, imageDoneCount); + else + printf("%s:\t", name); + + dumpHeader(pam); + if (wantComments) + dumpComments(comments); + } + strfree(comments); + + pnm_checkpam(&pam, PM_CHECK_BASIC, &checkRetval); + if (allimages) { + tuple * tuplerow; + unsigned int row; + + tuplerow = pnm_allocpamrow(&pam); + + for (row = 0; row < pam.height; ++row) + pnm_readpamrow(&pam, tuplerow); + + pnm_freepamrow(tuplerow); + + pnm_nextimage(fileP, eofP); + } +} + + + +static void +describeOneFile(const char * const name, + FILE * const fileP, + bool const allimages, + bool const justCount, + bool const wantComments) { +/*---------------------------------------------------------------------------- + Describe one image stream (file). Its name, for purposes of display, + is 'name'. The file is open as *fileP and positioned to the beginning. + + 'allimages' means report on every image in the stream and read all of + every image from it, as opposed to reading just the header of the first + image and reporting just on that. + + 'justCount' means don't tell anything about the stream except how + many images are in it. Pretty useless without 'allimages'. + + 'wantComments' means to show the comments from the image header. + Meaningless with 'justCount'. +-----------------------------------------------------------------------------*/ + unsigned int imageDoneCount; + /* Number of images we've processed so far */ + bool eof; + + eof = FALSE; + imageDoneCount = 0; + + while (!eof && (imageDoneCount < 1 || allimages)) { + doOneImage(name, imageDoneCount, fileP, + allimages, justCount, wantComments, + &eof); + ++imageDoneCount; + } + if (justCount) + printf("%s:\t%u images\n", name, imageDoneCount); +} + + + +int +main(int argc, char *argv[]) { + + struct cmdlineInfo cmdline; + + pnm_init(&argc, argv); + + parseCommandLine(argc, argv, &cmdline); + + if (cmdline.inputFileCount == 0) + describeOneFile("stdin", stdin, cmdline.allimages || cmdline.count, + cmdline.count, cmdline.comments); + else { + unsigned int i; + for (i = 0; i < cmdline.inputFileCount; ++i) { + FILE * ifP; + ifP = pm_openr(cmdline.inputFilespec[i]); + describeOneFile(cmdline.inputFilespec[i], ifP, + cmdline.allimages || cmdline.count, + cmdline.count, cmdline.comments); + pm_close(ifP); + } + } + + return 0; +} diff --git a/analyzer/pamsharpmap.c b/analyzer/pamsharpmap.c new file mode 100644 index 00000000..73923ab9 --- /dev/null +++ b/analyzer/pamsharpmap.c @@ -0,0 +1,188 @@ +/*---------------------------------------------------------------------------- + Pnmsharpness +------------------------------------------------------------------------------ + + Bryan Henderson derived this (January 2004) from the program Pnmsharp + by B.W. van Schooten and distributed to Bryan under the Perl + Artistic License, as part of the Photopnmtools package. Bryan placed + his modifications in the public domain. + + This is the copyright/license notice from the original: + + Copyright (c) 2002, 2003 by B.W. van Schooten. All rights reserved. + This software is distributed under the Perl Artistic License. + No warranty. See file 'artistic.license' for more details. + + boris@13thmonkey.org + www.13thmonkey.org/~boris/photopnmtools/ +-----------------------------------------------------------------------------*/ + +#include <stdio.h> +#include <math.h> + +#include "pam.h" +#include "shhopt.h" +#include "mallocvar.h" + +struct cmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char *inputFilespec; /* '-' if stdin */ + unsigned int context; +}; + + + +static void +parseCommandLine ( int argc, char ** argv, + struct cmdlineInfo *cmdlineP ) { +/*---------------------------------------------------------------------------- + parse program command line described in Unix standard form by argc + and argv. Return the information in the options as *cmdlineP. + + If command line is internally inconsistent (invalid options, etc.), + issue error message to stderr and abort program. + + Note that the strings we return are stored in the storage that + was passed to us as the argv array. We also trash *argv. +-----------------------------------------------------------------------------*/ + optEntry *option_def = malloc(100*sizeof(optEntry)); + /* Instructions to optParseOptions3 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + + unsigned int contextSpec; + + option_def_index = 0; /* incremented by OPTENT3 */ + OPTENT3(0, "context", OPT_UINT, &cmdlineP->context, + &contextSpec, 0 ); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ + opt.allowNegNum = FALSE; /* We have no parms that are negative numbers */ + + optParseOptions3( &argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others. */ + + if (!contextSpec) + cmdlineP->context = 1; + + if (cmdlineP->context < 1) + pm_error("-context must be at least 1"); + + + if (argc-1 > 1) + pm_error("The only argument is the input file name"); + else if (argc-1 < 1) + cmdlineP->inputFilespec = "-"; + else + cmdlineP->inputFilespec = argv[1]; +} + + + +static void +makeSharpnessPixel(struct pam * const pamP, + float * const sharpness, + tuple const sharpnessTuple) { + + unsigned int plane; + for (plane = 0; plane < pamP->depth; ++plane) + sharpnessTuple[plane] = + (sample)(sharpness[plane] * pamP->maxval + 0.5); +} + + + +static void +makeBlackTuplen(struct pam * const pamP, + tuplen const tuplen) { + + unsigned int plane; + for (plane = 0; plane < pamP->depth; ++plane) + tuplen[plane] = 0.0; +} + + + +static void +makeBlackRown(struct pam * const pamP, + tuplen * const tuplenrow) { + + unsigned int col; + for (col = 0; col < pamP->width; ++col) + makeBlackTuplen(pamP, tuplenrow[col]); +} + + + +int +main(int argc, char **argv) { + + struct cmdlineInfo cmdline; + FILE * ifP; + tuplen ** tuplenarray; + struct pam inpam; + struct pam mappam; + tuple ** map; + int row; + float * sharpness; + + pnm_init(&argc, argv); + + parseCommandLine(argc, argv, &cmdline); + + ifP = pm_openr(cmdline.inputFilespec); + + tuplenarray = pnm_readpamn(ifP, &inpam, sizeof(inpam)); + + mappam = inpam; + mappam.file = stdout; + mappam.maxval = 255; + + MALLOCARRAY_NOFAIL(sharpness, inpam.depth); + + map = pnm_allocpamarray(&mappam); + makeBlackRown(&inpam, tuplenarray[0]); + for (row = 1; row < inpam.height-1; ++row) { + int col; + makeBlackTuplen(&inpam, tuplenarray[row][0]); + for (col = 1; col < inpam.width-1; ++col) { + int dy; + unsigned int plane; + + for (plane = 0; plane < inpam.depth; ++plane) + sharpness[plane] = 0.0; + + for (dy = -1; dy <= 1; ++dy) { + int dx; + for (dx = -1; dx <= 1; ++dx) { + if (dx != 0 || dy != 0) { + unsigned int plane; + for (plane = 0; plane < inpam.depth; ++plane) { + samplen const sampleval = + tuplenarray[row][col][plane]; + samplen const sampleval2 = + tuplenarray[row+dy][col+dx][plane]; + sharpness[plane] += fabs(sampleval - sampleval2); + } + } + } + } + makeSharpnessPixel(&mappam, sharpness, map[row][col]); + } + makeBlackTuplen(&inpam, tuplenarray[row][inpam.width-1]); + } + makeBlackRown(&inpam, tuplenarray[inpam.height-1]); + free(sharpness); + + pnm_writepam(&mappam, map); + + pnm_freepamarray(map, &mappam); + pnm_freepamarrayn(tuplenarray, &inpam); + + return 0; +} diff --git a/analyzer/pamsharpness.c b/analyzer/pamsharpness.c new file mode 100644 index 00000000..7e52a9ba --- /dev/null +++ b/analyzer/pamsharpness.c @@ -0,0 +1,153 @@ +/*---------------------------------------------------------------------------- + Pnmsharpness +------------------------------------------------------------------------------ + + Bryan Henderson derived this (January 2004) from the program of the + same name by B.W. van Schooten and distributed to Bryan under the Perl + Artistic License, as part of the Photopnmtools package. Bryan placed + his modifications in the public domain. + + This is the copyright/license notice from the original: + + Copyright (c) 2002, 2003 by B.W. van Schooten. All rights reserved. + This software is distributed under the Perl Artistic License. + No warranty. See file 'artistic.license' for more details. + + boris@13thmonkey.org + www.13thmonkey.org/~boris/photopnmtools/ +-----------------------------------------------------------------------------*/ + +#include <stdio.h> +#include <math.h> + +#include "pam.h" +#include "shhopt.h" + +struct cmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char *inputFilespec; /* '-' if stdin */ + unsigned int context; +}; + + + +static void +parseCommandLine ( int argc, char ** argv, + struct cmdlineInfo *cmdlineP ) { +/*---------------------------------------------------------------------------- + parse program command line described in Unix standard form by argc + and argv. Return the information in the options as *cmdlineP. + + If command line is internally inconsistent (invalid options, etc.), + issue error message to stderr and abort program. + + Note that the strings we return are stored in the storage that + was passed to us as the argv array. We also trash *argv. +-----------------------------------------------------------------------------*/ + optEntry *option_def = malloc(100*sizeof(optEntry)); + /* Instructions to optParseOptions3 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + + unsigned int contextSpec; + + option_def_index = 0; /* incremented by OPTENT3 */ + OPTENT3(0, "context", OPT_UINT, &cmdlineP->context, + &contextSpec, 0 ); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ + opt.allowNegNum = FALSE; /* We have no parms that are negative numbers */ + + optParseOptions3( &argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others. */ + + if (!contextSpec) + cmdlineP->context = 1; + + if (cmdlineP->context < 1) + pm_error("-context must be at least 1"); + + + if (argc-1 > 1) + pm_error("The only argument is the input file name"); + else if (argc-1 < 1) + cmdlineP->inputFilespec = "-"; + else + cmdlineP->inputFilespec = argv[1]; +} + + + +static void +computeSharpness(struct pam * const inpamP, + tuplen ** const tuplenarray, + double * const sharpnessP) { + + int row; + double totsharp; + + totsharp = 0.0; + + for (row = 1; row < inpamP->height-1; ++row) { + int col; + for (col = 1; col < inpamP->width-1; ++col) { + int dy; + for (dy = -1; dy <= 1; ++dy) { + int dx; + for (dx = -1; dx <= 1; ++dx) { + if (dx != 0 || dy != 0) { + unsigned int plane; + for (plane = 0; plane < inpamP->depth; ++plane) { + samplen const sampleval = + tuplenarray[row][col][plane]; + samplen const sampleval2 = + tuplenarray[row+dy][col+dx][plane]; + totsharp += fabs(sampleval - sampleval2); + } + } + } + } + } + } + *sharpnessP = + totsharp / (inpamP->width * inpamP->height * inpamP->depth * 8); + /* The 8 above is for the 8 neighbors to which we compare each pixel */ +} + + + +int +main(int argc, char **argv) { + + struct cmdlineInfo cmdline; + FILE * ifP; + tuplen ** tuplenarray; + struct pam inpam; + double sharpness; + + pnm_init(&argc, argv); + + parseCommandLine(argc, argv, &cmdline); + + ifP = pm_openr(cmdline.inputFilespec); + + tuplenarray = pnm_readpamn(ifP, &inpam, sizeof(inpam)); + + if (inpam.height < 3 || inpam.width < 3) + pm_error("sharpness is undefined for an image less than 3 pixels " + "in all directions. This image is %d x %d", + inpam.width, inpam.height); + + computeSharpness(&inpam, tuplenarray, &sharpness); + + pm_message("Sharpness = %f\n", sharpness); + + pnm_freepamarrayn(tuplenarray, &inpam); + pm_close(ifP); + return 0; +} diff --git a/analyzer/pamslice.c b/analyzer/pamslice.c new file mode 100644 index 00000000..fc63a2cc --- /dev/null +++ b/analyzer/pamslice.c @@ -0,0 +1,199 @@ +/* + * + * This program (Pamslice) was derived by Bryan Henderson in July 2002 + * from Pgmslice by Jos Dingjan. Pgmslice did the same thing, but + * only on PGM images. Pamslice is a total rewrite. + * + * Copyright (C) 2000 Jos Dingjan <jos@tuatha.org> + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ + +enum orientation {ROW, COLUMN}; + +#include "pam.h" +#include "shhopt.h" + +struct cmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char *inputFilespec; /* Filespec of input file */ + enum orientation orientation; + union { + unsigned int row; + unsigned int col; + } u; + unsigned int onePlane; + unsigned int plane; + unsigned int xmgr; +}; + + +static void +parseCommandLine(int argc, char ** const argv, + struct cmdlineInfo * const cmdlineP) { +/*---------------------------------------------------------------------------- + Note that the file spec array we return is stored in the storage that + was passed to us as the argv array. +-----------------------------------------------------------------------------*/ + optEntry *option_def = malloc(100*sizeof(optEntry)); + /* Instructions to OptParseOptions2 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + + unsigned int rowSpec, colSpec; + + option_def_index = 0; /* incremented by OPTENTRY */ + OPTENT3(0, "row", OPT_UINT, &cmdlineP->u.row, &rowSpec, 0); + OPTENT3(0, "column", OPT_UINT, &cmdlineP->u.col, &colSpec, 0); + OPTENT3(0, "plane", OPT_UINT, &cmdlineP->plane, &cmdlineP->onePlane, 0); + OPTENT3(0, "xmgr", OPT_FLAG, NULL, &cmdlineP->xmgr, 0); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ + opt.allowNegNum = FALSE; /* We have no parms that are negative numbers */ + + optParseOptions3(&argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others. */ + + if (rowSpec && colSpec) + pm_error("You cannot specify both -row and -col"); + + if (argc-1 > 1) + pm_error("Too many arguments (%d). Only argument is filename.", + argc-1); + else if (argc-1 == 1) + cmdlineP->inputFilespec = argv[1]; + else + cmdlineP->inputFilespec = "-"; + + if (rowSpec) + cmdlineP->orientation = ROW; + else if (colSpec) + cmdlineP->orientation = COLUMN; + else + pm_error("You must specify either -column or -row"); +} + + + +static void +printSlice(FILE * const outfile, + tuple ** const tuples, + unsigned int const rowstart, + unsigned int const rowend, + unsigned int const colstart, + unsigned int const colend, + unsigned int const planestart, + unsigned int const planeend, + bool const xmgr) { + + unsigned int count; + unsigned int row; + + + if (xmgr) { + fprintf(outfile,"@ title \"Graylevel\"\n"); + fprintf(outfile,"@ subtitle \"from (%d,%d) to (%d,%d)\"\n", + colstart, rowstart, colend, rowend); + if (colend - colstart == 1) + fprintf(outfile,"@ xaxis label \"row\"\n"); + else + fprintf(outfile,"@ xaxis label \"column\"\n"); + fprintf(outfile,"@ yaxis label \"grayvalue\"\n"); + if (planeend - planestart == 1) + fprintf(stdout,"@ type xy\n"); + else + fprintf(stdout,"@ type nxy\n"); + } + + count = 0; + for (row = rowstart; row < rowend; ++row) { + unsigned int col; + for (col = colstart; col < colend; ++col) { + unsigned int plane; + fprintf(outfile, "%d ", count++); + for (plane = planestart; plane < planeend; ++plane) + fprintf(outfile, "%lu ", tuples[row][col][plane]); + fprintf(outfile, "\n"); + } + } +} + + + +int +main(int argc, char *argv[]) { + + struct cmdlineInfo cmdline; + FILE *ifP; + struct pam inpam; + tuple **tuples; + unsigned int colstart, colend, rowstart, rowend, planestart, planeend; + + pgm_init( &argc, argv ); + + parseCommandLine(argc, argv, &cmdline); + + ifP = pm_openr(cmdline.inputFilespec); + + tuples = pnm_readpam(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type)); + + switch (cmdline.orientation) { + case COLUMN: + if (cmdline.u.col >= inpam.width) + pm_error("You specified column %u, but there are only %u columns " + "in the image.", cmdline.u.col, inpam.width); + + colstart = cmdline.u.col; + colend = colstart + 1; + rowstart = 0; + rowend = inpam.height; + + break; + case ROW: + if (cmdline.u.row >= inpam.height) + pm_error("You specified row %u, but there are only %u rows " + "in the image.", cmdline.u.row, inpam.height); + colstart = 0; + colend = inpam.width; + rowstart = cmdline.u.row; + rowend = rowstart + 1; + + break; + } + + if (cmdline.onePlane) { + if (cmdline.plane >= inpam.depth) + pm_error("You specified plane %u, but there are only %u planes " + "in the image.", cmdline.plane, inpam.depth); + planestart = cmdline.plane; + planeend = planestart + 1; + } else { + planestart = 0; + planeend = inpam.depth; + } + + printSlice(stdout, tuples, + rowstart, rowend, colstart, colend, planestart, planeend, + cmdline.xmgr); + + pm_close(ifP); + pm_close(stdout); + + exit(0); +} diff --git a/analyzer/pamsumm.c b/analyzer/pamsumm.c new file mode 100644 index 00000000..390b8ebb --- /dev/null +++ b/analyzer/pamsumm.c @@ -0,0 +1,231 @@ +/****************************************************************************** + pamsumm +******************************************************************************* + Summarize all the samples of a PAM image with various functions. + + By Bryan Henderson, San Jose CA 2004.02.07. + + Contributed to the public domain + + +******************************************************************************/ + +#include "pam.h" +#include "shhopt.h" +#include "mallocvar.h" + +enum function {FN_ADD, FN_MEAN, FN_MIN, FN_MAX}; + +struct cmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char *inputFilespec; /* Filespec of input file */ + enum function function; + unsigned int normalize; + unsigned int brief; + unsigned int verbose; +}; + + +static void +parseCommandLine(int argc, char ** const argv, + struct cmdlineInfo * const cmdlineP) { +/*---------------------------------------------------------------------------- + Note that the file spec array we return is stored in the storage that + was passed to us as the argv array. +-----------------------------------------------------------------------------*/ + optEntry *option_def = malloc(100*sizeof(optEntry)); + /* Instructions to OptParseOptions2 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + + unsigned int sumSpec, meanSpec, minSpec, maxSpec; + + option_def_index = 0; /* incremented by OPTENTRY */ + OPTENT3(0, "sum", OPT_FLAG, NULL, &sumSpec, 0); + OPTENT3(0, "mean", OPT_FLAG, NULL, &meanSpec, 0); + OPTENT3(0, "min", OPT_FLAG, NULL, &minSpec, 0); + OPTENT3(0, "max", OPT_FLAG, NULL, &maxSpec, 0); + OPTENT3(0, "normalize", OPT_FLAG, NULL, &cmdlineP->normalize, 0); + OPTENT3(0, "brief", OPT_FLAG, NULL, &cmdlineP->brief, 0); + OPTENT3(0, "verbose", OPT_FLAG, NULL, &cmdlineP->verbose, 0); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ + opt.allowNegNum = FALSE; /* We have no parms that are negative numbers */ + + optParseOptions3(&argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others. */ + + if (sumSpec + minSpec + maxSpec > 1) + pm_error("You may specify at most one of -sum, -min, and -max"); + + if (sumSpec) { + cmdlineP->function = FN_ADD; + } else if (meanSpec) { + cmdlineP->function = FN_MEAN; + } else if (minSpec) { + cmdlineP->function = FN_MIN; + } else if (maxSpec) { + cmdlineP->function = FN_MAX; + } else + pm_error("You must specify one of -sum, -min, or -max"); + + if (argc-1 > 1) + pm_error("Too many arguments (%d). File spec is the only argument.", + argc-1); + + if (argc-1 < 1) + cmdlineP->inputFilespec = "-"; + else + cmdlineP->inputFilespec = argv[1]; + +} + + +struct accum { + union { + double sum; + unsigned int min; + unsigned int max; + } u; +}; + + + +static void +initAccumulator(struct accum * const accumulatorP, + enum function const function) { + + switch(function) { + case FN_ADD: accumulatorP->u.sum = 0.0; break; + case FN_MEAN: accumulatorP->u.sum = 0.0; break; + case FN_MIN: accumulatorP->u.min = UINT_MAX; break; + case FN_MAX: accumulatorP->u.max = 0; break; + } +} + + + +static void +aggregate(struct pam * const inpamP, + tuple * const tupleRow, + enum function const function, + struct accum * const accumulatorP) { + + unsigned int col; + + for (col = 0; col < inpamP->width; ++col) { + unsigned int plane; + for (plane = 0; plane < inpamP->depth; ++plane) { + switch(function) { + case FN_ADD: + case FN_MEAN: + accumulatorP->u.sum += tupleRow[col][plane]; + break; + case FN_MIN: + if (tupleRow[col][plane] < accumulatorP->u.min) + accumulatorP->u.min = tupleRow[col][plane]; + break; + case FN_MAX: + if (tupleRow[col][plane] > accumulatorP->u.min) + accumulatorP->u.min = tupleRow[col][plane]; + break; + } + } + } +} + + + +static void +printSummary(struct accum const accumulator, + unsigned int const scale, + unsigned int const count, + enum function const function, + bool const normalize, + bool const brief) { + + switch(function) { + case FN_ADD: { + const char * const intro = brief ? "" : "the sum of all samples is "; + + if (normalize) + printf("%s%f\n", intro, accumulator.u.sum/scale); + else + printf("%s%u\n", intro, (unsigned int)accumulator.u.sum); + } + break; + case FN_MEAN: { + const char * const intro = brief ? "" : "the mean of all samples is "; + + if (normalize) + printf("%s%f\n", intro, accumulator.u.sum/count/scale); + else + printf("%s%f\n", intro, accumulator.u.sum/count); + } + break; + case FN_MIN: { + const char * const intro = + brief ? "" : "the minimum of all samples is "; + + if (normalize) + printf("%s%f\n", intro, (double)accumulator.u.min/scale); + else + printf("%s%u\n", intro, accumulator.u.min); + } + break; + case FN_MAX: { + const char * const intro = + brief ? "" : "the maximum of all samples is "; + + if (normalize) + printf("%s%f\n", intro, (double)accumulator.u.max/scale); + else + printf("%s%u\n", intro, accumulator.u.max); + } + break; + } +} + + + +int +main(int argc, char *argv[]) { + + FILE* ifP; + tuple* inputRow; /* Row from input image */ + int row; + struct cmdlineInfo cmdline; + struct pam inpam; /* Input PAM image */ + struct accum accumulator; + + pnm_init(&argc, argv); + + parseCommandLine(argc, argv, &cmdline); + + ifP = pm_openr(cmdline.inputFilespec); + + pnm_readpaminit(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type)); + + inputRow = pnm_allocpamrow(&inpam); + + initAccumulator(&accumulator, cmdline.function); + + for (row = 0; row < inpam.height; row++) { + pnm_readpamrow(&inpam, inputRow); + + aggregate(&inpam, inputRow, cmdline.function, &accumulator); + } + printSummary(accumulator, (unsigned)inpam.maxval, + inpam.height * inpam.width * inpam.depth, + cmdline.function, cmdline.normalize, cmdline.brief); + + pnm_freepamrow(inputRow); + pm_close(inpam.file); + + return 0; +} diff --git a/analyzer/pamtilt.c b/analyzer/pamtilt.c new file mode 100644 index 00000000..3753955b --- /dev/null +++ b/analyzer/pamtilt.c @@ -0,0 +1,437 @@ +/*============================================================================= + pgmtilt +=============================================================================== + Print the tilt angle of a PGM file + + Based on pgmskew by Gregg Townsend, August 2005. + + Adapted to Netpbm by Bryan Henderson, August 2005. + + All work has been contributed to the public domain by its authors. +=============================================================================*/ + +#include <math.h> + +#include "pam.h" +#include "mallocvar.h" +#include "shhopt.h" + +/* program constants */ +#define DIFFMAX 255 /* maximum scaling for differences */ +#define BARLENGTH 60 /* length of bars printed with -v */ + +struct cmdlineInfo { + /* command parameters */ + const char * inputFilename; + float maxangle; /* maximum angle attempted */ + float astep; /* initial angle increment */ + float qmin; /* minimum quality of solution */ + unsigned int hstep; /* horizontal step size */ + unsigned int vstep; /* vertical step size */ + unsigned int dstep; /* difference distance */ + unsigned int fast; /* skip third iteration */ + unsigned int verbose; /* generate commentary */ +}; + +static void +abandon(void) { + + printf("00.00\n"); + exit(0); +} + + + +static void +parseCommandLine(int argc, char *argv[], + struct cmdlineInfo * const cmdlineP) { + + static optEntry option_def[50]; + static optStruct3 opt; + unsigned int option_def_index; + + /* set defaults */ + cmdlineP->hstep = 11; /* read only every 11th column */ + cmdlineP->vstep = 5; /* calc differences every 5th row */ + cmdlineP->dstep = 2; /* check for differences two rows down */ + cmdlineP->maxangle = 10.0; /* assume skew is less than +/- ten degrees */ + cmdlineP->astep = 1.0; /* initially check by one-degree increments */ + cmdlineP->qmin = 1.0; /* don't require S/N better than 1.0 */ + + /* initialize option table */ + option_def_index = 0; /* incremented by OPTENT3 */ + OPTENT3(0, "fast", OPT_FLAG, NULL, &cmdlineP->fast, 0); + OPTENT3(0, "verbose", OPT_FLAG, NULL, &cmdlineP->verbose, 0); + OPTENT3(0, "angle", OPT_FLOAT, &cmdlineP->maxangle, NULL, 0); + OPTENT3(0, "quality", OPT_FLOAT, &cmdlineP->qmin, NULL, 0); + OPTENT3(0, "astep", OPT_FLOAT, &cmdlineP->astep, NULL, 0); + OPTENT3(0, "hstep", OPT_UINT, &cmdlineP->hstep, NULL, 0); + OPTENT3(0, "vstep", OPT_UINT, &cmdlineP->vstep, NULL, 0); + OPTENT3(0, "dstep", OPT_UINT, &cmdlineP->dstep, NULL, 0); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* no short options used */ + opt.allowNegNum = FALSE; /* don't allow negative values */ + optParseOptions3(&argc, argv, opt, sizeof(opt), 0); + + if (cmdlineP->hstep < 1) + pm_error("-hstep must be at least 1 column."); + if (cmdlineP->vstep < 1) + pm_error("-vstep must be at least 1 row."); + if (cmdlineP->dstep < 1) + pm_error("-dstep must be at least 1 row."); + if (cmdlineP->maxangle < 1 || cmdlineP->maxangle > 45) + pm_error("-maxangle must be between 1 and 45 degrees."); + + if (argc-1 < 1) /* if input file name given */ + cmdlineP->inputFilename = "-"; + else { + cmdlineP->inputFilename = argv[1]; + + if (argc-1 > 1) + pm_error("There is at most one argument. You specified %d", + argc-1); + } +} + + + +static void +computeSteps(const struct pam * const pamP, + unsigned int const hstepReq, + unsigned int const vstepReq, + unsigned int * const hstepP, + unsigned int * const vstepP) { +/*---------------------------------------------------------------------------- + Adjust parameters if necessary now that we have the image size +-----------------------------------------------------------------------------*/ + if (pamP->width < 10 || pamP->height < 10) + abandon(); + + if (pamP->width < 10 * hstepReq) + *hstepP = pamP->width / 10; + else + *hstepP = hstepReq; + + if (pamP->height < 10 * vstepReq) + *vstepP = pamP->height / 10; + else + *vstepP = vstepReq; +} + + + +static void +load(const struct pam * const pamP, + unsigned int const hstep, + sample *** const pixelsP, + unsigned int * const hsamplesP) { +/*---------------------------------------------------------------------------- + read file into memory, returning array of rows +-----------------------------------------------------------------------------*/ + unsigned int const hsamples = 1 + (pamP->width - 1) / hstep; + /* use only this many cols */ + + unsigned int row; + tuple * tuplerow; + sample ** pixels; + + tuplerow = pnm_allocpamrow(pamP); + + MALLOCARRAY(pixels, pamP->height); + if (pixels == NULL) + pm_error("Unable to allocate array of %u pixel rows", + pamP->height); + + if (pixels == NULL) + pm_error("Unable to allocate array of %u rows", pamP->height); + + for (row = 0; row < pamP->height; ++row) { + unsigned int i; + unsigned int col; + + pnm_readpamrow(pamP, tuplerow); + + MALLOCARRAY(pixels[row], hsamples); + if (pixels[row] == NULL) + pm_error("Unable to allocate %u-sample array for Row %u", + hsamples, row); + + /* save every hstep'th column */ + for (i = col = 0; i < hsamples; ++i, col += hstep) + pixels[row][i] = tuplerow[col][0]; + } + + pnm_freepamrow(tuplerow); + + *hsamplesP = hsamples; + *pixelsP = pixels; +} + + + +static void +freePixels(sample ** const samples, + unsigned int const rows) { + + unsigned int row; + + for (row = 0; row < rows; ++row) + free(samples[row]); + + free(samples); +} + + + +static void +replacePixelValuesWithScaledDiffs( + const struct pam * const pamP, + sample ** const pixels, + unsigned int const hsamples, + unsigned int const dstep) { +/*---------------------------------------------------------------------------- + Replace pixel values with scaled differences used for all + calculations +-----------------------------------------------------------------------------*/ + float const m = DIFFMAX / (float) pamP->maxval; /* scale multiplier */ + + unsigned int row; + + for (row = dstep; row < pamP->height; ++row) { + unsigned int col; + for (col = 0; col < hsamples; ++col) { + int const d = pixels[row - dstep][col] - pixels[row][col]; + unsigned int const absd = d < 0 ? -d : d; + pixels[row - dstep][col] = m * absd; /* scale the difference */ + } + } +} + + + +static void +scoreAngle(const struct pam * const pamP, + sample ** const pixels, + unsigned int const hstep, + unsigned int const vstep, + unsigned int const hsamples, + float const deg, + float * const scoreP) { +/*---------------------------------------------------------------------------- + calculate score for a given angle +-----------------------------------------------------------------------------*/ + float const radians = (float)deg/360 * 2 * M_PI; + float const dy = hstep * tan(radians); + int const dtotal = pamP->width * tan(radians); + double const tscale = 1.0 / hsamples; + + double total; + int first; + int last; + + unsigned int row; + + total = 0.0; /* initial value */ + + if (dtotal > 0) { + first = 0; + last = pamP->height - 1 - (dtotal + 1); + } else { + first = -(dtotal - 1); + last = pamP->height - 1; + } + for (row = first; row < last; row += vstep) { + float o; + long t; + double dt; + + unsigned int i; + + for (i = 0, t = 0, o = 0.5; + i < hsamples; + ++i, t += pixels[(int)(row + o)][i], o += dy) { + } + dt = tscale * t; + total += dt * dt; + } + *scoreP = total / (last - first); +} + + + +static void +getBestAngleLocal( + const struct pam * const pamP, + sample ** const pixels, + unsigned int const hstep, + unsigned int const vstep, + unsigned int const hsamples, + float const minangle, + float const maxangle, + float const incr, + bool const verbose, + float * const bestAngleP, + float * const qualityP) { +/*---------------------------------------------------------------------------- + find angle of highest score within a range +-----------------------------------------------------------------------------*/ + int const nsamples = ((maxangle - minangle) / incr + 1.5); + + float score; + float quality; /* signal/noise ratio of the best angle */ + float angle; + float bestangle; + float bestscore; + float total; + float others; + float * results; + int i; + + MALLOCARRAY_NOFAIL(results, nsamples); /* allocate array of results */ + + /* try all angles within the given range, stepping by incr */ + bestangle = minangle; + bestscore = 0; + total = 0; + for (i = 0; i < nsamples; i++) { + angle = minangle + i * incr; + scoreAngle(pamP, pixels, hstep, vstep, hsamples, angle, &score); + results[i] = score; + if (score > bestscore || + (score == bestscore && fabs(angle) < fabs(bestangle))) { + bestscore = score; + bestangle = angle; + } + total += score; + } + + others = (total-bestscore) / (nsamples-1); /* get mean of other scores */ + quality = bestscore / others; + + if (verbose) { + fprintf(stderr, + "\n%2d angles from %6.2f to %5.2f by %4.2f: " + "best = %6.2f, S:N = %4.2f\n", + nsamples, minangle, maxangle, incr, bestangle, quality); + for (i = 0; i < nsamples; ++i) { + float const angle = minangle + i * incr; + int const n = (int) (BARLENGTH * results[i] / bestscore + 0.5); + fprintf(stderr, "%6.2f: %8.2f %0*d\n", angle, results[i], n, 0); + } + } + free(results); + + *qualityP = quality; + *bestAngleP = bestangle; +} + + + +static void +readRelevantPixels(const char * const inputFilename, + unsigned int const hstepReq, + unsigned int const vstepReq, + unsigned int * const hstepP, + unsigned int * const vstepP, + sample *** const pixelsP, + struct pam * const pamP, + unsigned int * const hsamplesP) { +/*---------------------------------------------------------------------------- + load the image, saving only the pixels we might actually inspect +-----------------------------------------------------------------------------*/ + FILE * ifP; + unsigned int hstep; + unsigned int vstep; + + ifP = pm_openr(inputFilename); + pnm_readpaminit(ifP, pamP, PAM_STRUCT_SIZE(tuple_type)); + computeSteps(pamP, hstepReq, vstepReq, &hstep, &vstep); + + load(pamP, hstep, pixelsP, hsamplesP); + + *hstepP = hstep; + *vstepP = vstep; + + pm_close(ifP); +} + + + +static void +getAngle(const struct pam * const pamP, + sample ** const pixels, + unsigned int const hstep, + unsigned int const vstep, + unsigned int const hsamples, + float const maxangle, + float const astep, + float const qmin, + bool const fast, + bool const verbose, + float * const angleP) { + + float a; + float da; + float lastq; /* quality (s/n ratio) of last measurement */ + + getBestAngleLocal(pamP, pixels, hstep, vstep, hsamples, + -maxangle, maxangle, astep, verbose, + &a, &lastq); + + if ((a < -maxangle + astep / 2) || (a > maxangle - astep / 2)) + /* extreme val almost certainly wrong */ + abandon(); + if (lastq < qmin) + /* insufficient s/n ratio */ + abandon(); + + /* make a finer search in the neighborhood */ + da = astep / 10; + getBestAngleLocal(pamP, pixels, hstep, vstep, hsamples, + a - 9 * da, a + 9 * da, da, verbose, + &a, &lastq); + + /* iterate once more unless we don't need that much accuracy */ + if (!fast) { + da /= 10; + getBestAngleLocal(pamP, pixels, hstep, vstep, hsamples, + a - 9 * da, a + 9 * da, da, verbose, + &a, &lastq); + } + *angleP = a; +} + + + +int +main(int argc, char *argv[]) { + + struct cmdlineInfo cmdline; + struct pam pam; + sample ** pixels; /* pixel data */ + unsigned int hsamples; /* horizontal samples used */ + unsigned int hstep; /* horizontal step size */ + unsigned int vstep; /* vertical step size */ + float angle; + + pgm_init(&argc, argv); /* initialize netpbm system */ + + parseCommandLine(argc, argv, &cmdline); + + readRelevantPixels(cmdline.inputFilename, cmdline.hstep, cmdline.vstep, + &hstep, &vstep, &pixels, &pam, &hsamples); + + replacePixelValuesWithScaledDiffs(&pam, pixels, hsamples, cmdline.dstep); + + getAngle(&pam, pixels, hstep, vstep, hsamples, + cmdline.maxangle, cmdline.astep, cmdline.qmin, + cmdline.fast, cmdline.verbose, &angle); + + /* report the result on stdout */ + printf("%.2f\n", angle); + + freePixels(pixels, pam.height); + + return 0; +} diff --git a/analyzer/pbmminkowski.c b/analyzer/pbmminkowski.c new file mode 100644 index 00000000..5edce506 --- /dev/null +++ b/analyzer/pbmminkowski.c @@ -0,0 +1,168 @@ +/* pbmminkowsky.c - read a portable bitmap and calculate the Minkowski Integrals +** +** Copyright (C) 2000 by Luuk van Dijk/Mind over Matter +** +** Based on pbmlife.c, +** Copyright (C) 1988,1 1991 by Jef Poskanzer. +** +** Permission to use, copy, modify, and distribute this software and its +** documentation for any purpose and without fee is hereby granted, provided +** that the above copyright notice appear in all copies and that both that +** copyright notice and this permission notice appear in supporting +** documentation. This software is provided "as is" without express or +** implied warranty. +*/ + +#include "pbm.h" + +#define ISWHITE(x) ( (x) == PBM_WHITE ) + + +int main( int argc, char** argv ){ + + FILE* ifp; + + bit* prevrow; + bit* thisrow; + bit* tmprow; + + int row; + int col; + + int countTile=0; + int countEdgeX=0; + int countEdgeY=0; + int countVertex=0; + + int rows; + int cols; + int format; + + int area, perimeter, eulerchi; + + + /* + * parse arg and initialize + */ + + pbm_init( &argc, argv ); + + if ( argc > 2 ) + pm_usage( "[pbmfile]" ); + + if ( argc == 2 ) + ifp = pm_openr( argv[1] ); + else + ifp = stdin; + + pbm_readpbminit( ifp, &cols, &rows, &format ); + + prevrow = pbm_allocrow( cols ); + thisrow = pbm_allocrow( cols ); + + + /* first row */ + + pbm_readpbmrow( ifp, thisrow, cols, format ); + + /* tiles */ + + for ( col = 0; col < cols; ++col ) + if( ISWHITE(thisrow[col]) ) ++countTile; + + /* shortcut: for the first row, edgeY == countTile */ + countEdgeY = countTile; + + /* x-edges */ + + if( ISWHITE(thisrow[0]) ) ++countEdgeX; + + for ( col = 0; col < cols-1; ++col ) + if( ISWHITE(thisrow[col]) || ISWHITE(thisrow[col+1]) ) ++countEdgeX; + + if( ISWHITE(thisrow[cols-1]) ) ++countEdgeX; + + /* shortcut: for the first row, countVertex == countEdgeX */ + + countVertex = countEdgeX; + + + for ( row = 1; row < rows; ++row ){ + + tmprow = prevrow; + prevrow = thisrow; + thisrow = tmprow; + + pbm_readpbmrow( ifp, thisrow, cols, format ); + + /* tiles */ + + for ( col = 0; col < cols; ++col ) + if( ISWHITE(thisrow[col]) ) ++countTile; + + /* y-edges */ + + for ( col = 0; col < cols; ++col ) + if( ISWHITE(thisrow[col]) || ISWHITE( prevrow[col] )) ++countEdgeY; + + /* x-edges */ + + if( ISWHITE(thisrow[0]) ) ++countEdgeX; + + for ( col = 0; col < cols-1; ++col ) + if( ISWHITE(thisrow[col]) || ISWHITE(thisrow[col+1]) ) ++countEdgeX; + + if( ISWHITE(thisrow[cols-1]) ) ++countEdgeX; + + /* vertices */ + + if( ISWHITE(thisrow[0]) || ISWHITE(prevrow[0]) ) ++countVertex; + + for ( col = 0; col < cols-1; ++col ) + if( ISWHITE(thisrow[col]) || ISWHITE(thisrow[col+1]) + || ISWHITE(prevrow[col]) || ISWHITE(prevrow[col+1]) ) ++countVertex; + + if( ISWHITE(thisrow[cols-1]) || ISWHITE(prevrow[cols-1]) ) ++countVertex; + + + } /* for row */ + + /* now thisrow contains the top row*/ + /* tiles and x-edges have been counted, now y-edges and top vertices remain */ + + + /* y-edges */ + + for ( col = 0; col < cols; ++col ) + if( ISWHITE(thisrow[col]) ) ++countEdgeY; + + /* vertices */ + + if( ISWHITE(thisrow[0]) ) ++countVertex; + + for ( col = 0; col < cols-1; ++col ) + if( ISWHITE(thisrow[col]) || ISWHITE(thisrow[col+1]) ) ++countVertex; + + if( ISWHITE(thisrow[cols-1]) ) ++countVertex; + + + /* cleanup */ + + pm_close( ifp ); + + /* print results */ + + printf( " tiles:\t%d\n x-edges:\t%d\n y-edges:\t%d\nvertices:\t%d\n", + countTile, countEdgeX, countEdgeY,countVertex ); + + area = countTile; + perimeter = 2*countEdgeX + 2*countEdgeY - 4*countTile; + eulerchi = countTile - countEdgeX - countEdgeY + countVertex; + + printf( " area:\t%d\nperimeter:\t%d\n eulerchi:\t%d\n", + area, perimeter, eulerchi ); + + exit( 0 ); + +} /* main */ + diff --git a/analyzer/pgmhist.c b/analyzer/pgmhist.c new file mode 100644 index 00000000..8f4e512e --- /dev/null +++ b/analyzer/pgmhist.c @@ -0,0 +1,87 @@ +/* pgmhist.c - print a histogram of the values in a portable graymap +** +** Copyright (C) 1989 by Jef Poskanzer. +** +** Permission to use, copy, modify, and distribute this software and its +** documentation for any purpose and without fee is hereby granted, provided +** that the above copyright notice appear in all copies and that both that +** copyright notice and this permission notice appear in supporting +** documentation. This software is provided "as is" without express or +** implied warranty. +*/ + +#include "pgm.h" +#include "mallocvar.h" + +int +main( argc, argv ) + int argc; + char *argv[]; +{ + FILE *ifp; + gray maxval, *grayrow; + register gray *gP; + int argn, rows, cols, format, row; + int i, *hist, *rcount, count, size; + register int col; + const char * const usage = "[pgmfile]"; + + pgm_init( &argc, argv ); + + argn = 1; + + if ( argn < argc ) + { + ifp = pm_openr( argv[argn] ); + argn++; + } + else + ifp = stdin; + + if ( argn != argc ) + pm_usage( usage ); + + pgm_readpgminit( ifp, &cols, &rows, &maxval, &format ); + grayrow = pgm_allocrow( cols ); + + /* Build histogram. */ + MALLOCARRAY(hist, maxval + 1); + MALLOCARRAY(rcount, maxval + 1); + if ( hist == NULL || rcount == NULL ) + pm_error( "out of memory" ); + for ( i = 0; i <= maxval; i++ ) + hist[i] = 0; + for ( row = 0; row < rows; row++ ) + { + pgm_readpgmrow( ifp, grayrow, cols, maxval, format ); + for ( col = 0, gP = grayrow; col < cols; col++, gP++ ) + hist[(int) *gP]++; + } + + pm_close( ifp ); + + /* Compute count-down */ + count = 0; + for ( i = maxval; i >= 0; i-- ) + { + count += hist[i]; + rcount[i] = count; + } + + /* And print it. */ + printf( "value\tcount\tb%%\tw%%\n" ); + printf( "-----\t-----\t--\t--\n" ); + count = 0; + size = rows * cols; + for ( i = 0; i <= maxval; i++ ) + if ( hist[i] > 0 ) + { + count += hist[i]; + printf( + "%d\t%d\t%5.3g%%\t%5.3g%%\n", i, hist[i], + (float) count * 100.0 / size, + (float) rcount[i] * 100.0 / size ); + } + + exit( 0 ); +} diff --git a/analyzer/pgmminkowski.c b/analyzer/pgmminkowski.c new file mode 100644 index 00000000..dfb08429 --- /dev/null +++ b/analyzer/pgmminkowski.c @@ -0,0 +1,222 @@ +/* pgmminkowsky.c - read a portable graymap and calculate the Minkowski +** Integrals as a function of the threshold. +** +** Copyright (C) 2000 by Luuk van Dijk/Mind over Matter +** +** Based on pgmhist.c, +** Copyright (C) 1989 by Jef Poskanzer. +** +** Permission to use, copy, modify, and distribute this software and its +** documentation for any purpose and without fee is hereby granted, provided +** that the above copyright notice appear in all copies and that both that +** copyright notice and this permission notice appear in supporting +** documentation. This software is provided "as is" without express or +** implied warranty. +*/ + +#include "pgm.h" +#include "mallocvar.h" + + +#define MAX2(a,b) ( ( (a)>(b) ) ? (a) : (b) ) +#define MAX4(a,b,c,d) MAX2( MAX2((a),(b)), MAX2((c),(d)) ) + +int main( int argc, char** argv ){ + + FILE *ifp; + + gray maxval; + int cols, rows, format; + + gray* prevrow; + gray* thisrow; + gray* tmprow; + + int* countTile; + int* countEdgeX; + int* countEdgeY; + int* countVertex; + + int i, col, row; + + int maxtiles, maxedgex, maxedgey, maxvertex; + int area, perimeter, eulerchi; + + double l2inv, linv; + + /* + * parse arg and initialize + */ + + pgm_init( &argc, argv ); + + if ( argc > 2 ) pm_usage( "[pgmfile]" ); + + if ( argc == 2 ) + ifp = pm_openr( argv[1] ); + else + ifp = stdin; + + /* + * initialize + */ + + pgm_readpgminit( ifp, &cols, &rows, &maxval, &format ); + + prevrow = pgm_allocrow( cols ); + thisrow = pgm_allocrow( cols ); + + MALLOCARRAY(countTile , maxval + 1 ); + MALLOCARRAY(countEdgeX , maxval + 1 ); + MALLOCARRAY(countEdgeY , maxval + 1 ); + MALLOCARRAY(countVertex , maxval + 1 ); + + if (countTile == NULL || countEdgeX == NULL || countEdgeY == NULL || + countVertex == NULL) + pm_error( "out of memory" ); + + for ( i = 0; i <= maxval; i++ ) countTile[i] = 0; + for ( i = 0; i <= maxval; i++ ) countEdgeX[i] = 0; + for ( i = 0; i <= maxval; i++ ) countEdgeY[i] = 0; + for ( i = 0; i <= maxval; i++ ) countVertex[i] = 0; + + + + + /* first row */ + + pgm_readpgmrow( ifp, thisrow, cols, maxval, format ); + + /* tiles */ + + for ( col = 0; col < cols; ++col ) ++countTile[thisrow[col]]; + + /* y-edges */ + + for ( col = 0; col < cols; ++col ) ++countEdgeY[thisrow[col]]; + + /* x-edges */ + + ++countEdgeX[thisrow[0]]; + + for ( col = 0; col < cols-1; ++col ) + ++countEdgeX[ MAX2(thisrow[col], thisrow[col+1]) ]; + + ++countEdgeX[thisrow[cols-1]]; + + /* shortcut: for the first row, countVertex == countEdgeX */ + + ++countVertex[thisrow[0]]; + + for ( col = 0; col < cols-1; ++col ) + ++countVertex[ MAX2(thisrow[col], thisrow[col+1]) ]; + + ++countVertex[thisrow[cols-1]]; + + + + for ( row = 1; row < rows; ++row ){ + + tmprow = prevrow; + prevrow = thisrow; + thisrow = tmprow; + + pgm_readpgmrow( ifp, thisrow, cols, maxval, format ); + + /* tiles */ + + for ( col = 0; col < cols; ++col ) ++countTile[thisrow[col]]; + + /* y-edges */ + + for ( col = 0; col < cols; ++col ) + ++countEdgeY[ MAX2(thisrow[col], prevrow[col]) ]; + /* x-edges */ + + ++countEdgeX[thisrow[0]]; + + for ( col = 0; col < cols-1; ++col ) + ++countEdgeX[ MAX2(thisrow[col], thisrow[col+1]) ]; + + ++countEdgeX[thisrow[cols-1]]; + + /* vertices */ + + ++countVertex[ MAX2(thisrow[0],prevrow[0]) ]; + + for ( col = 0; col < cols-1; ++col ) + ++countVertex[ + MAX4(thisrow[col], thisrow[col+1], prevrow[col], prevrow[col+1]) + ]; + + ++countVertex[ MAX2(thisrow[cols-1],prevrow[cols-1]) ]; + + } /* for row */ + + /* now thisrow contains the top row*/ + + /* tiles and x-edges have been counted, now upper + y-edges and top vertices remain */ + + /* y-edges */ + + for ( col = 0; col < cols; ++col ) ++countEdgeY[ thisrow[col] ]; + + /* vertices */ + + ++countVertex[thisrow[0]]; + + for ( col = 0; col < cols-1; ++col ) + ++countVertex[ MAX2(thisrow[col],thisrow[col+1]) ]; + + ++countVertex[ thisrow[cols-1] ]; + + + /* cleanup */ + + maxtiles = rows * cols; + maxedgex = rows * (cols+1); + maxedgey = (rows+1) * cols; + maxvertex= (rows+1) * (cols+1); + + l2inv = 1.0/maxtiles; + linv = 0.5/(rows+cols); + + /* And print it. */ + printf( "#threshold\t tiles\tx-edges\ty-edges\tvertices\n" ); + printf( "#---------\t -----\t-------\t-------\t--------\n" ); + for ( i = 0; i <= maxval; i++ ){ + + if( !(countTile[i] || countEdgeX[i] || countEdgeY[i] || countVertex[i] ) ) + continue; /* skip empty slots */ + + area = maxtiles; + perimeter = 2*maxedgex + 2*maxedgey - 4*maxtiles; + eulerchi = maxtiles - maxedgex - maxedgey + maxvertex; + + printf( "%f\t%6d\t%7d\t%7d\t%8d\t%g\t%g\t%6d\n", (float) i/(1.0*maxval), + maxtiles, maxedgex, maxedgey, maxvertex, + area*l2inv, perimeter*linv, eulerchi + ); + + + maxtiles -= countTile[i]; + maxedgex -= countEdgeX[i]; + maxedgey -= countEdgeY[i]; + maxvertex-= countVertex[i]; + + /* i, countTile[i], countEdgeX[i], countEdgeY[i], countVertex[i] */ + + } + + /* these should be zero: */ + printf( "# check:\t%6d\t%7d\t%7d\t%8d\n", + maxtiles, maxedgex, maxedgey, maxvertex ); + + pm_close( ifp ); + + exit( 0 ); + +} /*main*/ + + diff --git a/analyzer/pgmtexture.c b/analyzer/pgmtexture.c new file mode 100644 index 00000000..38eab114 --- /dev/null +++ b/analyzer/pgmtexture.c @@ -0,0 +1,1047 @@ +/* pgmtexture.c - calculate textural features on a portable graymap +** +** Author: James Darrell McCauley +** Texas Agricultural Experiment Station +** Department of Agricultural Engineering +** Texas A&M University +** College Station, Texas 77843-2117 USA +** +** Code written partially taken from pgmtofs.c in the PBMPLUS package +** by Jef Poskanzer. +** +** Algorithms for calculating features (and some explanatory comments) are +** taken from: +** +** Haralick, R.M., K. Shanmugam, and I. Dinstein. 1973. Textural features +** for image classification. IEEE Transactions on Systems, Man, and +** Cybertinetics, SMC-3(6):610-621. +** +** Copyright (C) 1991 Texas Agricultural Experiment Station, employer for +** hire of James Darrell McCauley +** +** Permission to use, copy, modify, and distribute this software and its +** documentation for any purpose and without fee is hereby granted, provided +** that the above copyright notice appear in all copies and that both that +** copyright notice and this permission notice appear in supporting +** documentation. This software is provided "as is" without express or +** implied warranty. +** +** THE TEXAS AGRICULTURAL EXPERIMENT STATION (TAES) AND THE TEXAS A&M +** UNIVERSITY SYSTEM (TAMUS) MAKE NO EXPRESS OR IMPLIED WARRANTIES +** (INCLUDING BY WAY OF EXAMPLE, MERCHANTABILITY) WITH RESPECT TO ANY +** ITEM, AND SHALL NOT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL +** OR CONSEQUENTAL DAMAGES ARISING OUT OF THE POSESSION OR USE OF +** ANY SUCH ITEM. LICENSEE AND/OR USER AGREES TO INDEMNIFY AND HOLD +** TAES AND TAMUS HARMLESS FROM ANY CLAIMS ARISING OUT OF THE USE OR +** POSSESSION OF SUCH ITEMS. +** +** Modification History: +** 24 Jun 91 - J. Michael Carstensen <jmc@imsor.dth.dk> supplied fix for +** correlation function. +** +** 05 Oct 05 - Marc Breithecker <Marc.Breithecker@informatik.uni-erlangen.de> +** Fix calculation or normalizing constants for d > 1. +*/ + +#include <math.h> + +#include "pm_c_util.h" +#include "pgm.h" +#include "mallocvar.h" + +#define RADIX 2.0 +#define EPSILON 0.000000001 +#define BL "Angle " +#define F1 "Angular Second Moment " +#define F2 "Contrast " +#define F3 "Correlation " +#define F4 "Variance " +#define F5 "Inverse Diff Moment " +#define F6 "Sum Average " +#define F7 "Sum Variance " +#define F8 "Sum Entropy " +#define F9 "Entropy " +#define F10 "Difference Variance " +#define F11 "Difference Entropy " +#define F12 "Meas of Correlation-1 " +#define F13 "Meas of Correlation-2 " +#define F14 "Max Correlation Coeff " + +#define SIGN(x,y) ((y)<0 ? -fabs(x) : fabs(x)) +#define DOT fprintf(stderr,".") +#define SWAP(a,b) {y=(a);(a)=(b);(b)=y;} + + +static bool sortit = FALSE; + +static float * +vector (int nl, int nh) +{ + float *v; + + MALLOCARRAY(v, (unsigned) (nh - nl + 1)); + if (v == NULL) + pm_error("Unable to allocate memory for a vector."); + return v - nl; +} + + +static float ** +matrix (int nrl, int nrh, int ncl, int nch) + +/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */ +{ + int i; + float **m; + + /* allocate pointers to rows */ + MALLOCARRAY(m, (unsigned) (nrh - nrl + 1)); + if (m == NULL) + pm_error("Unable to allocate memory for a matrix."); + + m -= ncl; + + /* allocate rows and set pointers to them */ + for (i = nrl; i <= nrh; i++) + { + MALLOCARRAY(m[i], (unsigned) (nch - ncl + 1)); + if (m[i] == NULL) + pm_error("Unable to allocate memory for a matrix row."); + m[i] -= ncl; + } + /* return pointer to array of pointers to rows */ + return m; +} + +static void +results (const char * const c, const float * const a) +{ + int i; + + DOT; + fprintf (stdout, "%s", c); + for (i = 0; i < 4; ++i) + fprintf (stdout, "% 1.3e ", a[i]); + fprintf (stdout, "% 1.3e\n", (a[0] + a[1] + a[2] + a[3]) / 4); +} + +static void +simplesrt (int n, float arr[]) +{ + int i, j; + float a; + + for (j = 2; j <= n; j++) + { + a = arr[j]; + i = j - 1; + while (i > 0 && arr[i] > a) + { + arr[i + 1] = arr[i]; + i--; + } + arr[i + 1] = a; + } +} + +static void +mkbalanced (float **a, int n) +{ + int last, j, i; + float s, r, g, f, c, sqrdx; + + sqrdx = RADIX * RADIX; + last = 0; + while (last == 0) + { + last = 1; + for (i = 1; i <= n; i++) + { + r = c = 0.0; + for (j = 1; j <= n; j++) + if (j != i) + { + c += fabs (a[j][i]); + r += fabs (a[i][j]); + } + if (c && r) + { + g = r / RADIX; + f = 1.0; + s = c + r; + while (c < g) + { + f *= RADIX; + c *= sqrdx; + } + g = r * RADIX; + while (c > g) + { + f /= RADIX; + c /= sqrdx; + } + if ((c + r) / f < 0.95 * s) + { + last = 0; + g = 1.0 / f; + for (j = 1; j <= n; j++) + a[i][j] *= g; + for (j = 1; j <= n; j++) + a[j][i] *= f; + } + } + } + } +} + + +static void +reduction (float **a, int n) +{ + int m, j, i; + float y, x; + + for (m = 2; m < n; m++) + { + x = 0.0; + i = m; + for (j = m; j <= n; j++) + { + if (fabs (a[j][m - 1]) > fabs (x)) + { + x = a[j][m - 1]; + i = j; + } + } + if (i != m) + { + for (j = m - 1; j <= n; j++) + SWAP (a[i][j], a[m][j]) + for (j = 1; j <= n; j++) + SWAP (a[j][i], a[j][m]) + a[j][i] = a[j][i]; + } + if (x) + { + for (i = m + 1; i <= n; i++) + { + if ((y = a[i][m - 1])) + { + y /= x; + a[i][m - 1] = y; + for (j = m; j <= n; j++) + a[i][j] -= y * a[m][j]; + for (j = 1; j <= n; j++) + a[j][m] += y * a[j][i]; + } + } + } + } +} + + + +static void +hessenberg (float **a, int n, float wr[], float wi[]) + +{ + int nn, m, l, k, j, its, i, mmin; + float z, y, x, w, v, u, t, s, r, q, p, anorm; + + anorm = fabs (a[1][1]); + for (i = 2; i <= n; i++) + for (j = (i - 1); j <= n; j++) + anorm += fabs (a[i][j]); + nn = n; + t = 0.0; + while (nn >= 1) + { + its = 0; + do + { + for (l = nn; l >= 2; l--) + { + s = fabs (a[l - 1][l - 1]) + fabs (a[l][l]); + if (s == 0.0) + s = anorm; + if ((float) (fabs (a[l][l - 1]) + s) == s) + break; + } + x = a[nn][nn]; + if (l == nn) + { + wr[nn] = x + t; + wi[nn--] = 0.0; + } + else + { + y = a[nn - 1][nn - 1]; + w = a[nn][nn - 1] * a[nn - 1][nn]; + if (l == (nn - 1)) + { + p = 0.5 * (y - x); + q = p * p + w; + z = sqrt (fabs (q)); + x += t; + if (q >= 0.0) + { + z = p + SIGN (z, p); + wr[nn - 1] = wr[nn] = x + z; + if (z) + wr[nn] = x - w / z; + wi[nn - 1] = wi[nn] = 0.0; + } + else + { + wr[nn - 1] = wr[nn] = x + p; + wi[nn - 1] = -(wi[nn] = z); + } + nn -= 2; + } + else + { + if (its == 30) + pm_error("Too many iterations to required " + "to find %s. Giving up", F14); + if (its == 10 || its == 20) + { + t += x; + for (i = 1; i <= nn; i++) + a[i][i] -= x; + s = fabs (a[nn][nn - 1]) + fabs (a[nn - 1][nn - 2]); + y = x = 0.75 * s; + w = -0.4375 * s * s; + } + ++its; + for (m = (nn - 2); m >= l; m--) + { + z = a[m][m]; + r = x - z; + s = y - z; + p = (r * s - w) / a[m + 1][m] + a[m][m + 1]; + q = a[m + 1][m + 1] - z - r - s; + r = a[m + 2][m + 1]; + s = fabs (p) + fabs (q) + fabs (r); + p /= s; + q /= s; + r /= s; + if (m == l) + break; + u = fabs (a[m][m - 1]) * (fabs (q) + fabs (r)); + v = fabs (p) * (fabs (a[m - 1][m - 1]) + fabs (z) + + fabs (a[m + 1][m + 1])); + if ((float) (u + v) == v) + break; + } + for (i = m + 2; i <= nn; i++) + { + a[i][i - 2] = 0.0; + if (i != (m + 2)) + a[i][i - 3] = 0.0; + } + for (k = m; k <= nn - 1; k++) + { + if (k != m) + { + p = a[k][k - 1]; + q = a[k + 1][k - 1]; + r = 0.0; + if (k != (nn - 1)) + r = a[k + 2][k - 1]; + if ((x = fabs (p) + fabs (q) + fabs (r))) + { + p /= x; + q /= x; + r /= x; + } + } + if ((s = SIGN (sqrt (p * p + q * q + r * r), p))) + { + if (k == m) + { + if (l != m) + a[k][k - 1] = -a[k][k - 1]; + } + else + a[k][k - 1] = -s * x; + p += s; + x = p / s; + y = q / s; + z = r / s; + q /= p; + r /= p; + for (j = k; j <= nn; j++) + { + p = a[k][j] + q * a[k + 1][j]; + if (k != (nn - 1)) + { + p += r * a[k + 2][j]; + a[k + 2][j] -= p * z; + } + a[k + 1][j] -= p * y; + a[k][j] -= p * x; + } + mmin = nn < k + 3 ? nn : k + 3; + for (i = l; i <= mmin; i++) + { + p = x * a[i][k] + y * a[i][k + 1]; + if (k != (nn - 1)) + { + p += z * a[i][k + 2]; + a[i][k + 2] -= p * r; + } + a[i][k + 1] -= p * q; + a[i][k] -= p; + } + } + } + } + } + } while (l < nn - 1); + } +} + + + +static float +f1_asm (float **P, int Ng) + +/* Angular Second Moment */ +{ + int i, j; + float sum = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + sum += P[i][j] * P[i][j]; + + return sum; + + /* + * The angular second-moment feature (ASM) f1 is a measure of homogeneity + * of the image. In a homogeneous image, there are very few dominant + * gray-tone transitions. Hence the P matrix for such an image will have + * fewer entries of large magnitude. + */ +} + + +static float +f2_contrast (float **P, int Ng) + +/* Contrast */ +{ + int i, j, n; + float sum = 0, bigsum = 0; + + for (n = 0; n < Ng; ++n) + { + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + if ((i - j) == n || (j - i) == n) + sum += P[i][j]; + bigsum += n * n * sum; + + sum = 0; + } + return bigsum; + + /* + * The contrast feature is a difference moment of the P matrix and is a + * measure of the contrast or the amount of local variations present in an + * image. + */ +} + +static float +f3_corr (float **P, int Ng) + +/* Correlation */ +{ + int i, j; + float sum_sqrx = 0, sum_sqry = 0, tmp, *px; + float meanx =0 , meany = 0 , stddevx, stddevy; + + px = vector (0, Ng); + for (i = 0; i < Ng; ++i) + px[i] = 0; + + /* + * px[i] is the (i-1)th entry in the marginal probability matrix obtained + * by summing the rows of p[i][j] + */ + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + px[i] += P[i][j]; + + + /* Now calculate the means and standard deviations of px and py */ + /*- fix supplied by J. Michael Christensen, 21 Jun 1991 */ + /*- further modified by James Darrell McCauley, 16 Aug 1991 + * after realizing that meanx=meany and stddevx=stddevy + */ + for (i = 0; i < Ng; ++i) + { + meanx += px[i]*i; + sum_sqrx += px[i]*i*i; + } + meany = meanx; + sum_sqry = sum_sqrx; + stddevx = sqrt (sum_sqrx - (meanx * meanx)); + stddevy = stddevx; + + /* Finally, the correlation ... */ + for (tmp = 0, i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + tmp += i*j*P[i][j]; + + return (tmp - meanx * meany) / (stddevx * stddevy); + /* + * This correlation feature is a measure of gray-tone linear-dependencies + * in the image. + */ +} + + +static float +f4_var (float **P, int Ng) + +/* Sum of Squares: Variance */ +{ + int i, j; + float mean = 0, var = 0; + + /*- Corrected by James Darrell McCauley, 16 Aug 1991 + * calculates the mean intensity level instead of the mean of + * cooccurrence matrix elements + */ + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + mean += i * P[i][j]; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + var += (i + 1 - mean) * (i + 1 - mean) * P[i][j]; + + return var; +} + +static float +f5_idm (float **P, int Ng) + +/* Inverse Difference Moment */ +{ + int i, j; + float idm = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + idm += P[i][j] / (1 + (i - j) * (i - j)); + + return idm; +} + +static float +Pxpy[2 * PGM_MAXMAXVAL]; + +static float +f6_savg (float **P, int Ng) + +/* Sum Average */ +{ + int i, j; + float savg = 0; + + for (i = 0; i <= 2 * Ng; ++i) + Pxpy[i] = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + Pxpy[i + j + 2] += P[i][j]; + for (i = 2; i <= 2 * Ng; ++i) + savg += i * Pxpy[i]; + + return savg; +} + + +static float +f7_svar (float **P, int Ng, float S) { +/* Sum Variance */ + int i, j; + float var = 0; + + for (i = 0; i <= 2 * Ng; ++i) + Pxpy[i] = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + Pxpy[i + j + 2] += P[i][j]; + + for (i = 2; i <= 2 * Ng; ++i) + var += (i - S) * (i - S) * Pxpy[i]; + + return var; +} + +static float +f8_sentropy (float **P, int Ng) + +/* Sum Entropy */ +{ + int i, j; + float sentropy = 0; + + for (i = 0; i <= 2 * Ng; ++i) + Pxpy[i] = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + Pxpy[i + j + 2] += P[i][j]; + + for (i = 2; i <= 2 * Ng; ++i) + sentropy -= Pxpy[i] * log10 (Pxpy[i] + EPSILON); + + return sentropy; +} + + +static float +f9_entropy (float **P, int Ng) + +/* Entropy */ +{ + int i, j; + float entropy = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + entropy += P[i][j] * log10 (P[i][j] + EPSILON); + + return -entropy; +} + + +static float +f10_dvar (float **P, int Ng) + +/* Difference Variance */ +{ + int i, j, tmp; + float sum = 0, sum_sqr = 0, var = 0; + + for (i = 0; i <= 2 * Ng; ++i) + Pxpy[i] = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + Pxpy[abs (i - j)] += P[i][j]; + + /* Now calculate the variance of Pxpy (Px-y) */ + for (i = 0; i < Ng; ++i) + { + sum += Pxpy[i]; + sum_sqr += Pxpy[i] * Pxpy[i]; + } + tmp = Ng * Ng; + var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp); + + return var; +} + +static float +f11_dentropy (float **P, int Ng) + +/* Difference Entropy */ +{ + int i, j; + float sum = 0; + + for (i = 0; i <= 2 * Ng; ++i) + Pxpy[i] = 0; + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + Pxpy[abs (i - j)] += P[i][j]; + + for (i = 0; i < Ng; ++i) + sum += Pxpy[i] * log10 (Pxpy[i] + EPSILON); + + return -sum; +} + +static float +f12_icorr (float **P, int Ng) + +/* Information Measures of Correlation */ +{ + int i, j; + float *px, *py; + float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0; + + px = vector (0, Ng); + py = vector (0, Ng); + + /* + * px[i] is the (i-1)th entry in the marginal probability matrix obtained + * by summing the rows of p[i][j] + */ + for (i = 0; i < Ng; ++i) + { + for (j = 0; j < Ng; ++j) + { + px[i] += P[i][j]; + py[j] += P[i][j]; + } + } + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + { + hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON); + hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON); + hxy -= P[i][j] * log10 (P[i][j] + EPSILON); + } + + /* Calculate entropies of px and py - is this right? */ + for (i = 0; i < Ng; ++i) + { + hx -= px[i] * log10 (px[i] + EPSILON); + hy -= py[i] * log10 (py[i] + EPSILON); + } +/* fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */ + return ((hxy - hxy1) / (hx > hy ? hx : hy)); +} + +static float +f13_icorr (float **P, int Ng) + +/* Information Measures of Correlation */ +{ + int i, j; + float *px, *py; + float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0; + + px = vector (0, Ng); + py = vector (0, Ng); + + /* + * px[i] is the (i-1)th entry in the marginal probability matrix obtained + * by summing the rows of p[i][j] + */ + for (i = 0; i < Ng; ++i) + { + for (j = 0; j < Ng; ++j) + { + px[i] += P[i][j]; + py[j] += P[i][j]; + } + } + + for (i = 0; i < Ng; ++i) + for (j = 0; j < Ng; ++j) + { + hxy1 -= P[i][j] * log10 (px[i] * py[j] + EPSILON); + hxy2 -= px[i] * py[j] * log10 (px[i] * py[j] + EPSILON); + hxy -= P[i][j] * log10 (P[i][j] + EPSILON); + } + + /* Calculate entropies of px and py */ + for (i = 0; i < Ng; ++i) + { + hx -= px[i] * log10 (px[i] + EPSILON); + hy -= py[i] * log10 (py[i] + EPSILON); + } + /* fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */ + return (sqrt (fabs (1 - exp (-2.0 * (hxy2 - hxy))))); +} + +static float +f14_maxcorr (float **P, int Ng) + +/* Returns the Maximal Correlation Coefficient */ +{ + int i, j, k; + float *px, *py, **Q; + float *x, *iy, tmp; + + px = vector (0, Ng); + py = vector (0, Ng); + Q = matrix (1, Ng + 1, 1, Ng + 1); + x = vector (1, Ng); + iy = vector (1, Ng); + + /* + * px[i] is the (i-1)th entry in the marginal probability matrix obtained + * by summing the rows of p[i][j] + */ + for (i = 0; i < Ng; ++i) + { + for (j = 0; j < Ng; ++j) + { + px[i] += P[i][j]; + py[j] += P[i][j]; + } + } + + /* Find the Q matrix */ + for (i = 0; i < Ng; ++i) + { + for (j = 0; j < Ng; ++j) + { + Q[i + 1][j + 1] = 0; + for (k = 0; k < Ng; ++k) + Q[i + 1][j + 1] += P[i][k] * P[j][k] / px[i] / py[k]; + } + } + + /* Balance the matrix */ + mkbalanced (Q, Ng); + /* Reduction to Hessenberg Form */ + reduction (Q, Ng); + /* Finding eigenvalue for nonsymetric matrix using QR algorithm */ + hessenberg (Q, Ng, x, iy); + if (sortit) + simplesrt(Ng,x); + /* Returns the sqrt of the second largest eigenvalue of Q */ + for (i = 2, tmp = x[1]; i <= Ng; ++i) + tmp = (tmp > x[i]) ? tmp : x[i]; + return sqrt (x[Ng - 1]); +} + +int +main (int argc, char *argv[]) { + FILE *ifp; + register gray **grays; + int tone[PGM_MAXMAXVAL], R0, R45, R90, angle, d = 1, x, y; + int argn, rows, cols, row, col; + int itone, jtone, tones; + float **P_matrix0, **P_matrix45, **P_matrix90, **P_matrix135; + float ASM[4], contrast[4], corr[4], var[4], idm[4], savg[4]; + float sentropy[4], svar[4], entropy[4], dvar[4], dentropy[4]; + float icorr[4], maxcorr[4]; + gray maxval; + const char * const usage = "[-d <d>] [pgmfile]"; + + + pgm_init( &argc, argv ); + + argn = 1; + + /* Check for flags. */ + if ( argn < argc && argv[argn][0] == '-' ) + { + if ( argv[argn][1] == 'd' ) + { + ++argn; + if ( argn == argc || sscanf( argv[argn], "%d", &d ) != 1 ) + pm_usage( usage ); + } + else + pm_usage( usage ); + ++argn; + } + + if ( argn < argc ) + { + ifp = pm_openr( argv[argn] ); + ++argn; + } + else + ifp = stdin; + + if ( argn != argc ) + pm_usage( usage ); + + grays = pgm_readpgm (ifp, &cols, &rows, &maxval); + pm_close (ifp); + + if (maxval > PGM_MAXMAXVAL) + pm_error("The maxval of the image (%d) is too high. \n" + "This program's maximum is %d.", maxval, PGM_MAXMAXVAL); + + /* Determine the number of different gray scales (not maxval) */ + for (row = PGM_MAXMAXVAL; row >= 0; --row) + tone[row] = -1; + for (row = rows - 1; row >= 0; --row) + for (col = 0; col < cols; ++col) + tone[grays[row][col]] = grays[row][col]; + for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row) + if (tone[row] != -1) + tones++; + pm_message("(Image has %d graylevels.)", tones); + + /* Collapse array, taking out all zero values */ + for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++) + if (tone[row] != -1) + tone[itone++] = tone[row]; + /* Now array contains only the gray levels present (in ascending order) */ + + /* Allocate memory for gray-tone spatial dependence matrix */ + P_matrix0 = matrix (0, tones, 0, tones); + P_matrix45 = matrix (0, tones, 0, tones); + P_matrix90 = matrix (0, tones, 0, tones); + P_matrix135 = matrix (0, tones, 0, tones); + for (row = 0; row < tones; ++row) + for (col = 0; col < tones; ++col) + { + P_matrix0[row][col] = P_matrix45[row][col] = 0; + P_matrix90[row][col] = P_matrix135[row][col] = 0; + } + + /* Find gray-tone spatial dependence matrix */ + fprintf (stderr, "(Computing spatial dependence matrix..."); + for (row = 0; row < rows; ++row) + for (col = 0; col < cols; ++col) + for (x = 0, angle = 0; angle <= 135; angle += 45) + { + while (tone[x] != grays[row][col]) + x++; + if (angle == 0 && col + d < cols) + { + y = 0; + while (tone[y] != grays[row][col + d]) + y++; + P_matrix0[x][y]++; + P_matrix0[y][x]++; + } + if (angle == 90 && row + d < rows) + { + y = 0; + while (tone[y] != grays[row + d][col]) + y++; + P_matrix90[x][y]++; + P_matrix90[y][x]++; + } + if (angle == 45 && row + d < rows && col - d >= 0) + { + y = 0; + while (tone[y] != grays[row + d][col - d]) + y++; + P_matrix45[x][y]++; + P_matrix45[y][x]++; + } + if (angle == 135 && row + d < rows && col + d < cols) + { + y = 0; + while (tone[y] != grays[row + d][col + d]) + y++; + P_matrix135[x][y]++; + P_matrix135[y][x]++; + } + } + /* Gray-tone spatial dependence matrices are complete */ + + /* Find normalizing constants */ + R0 = 2 * rows * (cols - d); + R45 = 2 * (rows - d) * (cols - d); + R90 = 2 * (rows - d) * cols; + + /* Normalize gray-tone spatial dependence matrix */ + for (itone = 0; itone < tones; ++itone) + for (jtone = 0; jtone < tones; ++jtone) + { + P_matrix0[itone][jtone] /= R0; + P_matrix45[itone][jtone] /= R45; + P_matrix90[itone][jtone] /= R90; + P_matrix135[itone][jtone] /= R45; + } + + fprintf (stderr, " done.)\n"); + fprintf (stderr, "(Computing textural features"); + fprintf (stdout, "\n"); + DOT; + fprintf (stdout, + "%s 0 45 90 135 Avg\n", + BL); + + ASM[0] = f1_asm (P_matrix0, tones); + ASM[1] = f1_asm (P_matrix45, tones); + ASM[2] = f1_asm (P_matrix90, tones); + ASM[3] = f1_asm (P_matrix135, tones); + results (F1, ASM); + + contrast[0] = f2_contrast (P_matrix0, tones); + contrast[1] = f2_contrast (P_matrix45, tones); + contrast[2] = f2_contrast (P_matrix90, tones); + contrast[3] = f2_contrast (P_matrix135, tones); + results (F2, contrast); + + + corr[0] = f3_corr (P_matrix0, tones); + corr[1] = f3_corr (P_matrix45, tones); + corr[2] = f3_corr (P_matrix90, tones); + corr[3] = f3_corr (P_matrix135, tones); + results (F3, corr); + + var[0] = f4_var (P_matrix0, tones); + var[1] = f4_var (P_matrix45, tones); + var[2] = f4_var (P_matrix90, tones); + var[3] = f4_var (P_matrix135, tones); + results (F4, var); + + + idm[0] = f5_idm (P_matrix0, tones); + idm[1] = f5_idm (P_matrix45, tones); + idm[2] = f5_idm (P_matrix90, tones); + idm[3] = f5_idm (P_matrix135, tones); + results (F5, idm); + + savg[0] = f6_savg (P_matrix0, tones); + savg[1] = f6_savg (P_matrix45, tones); + savg[2] = f6_savg (P_matrix90, tones); + savg[3] = f6_savg (P_matrix135, tones); + results (F6, savg); + + sentropy[0] = f8_sentropy (P_matrix0, tones); + sentropy[1] = f8_sentropy (P_matrix45, tones); + sentropy[2] = f8_sentropy (P_matrix90, tones); + sentropy[3] = f8_sentropy (P_matrix135, tones); + svar[0] = f7_svar (P_matrix0, tones, sentropy[0]); + svar[1] = f7_svar (P_matrix45, tones, sentropy[1]); + svar[2] = f7_svar (P_matrix90, tones, sentropy[2]); + svar[3] = f7_svar (P_matrix135, tones, sentropy[3]); + results (F7, svar); + results (F8, sentropy); + + entropy[0] = f9_entropy (P_matrix0, tones); + entropy[1] = f9_entropy (P_matrix45, tones); + entropy[2] = f9_entropy (P_matrix90, tones); + entropy[3] = f9_entropy (P_matrix135, tones); + results (F9, entropy); + + dvar[0] = f10_dvar (P_matrix0, tones); + dvar[1] = f10_dvar (P_matrix45, tones); + dvar[2] = f10_dvar (P_matrix90, tones); + dvar[3] = f10_dvar (P_matrix135, tones); + results (F10, dvar); + + dentropy[0] = f11_dentropy (P_matrix0, tones); + dentropy[1] = f11_dentropy (P_matrix45, tones); + dentropy[2] = f11_dentropy (P_matrix90, tones); + dentropy[3] = f11_dentropy (P_matrix135, tones); + results (F11, dentropy); + + icorr[0] = f12_icorr (P_matrix0, tones); + icorr[1] = f12_icorr (P_matrix45, tones); + icorr[2] = f12_icorr (P_matrix90, tones); + icorr[3] = f12_icorr (P_matrix135, tones); + results (F12, icorr); + + icorr[0] = f13_icorr (P_matrix0, tones); + icorr[1] = f13_icorr (P_matrix45, tones); + icorr[2] = f13_icorr (P_matrix90, tones); + icorr[3] = f13_icorr (P_matrix135, tones); + results (F13, icorr); + + maxcorr[0] = f14_maxcorr (P_matrix0, tones); + maxcorr[1] = f14_maxcorr (P_matrix45, tones); + maxcorr[2] = f14_maxcorr (P_matrix90, tones); + maxcorr[3] = f14_maxcorr (P_matrix135, tones); + results (F14, maxcorr); + + + fprintf (stderr, " done.)\n"); + + return 0; +} diff --git a/analyzer/pnmhistmap.c b/analyzer/pnmhistmap.c new file mode 100644 index 00000000..0b1eeb1f --- /dev/null +++ b/analyzer/pnmhistmap.c @@ -0,0 +1,483 @@ +/* pnmhistmap.c - + * Draw a histogram for a PGM or PPM file + * + * Options: -verbose: the usual + * -max N: force scaling value to N + * -black: ignore all-black count + * -white: ignore all-white count + * + * - PGM histogram is a PBM file, PPM histogram is a PPM file + * - No conditional code - assumes all three: PBM, PGM, PPM + * + * Copyright (C) 1993 by Wilson H. Bent, Jr (whb@usc.edu) + * + * 2004-12-11 john h. dubois iii (john@armory.com) + * - Added options: + * -dots, -nmax, -red, -green, -blue, -width, -height, -lval, -rval + * - Deal properly with maxvals other than 256 + */ + +#include <string.h> + +#include "pnm.h" +#include "shhopt.h" +#include "mallocvar.h" + +static double const epsilon = .00001; + +#define SCALE_H(value) (hscale_unity ? (value) : (int)((value) * hscale)) + +enum wantedColor {WANT_RED=0, WANT_GRN=1, WANT_BLU=2}; + +struct cmdlineInfo { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char *inputFilespec; /* Filespecs of input files */ + unsigned int black; + unsigned int white; + unsigned int dots; + bool colorWanted[3]; + /* subscript is enum wantedColor */ + unsigned int verbose; + unsigned int nmaxSpec; + float nmax; + unsigned int lval; + unsigned int rval; + unsigned int widthSpec; + unsigned int width; + unsigned int height; +}; + + + +static void +parseCommandLine(int argc, char ** argv, + struct cmdlineInfo *cmdlineP) { +/*---------------------------------------------------------------------------- + Note that the file spec array we return is stored in the storage that + was passed to us as the argv array. +-----------------------------------------------------------------------------*/ + optEntry *option_def; + /* Instructions to optParseOptions3 on how to parse our options. + */ + optStruct3 opt; + + unsigned int option_def_index; + unsigned int lvalSpec, rvalSpec, heightSpec; + unsigned int redSpec, greenSpec, blueSpec; + + MALLOCARRAY_NOFAIL(option_def, 100); + + option_def_index = 0; /* incremented by OPTENT3 */ + OPTENT3(0, "black", OPT_FLAG, NULL, &cmdlineP->black, 0); + OPTENT3(0, "white", OPT_FLAG, NULL, &cmdlineP->white, 0); + OPTENT3(0, "dots", OPT_FLAG, NULL, &cmdlineP->dots, 0); + OPTENT3(0, "red", OPT_FLAG, NULL, &redSpec, 0); + OPTENT3(0, "green", OPT_FLAG, NULL, &greenSpec, 0); + OPTENT3(0, "blue", OPT_FLAG, NULL, &blueSpec, 0); + OPTENT3(0, "verbose", OPT_FLAG, NULL, &cmdlineP->verbose, 0); + OPTENT3(0, "nmax", OPT_FLOAT, &cmdlineP->nmax, + &cmdlineP->nmaxSpec, 0); + OPTENT3(0, "lval", OPT_UINT, &cmdlineP->lval, + &lvalSpec, 0); + OPTENT3(0, "rval", OPT_UINT, &cmdlineP->rval, + &rvalSpec, 0); + OPTENT3(0, "width", OPT_UINT, &cmdlineP->width, + &cmdlineP->widthSpec, 0); + OPTENT3(0, "height", OPT_UINT, &cmdlineP->height, + &heightSpec, 0); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ + opt.allowNegNum = FALSE; /* We may have parms that are negative numbers */ + + optParseOptions3(&argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others. */ + + if (!lvalSpec) + cmdlineP->lval = 0; + if (!rvalSpec) + cmdlineP->rval = PNM_OVERALLMAXVAL; + + if (!redSpec && !greenSpec && !blueSpec) { + cmdlineP->colorWanted[WANT_RED] = TRUE; + cmdlineP->colorWanted[WANT_GRN] = TRUE; + cmdlineP->colorWanted[WANT_BLU] = TRUE; + } else { + cmdlineP->colorWanted[WANT_RED] = redSpec; + cmdlineP->colorWanted[WANT_GRN] = greenSpec; + cmdlineP->colorWanted[WANT_BLU] = blueSpec; + } + + if (!heightSpec) + cmdlineP->height = 200; + + if (argc-1 == 0) + cmdlineP->inputFilespec = "-"; + else if (argc-1 != 1) + pm_error("Program takes zero or one argument (filename). You " + "specified %d", argc-1); + else + cmdlineP->inputFilespec = argv[1]; +} + + + +static unsigned int +maxSlotCount(const unsigned int * const hist, + unsigned int const hist_width, + bool const no_white, + bool const no_black) { +/*---------------------------------------------------------------------------- + Return the maximum count among all the slots in hist[], not counting + the first and last as suggested by 'no_white' and 'no_black'. +-----------------------------------------------------------------------------*/ + unsigned int hmax; + unsigned int i; + + unsigned int const start = (no_black ? 1 : 0); + unsigned int const finish = (no_white ? hist_width - 1 : hist_width); + for (hmax = 0, i = start; i < finish; ++i) + if (hmax < hist[i]) + hmax = hist[i]; + + return hmax; +} + + + +static void +clipHistogram(unsigned int * const hist, + unsigned int const hist_width, + unsigned int const hmax) { + + unsigned int i; + + for (i = 0; i < hist_width; ++i) + hist[i] = MIN(hmax, hist[i]); +} + + + +static void +pgm_hist(FILE * const fp, + int const cols, + int const rows, + xelval const maxval, + int const format, + bool const dots, + bool const no_white, + bool const no_black, + bool const verbose, + xelval const startval, + xelval const endval, + unsigned int const hist_width, + unsigned int const hist_height, + bool const clipSpec, + unsigned int const clipCount, + double const hscale) { + + bool const hscale_unity = hscale - 1 < epsilon; + + gray * grayrow; + bit ** bits; + int i, j; + unsigned int * ghist; + double vscale; + unsigned int hmax; + + if ((ghist = calloc(hist_width, sizeof(int))) == NULL) + pm_error ("Not enough memory for histogram array (%d bytes)", + hist_width * sizeof(int)); + if ((bits = pbm_allocarray (hist_width, hist_height)) == NULL) + pm_error ("no space for output array (%d bits)", + hist_width * hist_height); + memset (ghist, 0, sizeof (ghist)); + + /* read the pixel values into the histogram arrays */ + grayrow = pgm_allocrow (cols); + /*XX error-check! */ + if (verbose) pm_message ("making histogram..."); + for (i = rows; i > 0; --i) { + pgm_readpgmrow (fp, grayrow, cols, maxval, format); + for (j = cols-1; j >= 0; --j) { + int value; + + if ((value = grayrow[j]) >= startval && value <= endval) + ghist[SCALE_H(value-startval)]++; + } + } + pgm_freerow (grayrow); + fclose (fp); + + /* find the highest-valued slot and set the vertical scale value */ + if (verbose) + pm_message ("finding max. slot height..."); + if (clipSpec) + hmax = clipCount; + else + hmax = maxSlotCount(ghist, hist_width, no_white, no_black); + + if (verbose) + pm_message ("Done: height = %u", hmax); + + clipHistogram(ghist, hist_width, hmax); + + vscale = (double) hist_height / hmax; + + for (i = 0; i < hist_width; ++i) { + int mark = hist_height - (int)(vscale * ghist[i]); + for (j = 0; j < mark; ++j) + bits[j][i] = PBM_BLACK; + if (j < hist_height) + bits[j++][i] = PBM_WHITE; + for ( ; j < hist_height; ++j) + bits[j][i] = dots ? PBM_BLACK : PBM_WHITE; + } + + pbm_writepbm (stdout, bits, hist_width, hist_height, 0); +} + + + +static unsigned int +maxSlotCountAll(unsigned int * const hist[3], + unsigned int const hist_width, + bool const no_white, + bool const no_black) { +/*---------------------------------------------------------------------------- + Return the maximum count among all the slots in hist[x] not + counting the first and last as suggested by 'no_white' and + 'no_black'. hist[x] may be NULL to indicate none. +-----------------------------------------------------------------------------*/ + unsigned int hmax; + unsigned int color; + + hmax = 0; + + for (color = 0; color < 3; ++color) + if (hist[color]) + hmax = MAX(hmax, + maxSlotCount(hist[color], + hist_width, no_white, no_black)); + + return hmax; +} + + + +static void +createHist(bool const colorWanted[3], + unsigned int const hist_width, + unsigned int * (* const histP)[3]) { +/*---------------------------------------------------------------------------- + Allocate the histogram arrays and set each slot count to zero. +-----------------------------------------------------------------------------*/ + unsigned int color; + + for (color = 0; color < 3; ++color) + if (colorWanted[color]) { + unsigned int * hist; + unsigned int i; + MALLOCARRAY(hist, hist_width); + if (hist == NULL) + pm_error ("Not enough memory for histogram arrays (%u bytes)", + hist_width * sizeof(int) * 3); + + for (i = 0; i < hist_width; ++i) + hist[i] = 0; + (*histP)[color] = hist; + } else + (*histP)[color] = NULL; +} + + + +static void +clipHistogramAll(unsigned int * const hist[3], + unsigned int const hist_width, + unsigned int const hmax) { + + unsigned int color; + + for (color = 0; color < 3; ++color) + if (hist[color]) + clipHistogram(hist[color], hist_width, hmax); +} + + + +static void +ppm_hist(FILE * const fp, + int const cols, + int const rows, + xelval const maxval, + int const format, + bool const dots, + bool const no_white, + bool const no_black, + bool const colorWanted[3], + bool const verbose, + xelval const startval, + xelval const endval, + unsigned int const hist_width, + unsigned int const hist_height, + bool const clipSpec, + unsigned int const clipCount, + double const hscale) { + + bool const hscale_unity = hscale - 1 < epsilon; + + pixel *pixrow; + pixel **pixels; + int i, j; + unsigned int * hist[3]; /* Subscript is enum wantedColor */ + double vscale; + unsigned int hmax; + + createHist(colorWanted, hist_width, &hist); + + if ((pixels = ppm_allocarray (hist_width, hist_height)) == NULL) + pm_error ("no space for output array (%d pixels)", + hist_width * hist_height); + for (i = 0; i < hist_height; ++i) + memset (pixels[i], 0, hist_width * sizeof (pixel)); + + /* read the pixel values into the histogram arrays */ + pixrow = ppm_allocrow (cols); + /*XX error-check! */ + if (verbose) pm_message ("making histogram..."); + for (i = rows; i > 0; --i) { + ppm_readppmrow (fp, pixrow, cols, maxval, format); + for (j = cols-1; j >= 0; --j) { + int value; + + if (colorWanted[WANT_RED] && + (value = PPM_GETR(pixrow[j])) >= startval && + value <= endval) + hist[WANT_RED][SCALE_H(value-startval)]++; + if (colorWanted[WANT_GRN] && + (value = PPM_GETG(pixrow[j])) >= startval && + value <= endval) + hist[WANT_GRN][SCALE_H(value-startval)]++; + if (colorWanted[WANT_BLU] && + (value = PPM_GETB(pixrow[j])) >= startval && + value <= endval) + hist[WANT_BLU][SCALE_H(value-startval)]++; + } + } + ppm_freerow (pixrow); + fclose (fp); + + /* find the highest-valued slot and set the vertical scale value */ + if (verbose) + pm_message ("finding max. slot height..."); + if (clipSpec) + hmax = clipCount; + else + hmax = maxSlotCountAll(hist, hist_width, no_white, no_black); + + clipHistogramAll(hist, hist_width, hmax); + + vscale = (double) hist_height / hmax; + if (verbose) + pm_message("Done: height = %d, vertical scale factor = %g", + hmax, vscale); + + for (i = 0; i < hist_width; ++i) { + if (hist[WANT_RED]) { + unsigned int j; + bool plotted; + plotted = FALSE; + for (j = hist_height - (int)(vscale * hist[WANT_RED][i]); + j < hist_height && !plotted; + ++j) { + PPM_PUTR(pixels[j][i], maxval); + plotted = dots; + } + } + if (hist[WANT_GRN]) { + unsigned int j; + bool plotted; + plotted = FALSE; + for (j = hist_height - (int)(vscale * hist[WANT_GRN][i]); + j < hist_height && !plotted; + ++j) { + PPM_PUTG(pixels[j][i], maxval); + plotted = dots; + } + } + if (hist[WANT_BLU]) { + unsigned int j; + bool plotted; + plotted = FALSE; + for (j = hist_height - (int)(vscale * hist[WANT_BLU][i]); + j < hist_height && !plotted; + ++j) { + PPM_PUTB(pixels[j][i], maxval); + plotted = dots; + } + } + } + ppm_writeppm (stdout, pixels, hist_width, hist_height, maxval, 0); +} + + + +int +main (int argc, char ** argv) { + + struct cmdlineInfo cmdline; + FILE* ifP; + int cols, rows; + xelval maxval; + int format; + unsigned int hist_width; + unsigned int range; + double hscale; + int hmax; + + pnm_init (&argc, argv); + + parseCommandLine(argc, argv, &cmdline); + + ifP = pm_openr(cmdline.inputFilespec); + + pnm_readpnminit(ifP, &cols, &rows, &maxval, &format); + + range = MIN(maxval, cmdline.rval) - cmdline.lval + 1; + + if (cmdline.widthSpec) + hist_width = cmdline.width; + else + hist_width = range; + + hscale = (float)hist_width / range; + if (hscale - 1.0 < epsilon && cmdline.verbose) + pm_message("Horizontal scale factor: %g (maxval = %u)", + hscale, maxval); + + if (cmdline.nmaxSpec) + hmax = cols * rows / hist_width * cmdline.nmax; + + switch (PNM_FORMAT_TYPE(format)) { + case PPM_TYPE: + ppm_hist(ifP, cols, rows, maxval, format, + cmdline.dots, cmdline.white, cmdline.black, + cmdline.colorWanted, + cmdline.verbose, cmdline.lval, cmdline.rval, + hist_width, cmdline.height, cmdline.nmaxSpec, hmax, hscale); + break; + case PGM_TYPE: + pgm_hist(ifP, cols, rows, maxval, format, + cmdline.dots, cmdline.white, cmdline.black, + cmdline.verbose, cmdline.lval, cmdline.rval, + hist_width, cmdline.height, cmdline.nmaxSpec, hmax, hscale); + break; + case PBM_TYPE: + pm_error("Cannot do a histogram of a a PBM file"); + break; + } + return 0; +} diff --git a/analyzer/pnmpsnr.c b/analyzer/pnmpsnr.c new file mode 100644 index 00000000..4a6bfe56 --- /dev/null +++ b/analyzer/pnmpsnr.c @@ -0,0 +1,209 @@ +/* + * pnmpsnr.c: Compute error (RMSE, PSNR) between images + * + * + * Derived from pnmpnsmr by Ulrich Hafner, part of his fiasco package, + * On 2001.03.04. + + * Copyright (C) 1994-2000 Ullrich Hafner <hafner@bigfoot.de> + */ + +#include <string.h> +#include <stdio.h> +#include <math.h> + +#include "pm_c_util.h" +#include "pam.h" + +#define MAXFILES 16 + +static int +udiff(unsigned int const subtrahend, unsigned int const subtractor) { + return subtrahend-subtractor; +} + + +static double +square(double const arg) { + return(arg*arg); +} + + +static void +validate_input(const struct pam pam1, const struct pam pam2) { + + if (pam1.width != pam2.width) + pm_error("images are not the same width, so can't be compared. " + "The first is %d columns wide, " + "while the second is %d columns wide.", + pam1.width, pam2.width); + if (pam1.height != pam2.height) + pm_error("images are not the same height, so can't be compared. " + "The first is %d rows high, " + "while the second is %d rows high.", + pam1.height, pam2.height); + + if (pam1.maxval != pam2.maxval) + pm_error("images do not have the same maxval. This programs works " + "only on like maxvals. " + "The first image has maxval %u, " + "while the second has %u. Use Pnmdepth to change the " + "maxval of one of them.", + (unsigned int) pam1.maxval, (unsigned int) pam2.maxval); + + if (strcmp(pam1.tuple_type, pam2.tuple_type) != 0) + pm_error("images are not of the same type. The tuple types are " + "'%s' and '%s', respectively.", + pam1.tuple_type, pam2.tuple_type); + + if (strcmp(pam1.tuple_type, PAM_PBM_TUPLETYPE) != 0 && + strcmp(pam1.tuple_type, PAM_PGM_TUPLETYPE) != 0 && + strcmp(pam1.tuple_type, PAM_PPM_TUPLETYPE) != 0) + pm_error("Images are not of a PNM type. Tuple type is '%s'", + pam1.tuple_type); +} + + + +static void +psnr_color(tuple const tuple1, tuple const tuple2, + double * const ySqDiffP, + double * const cbSqDiffP, double * const crSqDiffP) { + + double y1, y2, cb1, cb2, cr1, cr2; + + pnm_YCbCrtuple(tuple1, &y1, &cb1, &cr1); + pnm_YCbCrtuple(tuple2, &y2, &cb2, &cr2); + + *ySqDiffP = square(y1 - y2); + *cbSqDiffP = square(cb1 - cb2); + *crSqDiffP = square(cr1 - cr2); +} + + + +static void +reportPsnr(struct pam const pam1, struct pam const pam2, + double const ySumSqDiff, + double const crSumSqDiff, double const cbSumSqDiff, + const char filespec1[], const char filespec2[]) { + + bool const color = (strcmp(pam1.tuple_type, PAM_PPM_TUPLETYPE) == 0); + + /* The PSNR is the mean of the sum of squares of the differences, + normalized to the range 0..1 + */ + double const yPsnr = ySumSqDiff + / (pam1.width * pam1.height) + / square(pam1.maxval); + + if (color) { + double const cbPsnr = cbSumSqDiff + / (pam1.width * pam1.height) + / square(pam1.maxval); + double const crPsnr = crSumSqDiff + / (pam1.width * pam1.height) + / (pam1.maxval * pam2.maxval); + + pm_message("PSNR between %s and %s:", filespec1, filespec2); + if (yPsnr > 1e-9) + pm_message("Y color component: %.2f dB", 10 * log10(1/yPsnr)); + else + pm_message("Y color component does not differ."); + if (cbPsnr > 1e-9) + pm_message("Cb color component: %.2f dB", 10 * log10(1/cbPsnr)); + else + pm_message("Cb color component does not differ."); + if (crPsnr > 1e-9) + pm_message("Cr color component: %.2f dB", 10 * log10(1/crPsnr)); + else + pm_message("Cr color component does not differ."); + } else { + if (yPsnr > 1e-9) + pm_message("PSNR between %s and %s: %.2f dB", + filespec1, filespec2, 10 * log10(1/yPsnr)); + else + pm_message("Images %s and %s don't differ.", + filespec1, filespec2); + } +} + + + +int +main (int argc, char **argv) { + char *filespec1, *filespec2; /* specs of two files to compare */ + FILE *file1, *file2; + struct pam pam1, pam2; + bool color; + /* It's a color image */ + double ySumSqDiff, crSumSqDiff, cbSumSqDiff; + tuple *tuplerow1, *tuplerow2; /* malloc'ed */ + int row; + + pnm_init(&argc, argv); + + if (argc < 2) + pm_error("Takes two arguments: specifications of the two files."); + else { + filespec1 = argv[1]; + filespec2 = argv[2]; + } + + file1 = pm_openr(filespec1); + file2 = pm_openr(filespec2); + + pnm_readpaminit(file1, &pam1, PAM_STRUCT_SIZE(tuple_type)); + pnm_readpaminit(file2, &pam2, PAM_STRUCT_SIZE(tuple_type)); + + validate_input(pam1, pam2); + + if (strcmp(pam1.tuple_type, PAM_PPM_TUPLETYPE) == 0) + color = TRUE; + else + color = FALSE; + + tuplerow1 = pnm_allocpamrow(&pam1); + tuplerow2 = pnm_allocpamrow(&pam2); + + ySumSqDiff = 0.0; + cbSumSqDiff = 0.0; + crSumSqDiff = 0.0; + + for (row = 0; row < pam1.height; ++row) { + int col; + + pnm_readpamrow(&pam1, tuplerow1); + pnm_readpamrow(&pam2, tuplerow2); + + for (col = 0; col < pam1.width; ++col) { + if (color) { + double ySqDiff, cbSqDiff, crSqDiff; + psnr_color(tuplerow1[col], tuplerow2[col], + &ySqDiff, &cbSqDiff, &crSqDiff); + ySumSqDiff += ySqDiff; + cbSumSqDiff += cbSqDiff; + crSumSqDiff += crSqDiff; + + } else { + unsigned int yDiffSq; + yDiffSq = square(udiff(tuplerow1[col][0], tuplerow2[col][0])); + ySumSqDiff += yDiffSq; + } + } + } + + reportPsnr(pam1, pam2, ySumSqDiff, crSumSqDiff, cbSumSqDiff, + filespec1, filespec2); + + pnm_freepamrow(tuplerow1); + pnm_freepamrow(tuplerow2); + + return 0; +} + + + + + + diff --git a/analyzer/ppmhist.c b/analyzer/ppmhist.c new file mode 100644 index 00000000..4c4d3c55 --- /dev/null +++ b/analyzer/ppmhist.c @@ -0,0 +1,290 @@ +/* ppmhist.c - read a PPM image and compute a color histogram +** +** Copyright (C) 1989 by Jef Poskanzer. +** +** Permission to use, copy, modify, and distribute this software and its +** documentation for any purpose and without fee is hereby granted, provided +** that the above copyright notice appear in all copies and that both that +** copyright notice and this permission notice appear in supporting +** documentation. This software is provided "as is" without express or +** implied warranty. +*/ + +#include <assert.h> + +#include "ppm.h" +#include "shhopt.h" +#include "nstring.h" + +enum sort {SORT_BY_FREQUENCY, SORT_BY_RGB}; + +enum colorFmt {FMT_DECIMAL, FMT_HEX, FMT_FLOAT, FMT_PPMPLAIN}; + +struct cmdline_info { + /* All the information the user supplied in the command line, + in a form easy for the program to use. + */ + const char *input_filespec; /* Filespecs of input files */ + unsigned int noheader; /* -noheader option */ + enum colorFmt colorFmt; + unsigned int colorname; /* -colorname option */ + enum sort sort; /* -sort option */ +}; + + + +static void +parse_command_line(int argc, char ** argv, + struct cmdline_info * const cmdlineP) { +/*---------------------------------------------------------------------------- + Note that the file spec array we return is stored in the storage that + was passed to us as the argv array. +-----------------------------------------------------------------------------*/ + optStruct3 opt; /* set by OPTENT3 */ + optEntry *option_def = malloc(100*sizeof(optEntry)); + unsigned int option_def_index; + + unsigned int hexcolorOpt, floatOpt, mapOpt, nomapOpt; + const char * sort_type; + + option_def_index = 0; /* incremented by OPTENTRY */ + OPTENT3(0, "map", OPT_FLAG, NULL, &mapOpt, 0); + OPTENT3(0, "nomap", OPT_FLAG, NULL, &nomapOpt, 0); + OPTENT3(0, "noheader", OPT_FLAG, NULL, &cmdlineP->noheader, 0); + OPTENT3(0, "hexcolor", OPT_FLAG, NULL, &hexcolorOpt, 0); + OPTENT3(0, "float", OPT_FLAG, NULL, &floatOpt, 0); + OPTENT3(0, "colorname", OPT_FLAG, NULL, &cmdlineP->colorname, 0); + OPTENT3(0, "sort", OPT_STRING, &sort_type, NULL, 0); + + opt.opt_table = option_def; + opt.short_allowed = FALSE; /* We have no short (old-fashioned) options */ + opt.allowNegNum = FALSE; /* We have no parms that are negative numbers */ + + /* Set defaults */ + sort_type = "frequency"; + + optParseOptions3(&argc, argv, opt, sizeof(opt), 0); + /* Uses and sets argc, argv, and some of *cmdlineP and others. */ + + if (argc-1 == 0) + cmdlineP->input_filespec = "-"; + else if (argc-1 != 1) + pm_error("Program takes zero or one argument (filename). You " + "specified %d", argc-1); + else + cmdlineP->input_filespec = argv[1]; + + if (hexcolorOpt + floatOpt + mapOpt > 1) + pm_error("You can specify only one of -hexcolor, -float, and -map"); + if (hexcolorOpt) + cmdlineP->colorFmt = FMT_HEX; + else if (floatOpt) + cmdlineP->colorFmt = FMT_FLOAT; + else if (mapOpt) + cmdlineP->colorFmt = FMT_PPMPLAIN; + else + cmdlineP->colorFmt = FMT_DECIMAL; + + if (strcmp(sort_type, "frequency") == 0) + cmdlineP->sort = SORT_BY_FREQUENCY; + else if (strcmp(sort_type, "rgb") == 0) + cmdlineP->sort = SORT_BY_RGB; + else + pm_error("Invalid -sort value: '%s'. The valid values are " + "'frequency' and 'rgb'.", sort_type); +} + + + +static int +countcompare(const void *ch1, const void *ch2) { + return ((colorhist_vector)ch2)->value - ((colorhist_vector)ch1)->value; +} + + +static int +rgbcompare(const void * arg1, const void * arg2) { + + colorhist_vector const ch1 = (colorhist_vector) arg1; + colorhist_vector const ch2 = (colorhist_vector) arg2; + + int retval; + + retval = (PPM_GETR(ch1->color) - PPM_GETR(ch2->color)); + if (retval == 0) { + retval = (PPM_GETG(ch1->color) - PPM_GETG(ch2->color)); + if (retval == 0) + retval = (PPM_GETB(ch1->color) - PPM_GETB(ch2->color)); + } + return retval; +} + + + +static const char * +colornameLabel(pixel const color, + pixval const maxval, + unsigned int const nDictColor, + pixel const dictColors[], + const char * const dictColornames[]) { +/*---------------------------------------------------------------------------- + Return the name of the color 'color' or the closest color in the + dictionary to it. If the name returned is not the exact color, + prefix it with "*". Otherwise, prefix it with " ". + + 'nDictColor', dictColors[], and dictColorNames[] are the color + dictionary. + + Return the name in static storage within this subroutine. +-----------------------------------------------------------------------------*/ + static char retval[32]; + int colorIndex; + + pixel color255; + /* The color, normalized to a maxval of 255: the maxval of a color + dictionary. + */ + + PPM_DEPTH(color255, color, maxval, 255); + + colorIndex = ppm_findclosestcolor(dictColors, nDictColor, &color); + + assert(colorIndex >= 0 && colorIndex < nDictColor); + + if (PPM_EQUAL(dictColors[colorIndex], color)) + STRSCPY(retval, " "); + else + STRSCPY(retval, "*"); + + STRSCAT(retval, dictColornames[colorIndex]); + + return retval; +} + + + +static void +printColors(colorhist_vector const chv, + int const nColors, + pixval const maxval, + enum colorFmt const colorFmt, + unsigned int const nKnown, + pixel const knownColors[], + const char * const colornames[]) { + + int i; + + for (i = 0; i < nColors; i++) { + pixval const r = PPM_GETR(chv[i].color); + pixval const g = PPM_GETG(chv[i].color); + pixval const b = PPM_GETB(chv[i].color); + double const lum = PPM_LUMIN(chv[i].color); + unsigned int const intLum = lum + 0.5; + double const floatLum = lum / maxval; + unsigned int const count = chv[i].value; + + const char * colornameValue; + + if (colornames) + colornameValue = colornameLabel(chv[i].color, maxval, + nKnown, knownColors, colornames); + else + colornameValue = ""; + + switch(colorFmt) { + case FMT_FLOAT: + printf(" %1.3f %1.3f %1.3f\t%1.3f\t%7d %s\n", + (double)r / maxval, + (double)g / maxval, + (double)b / maxval, + floatLum, count, colornameValue); + break; + case FMT_HEX: + printf(" %04x %04x %04x\t%5d\t%7d %s\n", + r, g, b, intLum, count, colornameValue); + break; + case FMT_DECIMAL: + printf(" %5d %5d %5d\t%5d\t%7d %s\n", + r, g, b, intLum, count, colornameValue); + break; + case FMT_PPMPLAIN: + printf(" %5d %5d %5d#\t%5d\t%7d %s\n", + r, g, b, intLum, count, colornameValue); + break; + } + } +} + + + +int +main(int argc, char *argv[] ) { + struct cmdline_info cmdline; + FILE* ifP; + colorhist_vector chv; + int rows, cols; + pixval maxval; + int format; + int nColors; + int (*compare_function)(const void *, const void *); + /* The compare function to be used with qsort() to sort the + histogram for output + */ + unsigned int nDictColor; + const char ** dictColornames; + pixel * dictColors; + + ppm_init( &argc, argv ); + + parse_command_line(argc, argv, &cmdline); + + ifP = pm_openr(cmdline.input_filespec); + + ppm_readppminit(ifP, &cols, &rows, &maxval, &format); + + chv = ppm_computecolorhist2(ifP, cols, rows, maxval, format, 0, &nColors); + + pm_close(ifP); + + switch (cmdline.sort) { + case SORT_BY_FREQUENCY: + compare_function = countcompare; break; + case SORT_BY_RGB: + compare_function = rgbcompare; break; + } + + qsort((char*) chv, nColors, sizeof(struct colorhist_item), + compare_function); + + /* And print the histogram. */ + if (cmdline.colorFmt == FMT_PPMPLAIN) + printf("P3\n# color map\n%d 1\n%d\n", nColors, maxval); + + if (!cmdline.noheader) { + const char commentDelim = cmdline.colorFmt == FMT_PPMPLAIN ? '#' : ' '; + printf("%c r g b \t lum \t count %s\n", + commentDelim, cmdline.colorname ? "name" : ""); + printf("%c----- ----- ----- \t-----\t------- %s\n", + commentDelim, cmdline.colorname ? "----" : ""); + } + if (cmdline.colorname) { + bool mustOpenTrue = TRUE; + ppm_readcolordict(NULL, mustOpenTrue, + &nDictColor, &dictColornames, &dictColors, NULL); + } else { + dictColors = NULL; + dictColornames = NULL; + } + + printColors(chv, nColors, maxval, + cmdline.colorFmt, nDictColor, dictColors, dictColornames); + + if (dictColors) + free(dictColors); + if (dictColornames) + free(dictColornames); + + ppm_freecolorhist(chv); + + return 0; +} |