CDMS/src/CDMS/Partridge/FieldFlattening
diff -N BlankImageFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ BlankImageFinder.java 4 Oct 2010 22:45:30 -0000 1.1
@@ -0,0 +1,101 @@
+/**
+ * BlankImageFinder.java
+ */
+
+package CDMS.Partridge.FieldFlattening;
+
+import CDMS.Partridge.ImageUtils.ImageData;
+import CDMS.Partridge.ImageUtils.ImageDescriptor;
+import CDMS.Partridge.ImageUtils.PixelArray;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Map.Entry;
+
+/**
+ * Find blank images by identifying images with small standard deviations
+ * for the pixels in the image.
+ *
+ * @author Richard Partridge
+ */
+public class BlankImageFinder {
+
+ public static void main(String[] args) throws IOException {
+
+ // Define the inputs and outputs - eventually it would be good to package
+ // this application and invoke these parameters from the command line
+ String fbase = "c:/CDMS/IZipG47_Side1_PKG/IZipG47_Side1_PKG-";
+ int ifirst = 1;
+ int ilast = 8585;
+ String suffix = ".png";
+ String outfile = "c:/CDMS/IZipv4_side1_blanks.txt";
+ double imin = 100.;
+ int nblank = 100;
+
+ // Create maps to hold the mean and sd calculations
+ Map<String, Double> imean = new HashMap<String, Double>();
+ Map<String, Double> isd = new HashMap<String, Double>();
+
+ // Loop over the image files
+ for (int ifile=ifirst; ifile<=ilast; ifile++) {
+
+ // Show progress in processing the image files
+ if (ifile > 0 && (ifile % 100) == 0) System.out.println("Processing file "+ifile);
+
+ // Open the image file with arbitrary origin and pixel size
+ String fname = fbase + ifile + suffix;
+ ImageDescriptor ID = new ImageDescriptor(fname, ifile, 0., 0., 1., 1.);
+ ImageData image = ID.getImageData();
+ PixelArray pixarray = ID.getPixelArray();
+
+ // Loop over pixels and find the intensity mean and standard deviation
+ double itot = 0.;
+ double i2tot = 0.;
+ for (int row=0; row<ID.getNRow(); row++) {
+ for (int col=0; col<ID.getNCol(); col++) {
+ double intensity = image.getIntensity(row, col);
+ itot += intensity;
+ i2tot += intensity*intensity;
+ }
+ }
+ int npixel = ID.getNRow() * ID.getNCol();
+ double mean = itot / npixel;
+ double sd = Math.sqrt(i2tot/npixel - mean*mean);
+
+ // Save the mean and st dev in maps for later access
+ imean.put(fname, mean);
+ isd.put(fname, sd);
+ }
+
+ // Make a list of sd for blank candidates, excluding those that have
+ // low mean intensity (typically images outside the detector boundary)
+ List<Double> sdlist = new ArrayList<Double>();
+ for (Entry<String, Double> ent : imean.entrySet()) {
+ if (ent.getValue() > imin) sdlist.add(isd.get(ent.getKey()));
+ }
+
+ // Sort the list and find the cutoff sd for our list of blank images
+ Collections.sort(sdlist);
+ double sdmax = sdlist.get(nblank-1);
+
+ // Open the output file
+ PrintWriter out = new PrintWriter(new FileWriter(outfile));
+
+ // Output the desired number of images choosing those with the smallest sd
+ for (String fname : imean.keySet()) {
+ double mean = imean.get(fname);
+ double sd = isd.get(fname);
+ if (sd > sdmax || mean <= imin) continue;
+ System.out.println("File "+fname+" selected - mean: "+mean+" sd: "+sd);
+ out.println(fname);
+ }
+
+ // Close the output file
+ out.close();
+ }
+}
CDMS/src/CDMS/Partridge/FieldFlattening
diff -N FieldFlattener.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FieldFlattener.java 4 Oct 2010 22:45:30 -0000 1.1
@@ -0,0 +1,422 @@
+/**
+ * FieldFlattener.java
+ */
+
+package CDMS.Partridge.FieldFlattening;
+
+import CDMS.Partridge.ImageUtils.ImageData;
+import CDMS.Partridge.ImageUtils.ImageDescriptor;
+import CDMS.Partridge.ImageUtils.PixelArray;
+import java.io.FileNotFoundException;
+import java.io.FileReader;
+import java.io.StreamTokenizer;
+import java.util.ArrayList;
+import java.util.List;
+
+/**
+ * Takes a set of uniform blank images and creates a 2D field flattening map that
+ * eliminates variations in intensity among the different image pixels. These
+ * variations may arise from non-uniform illumination and/or camera response.
+ *
+ * The BlankImageFinder class can be used to identify blank images and create
+ * the blank image file used by this class.
+ *
+ * The FlattenImage method is of particular note - it applies the field
+ * flattening map to an image and returns a new image with field flattening.
+ *
+ * All blank images, as well as images processed for field flattening, must be
+ * identical in size.
+ *
+ * @author Richard Partridge
+ */
+public class FieldFlattener {
+
+ private List<String> _fnames;
+ private int _nrow;
+ private int _ncol;
+ private double[][] _rcor;
+ private double[][] _gcor;
+ private double[][] _bcor;
+ private int _max_iterations = 10;
+
+ /**
+ * Fully qualified constructor
+ *
+ * @param blankfile file containing the list of blank images
+ * @param cut number of std dev used for the outlier removal cut
+ */
+ public FieldFlattener(String blankfile, double cut) {
+
+ // Read in the list of blank images
+ ReadBlankFile(blankfile);
+
+ // Initialize the arrays for the image data and pixel descriptors
+ int nfiles = _fnames.size();
+ if (nfiles == 0) throw new RuntimeException("No blank images found");
+ ImageData[] _images = new ImageData[nfiles];
+ PixelArray[] _pixarray = new PixelArray[nfiles];
+
+ // Loop over files and retrieve the image data
+ for (int ifile=0; ifile<nfiles; ifile++) {
+
+ // Create an image descriptor with arbitrary position / pixel size
+ String fname = _fnames.get(ifile);
+ ImageDescriptor ID = new ImageDescriptor(fname, ifile, 0., 0., 1., 1.);
+
+ // Get the image data and pixel array descriptor for this file
+ _images[ifile] = ID.getImageData();
+ _pixarray[ifile] = ID.getPixelArray();
+
+ // If this is the first image, save the image size
+ if (ifile == 0) {
+ _nrow = _pixarray[0].getNRow();
+ _ncol = _pixarray[0].getNCol();
+ }
+
+ // Check that all images are the same size
+ if (_pixarray[ifile].getNRow() != _nrow ||
+ _pixarray[ifile].getNCol() != _ncol)
+ throw new RuntimeException("Blank image files have different sizes");
+ }
+
+ // Create arrays for the average color intensities of each pixel
+ double ared[][] = new double[_nrow][_ncol];
+ double agreen[][] = new double[_nrow][_ncol];
+ double ablue[][] = new double[_nrow][_ncol];
+
+ // Loop over all pixels
+ for (int row=0; row<_nrow; row++) {
+ for (int col=0; col<_ncol; col++) {
+
+ // Arrays for the color intensities from each blank image
+ int[] red = new int[nfiles];
+ int[] blue = new int[nfiles];
+ int[] green = new int[nfiles];
+
+ // Loop over the blank images
+ for (int ifile=0; ifile<nfiles; ifile++) {
+
+ // Save the color intensities found for this pixel
+ ImageData image = _images[ifile];
+ int rgb = image.getRGB(row, col);
+ red[ifile] = image.getRed(rgb);
+ blue[ifile] = image.getBlue(rgb);
+ green[ifile] = image.getGreen(rgb);
+ }
+
+ // Average the color intensities for this pixel using a
+ // truncated mean that discards outliers since our blank
+ // images may not be flawless
+ ared[row][col] = TruncatedMean(red, cut);
+ ablue[row][col] = TruncatedMean(blue, cut);
+ agreen[row][col] = TruncatedMean(green, cut);
+ }
+ }
+
+ // Find the global average values of each color intensity averaging over
+ // the entire image the color intensities that were averaged
+ double rave = 0.;
+ double gave = 0.;
+ double bave = 0.;
+ for (int row=0; row<_nrow; row++) {
+ for (int col=0; col<_ncol; col++) {
+ rave += ared[row][col] / (_nrow * _ncol);
+ gave += agreen[row][col] / (_nrow * _ncol);
+ bave += ablue[row][col] / (_nrow * _ncol);
+ }
+ }
+
+ // Create the arrays to hold the field flattening map
+ _rcor = new double[_nrow][_ncol];
+ _gcor = new double[_nrow][_ncol];
+ _bcor = new double[_nrow][_ncol];
+
+ // Calculate the correction factors to be applied to the
+ // color intensities for each image pixel
+ for (int row=0; row<_nrow; row++) {
+ for (int col=0; col<_ncol; col++) {
+
+ _rcor[row][col] = rave / ared[row][col];
+ _gcor[row][col] = gave / agreen[row][col];
+ _bcor[row][col] = bave / ablue[row][col];
+ }
+ }
+ }
+
+ /**
+ * Return the array of field flattening correction factors for
+ * the red color intensities.
+ *
+ * flatRed = RedCorrection[row][col] * originalRed
+ *
+ * @return correction factor indexed by [row][col]
+ */
+ public double[][] getRedCorrection() {
+ return _rcor;
+ }
+
+ /**
+ * Return the array of field flattening correction factors for
+ * the green color intensities.
+ *
+ * flatGreen = GreenCorrection[row][col] * originalGreen
+ *
+ * @return correction factor indexed by [row][col]
+ */
+ public double[][] getGreenCorrection() {
+ return _gcor;
+ }
+
+ /**
+ * Return the array of field flattening correction factors for
+ * the blue color intensities.
+ *
+ * flatBlue = BlueCorrection[row][col] * originalBlue
+ *
+ * @return correction factor indexed by [row][col]
+ */
+ public double[][] getBlueCorrection() {
+ return _bcor;
+ }
+
+ /**
+ * Return the list of blank image files used to construct the field
+ * flattening correction.
+ *
+ * @return
+ */
+ public List<String> getBlankFileNames() {
+ return _fnames;
+ }
+
+ /**
+ * Apply the field flattening corrections to an image, creating a new
+ * "flattened" image.
+ *
+ * @param image image to be flattened
+ * @return flattened image
+ */
+ public ImageData FlattenImage(ImageData image) {
+
+ // Get the pixel array descriptor and create an identical new empty image
+ PixelArray pixarray = image.getPixelArray();
+ ImageData flatimage = new ImageData(pixarray);
+
+ // Check that the dimensions of the image to be flattened are the
+ // same as the blank images used to develop the flattening correction
+ if (pixarray.getNRow() != _nrow || pixarray.getNCol() != _ncol)
+ throw new RuntimeException("Image sizes don't match");
+
+ // Loop over the image pixels
+ for (int row=0; row<_nrow; row++) {
+ for (int col=0; col<_ncol; col++) {
+
+ // Get the packed rgb color intensities for this pixel
+ int rgb = image.getRGB(row, col);
+
+ // Unpack the color intensities. apply the flattening
+ // correction, and round off to the nearest integer
+ int red = (int) Math.round(image.getRed(rgb) * _rcor[row][col]);
+ int green = (int) Math.round(image.getGreen(rgb) * _gcor[row][col]);
+ int blue = (int) Math.round(image.getBlue(rgb) * _bcor[row][col]);
+
+ // Save the corrected color intensities in the flattened image
+ flatimage.setRGB(row, col, red, green, blue);
+ }
+ }
+ return flatimage;
+ }
+
+ /**
+ * Read and parse the blank image file.
+ *
+ * @param blankfile file containing the list of blank images
+ */
+ private void ReadBlankFile(String blankfile) {
+
+ // Initialize the list of file names
+ _fnames = new ArrayList<String>();
+
+ // Open the file with the list of blank images
+ FileReader fr;
+ try {
+ fr = new FileReader(blankfile);
+ } catch (FileNotFoundException ex) {
+ throw new RuntimeException("Blank image index file not found: "+blankfile);
+ }
+
+ // Setup the stream tokenizer for decoding the file
+ StreamTokenizer tok = new StreamTokenizer(fr);
+ tok.resetSyntax();
+ tok.wordChars(32, 126);
+ tok.parseNumbers();
+ tok.whitespaceChars(0, 32);
+ tok.eolIsSignificant(true);
+
+ // Loop over the blank file entries
+ try {
+ while (tok.nextToken() != StreamTokenizer.TT_EOF) {
+ // Get the next file name
+ if (tok.ttype != StreamTokenizer.TT_WORD) ParseError();
+ String fname = tok.sval;
+ _fnames.add(fname);
+ if (tok.nextToken() != StreamTokenizer.TT_EOL) ParseError();
+ }
+ } catch (Exception ex) {
+ throw new RuntimeException("Empty file IO error");
+ }
+ }
+
+ /**
+ * Throw an exception if there is a parsing error for the blank image list.
+ */
+ private void ParseError() {
+ throw new RuntimeException("Invalid blank file format");
+ }
+
+ /**
+ * Calculated the truncated mean from an array of values. The algorithm
+ * is an iterative one, where initially all values are included in the
+ * calculation of the mean and stddev, but in subsequent iterations those
+ * values that are more than cut * stddev away from the last mean are not
+ * included in the calculation of an updated mean and stddev. This procedure
+ * is repeated until either the result is stable or a maximum number of
+ * iterations is reached (to handle cases where the result ping-pongs
+ * between two or more values).
+ *
+ * @param values array of values to be averaged
+ * @param cut truncation cut in units of the standard deviation
+ * @return truncated mean
+ */
+ private double TruncatedMean(int[] values, double cut) {
+
+ // Get the number of values to average
+ int nval = values.length;
+
+ // Create and initialize the array of values to be excluded
+ boolean[] exclude = new boolean[nval];
+ for (int i=0; i<nval; i++) {
+ exclude[i] = false;
+ }
+
+ // Calculate the initial mean and stddev using all values
+ double mean = CalcMean(values, exclude);
+ double sd = CalcStdDev(values, mean, exclude);
+
+ // Iterate the calculation of the truncated mean for up to _max_iterations
+ for (int i=0; i<_max_iterations; i++) {
+
+ // Update the array of excluded (truncated) values, stop looping
+ // when the result is stable (i.e., no changes to exclusion array)
+ if (!UpdateExclusions(mean, sd, cut, values, exclude)) break;
+
+ // Calculate an updated mean and standard deviation using
+ // the updated exclusion array
+ mean = CalcMean(values, exclude);
+ sd = CalcStdDev(values, mean, exclude);
+ }
+
+ // Finished - return the last calculation of the truncated mean
+ return mean;
+ }
+
+ /**
+ * Calculate the mean of an array of values, excluding those values flagged
+ * for exclusion.
+ *
+ * @param values array of values to be averaged
+ * @param exclude values to be excluded are flagged true
+ * @return mean
+ */
+ private double CalcMean(int[] values, boolean[]exclude) {
+
+ // Initialize
+ int nval = values.length;
+ int nsum = 0;
+ double sum = 0.;
+
+ // Count and sum the values not excluded
+ for (int i=0; i<nval; i++) {
+ if (exclude[i]) continue;
+ nsum++;
+ sum += values[i];
+ }
+
+ // Return the mean of the values not excluded
+ return sum / nsum;
+ }
+
+ /**
+ * Calculate the standard deviation of an array of values, excluding those
+ * values flagged for exclusion.
+ *
+ * @param values array of values
+ * @param mean mean value of the values not marked for exclusion
+ * @param exclude values to be excluded are flagged true
+ * @return standard deviation
+ */
+ private double CalcStdDev(int[] values, double mean, boolean[]exclude) {
+
+ // Initialize
+ int nval = values.length;
+ int nsum = 0;
+ double sum = 0.;
+
+ // Count and sum the squared deviation from the mean for values not excluded
+ for (int i=0; i<nval; i++) {
+ if (exclude[i]) continue;
+ nsum++;
+ double dif = values[i] - mean;
+ sum += dif*dif;
+ }
+
+ // Calculate the standard deviation estimator for values not excluded
+ double sd;
+ if (nsum <= 2) sd = Math.sqrt(sum);
+ else sd = Math.sqrt(sum / (nsum - 1));
+
+ return sd;
+ }
+
+ /**
+ * Update the list of values to exclude in the truncated mean calculation.
+ * On input, the array of excluded values should be those from the previous
+ * iteration upon which the mean and standard deviation inputs are calculated.
+ * On return from this method, the array of excluded values will be
+ * updated to identify which values should be excluded in the next iteration.
+ *
+ * @param mean average of non-excluded values from previous iteration
+ * @param sd standard deviation of non-excluded values from previous iteration
+ * @param cut truncation cut in units of standard deviations
+ * @param values array of values
+ * @param exclude values to be excluded are flagged true
+ * @return
+ */
+ private boolean UpdateExclusions(double mean, double sd, double cut,
+ int[] values, boolean[] exclude) {
+
+ // Initialize - don't bother with exclusions if there is only 1 value
+ boolean changed = false;
+ int nval = values.length;
+ if (nval < 2) return false;
+
+ // Calculate the minimum and maximum allowed values - make sure our
+ // allowed window is at least +- 1 bin wide and at least +- 1 sd
+ // wide to ensure we don't exclude everything
+ double delta = sd * cut;
+ if (cut < 1.0) delta = sd;
+ if (delta < 1.) delta = 1.;
+ int vmin = (int) Math.floor(mean - sd * cut);
+ int vmax = (int) Math.ceil(mean + sd * cut);
+
+ // Check if any exclusions have changed and mark new exclusions
+ for (int i=0; i<nval; i++) {
+ boolean bad = values[i] < vmin || values[i] > vmax;
+ if (bad != exclude[i]) changed = true;
+ exclude[i] = bad;
+ }
+
+ // Return true if there were changes in the exclusion array
+ return changed;
+ }
+}
\ No newline at end of file
CDMS/src/CDMS/Partridge/FieldFlattening
diff -N FieldFlattenerTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FieldFlattenerTest.java 4 Oct 2010 22:45:30 -0000 1.1
@@ -0,0 +1,54 @@
+/*
+ * FieldFlattenerTest.java
+ */
+
+package CDMS.Partridge.FieldFlattening;
+
+import CDMS.Partridge.ImageUtils.ImageData;
+import CDMS.Partridge.ImageUtils.ImageDescriptor;
+import CDMS.Partridge.ImageUtils.PixelArray;
+import java.io.File;
+import java.io.IOException;
+import javax.imageio.ImageIO;
+
+/**
+ *
+ * @author Richard Partridge
+ */
+public class FieldFlattenerTest {
+
+ public static void main(String[] args) throws IOException {
+
+ String blankindexfile = "c:/CDMS/IZipv4_side1_blanks.txt";
+ FieldFlattener flat = new FieldFlattener(blankindexfile, 2.0);
+ int ifile = 0;
+ for (String fname : flat.getBlankFileNames()) {
+ ImageDescriptor ID = new ImageDescriptor(fname, 1, 0., 0., 1., 1.);
+ ImageData image = ID.getImageData();
+ PixelArray pixarray = ID.getPixelArray();
+ ImageData flatimage = flat.FlattenImage(image);
+ double isum = 0.;
+ double i2sum = 0.;
+ int nrow = ID.getNRow();
+ int ncol = ID.getNCol();
+ for (int row=0; row<nrow; row++) {
+ for (int col=0; col<ncol; col++) {
+ double intensity = flatimage.getIntensity(row, col);
+ isum += intensity;
+ i2sum += intensity*intensity;
+ }
+ }
+ int npixel = nrow * ncol;
+ double mean = isum / npixel;
+ double sd = Math.sqrt(i2sum/npixel - mean*mean);
+ System.out.println("Flattened file: "+fname+" mean: "+mean+" sd: "+sd);
+ String otype = "png";
+ File outfile = new File("c:/CDMS/flattened"+ifile+"."+otype);
+ ifile++;
+ ImageIO.write(flatimage.getBufferedImage(), otype, outfile);
+
+ }
+
+ }
+
+}
