broken resolution enhancement
Kragen Sitaker
kragen@pobox.com
Fri, 31 Mar 2000 22:54:27 -0500
Well, actually, it sort of does work. It just doesn't work for what it
was written for.
It seems to have some overflow problems. When modified to work on real
Landsat data (by pretending the last pixel on each row doesn't exist,
and just writing out uninitialized data for it), it ends up pushing the
limits --- whether due to an undetected bug or just the data, I'm not
sure.
/*
* resolution enhancement.
* Landsat 7 ETM+, as well as other remote sensors, has a high-resolution
* panchromatic band and several lower-resolution multispectral bands. Here's
* a straightforward algorithm I cooked up to use the pan band to
* resolution-enhance a multi band.
*
* Here is the problem: how can I use the pan band in such a way that
* averaging the pixels of the res-enhanced multi band will yield the original
* multi band? Keeping the problem simple is that all multispectral pixel
* boundaries coincide exactly with pan pixel boundaries. (In the Landsat-7
* case, each multispectral pixel corresponds to exactly four panchromatic
* pixels.) [Yow, these are both frighteningly wrong; see below for details.]
*
* The answer was not immediately apparent, but I think it's the obvious
* solution if you think about it for a while. Compute the average of the pan
* pixels over the area of one multi pixel; divide the multi pixel value by
* that pan average value. Then, normalize the pan pixels by multiplying them
* by the quotient. The normalized values are the res-enhanced multi pixels.
* You can see that averaging them will give the original multi value.
*
* This works on the assumption that local variations in brightness between the
* pan pixels do a good job of predicting local variations in brightness in the
* multispectral band. We can validate this assumption by doing a linear
* regression between the multi pixel values and the pan average values. If
* the correlation is low, then the pan band does a bad job of predicting the
* brightness of multispectral pixels, and thus, presumably, a bad job of
* predicting variations in their brightnesses. If the correlation is high,
* then the pan band does a good job of predicting the brightness of
* multispectral pixels.
*
* This algorithm has the notable advantage of leaving the spectral content of
* the pixels essentially unchanged, except by linear scaling.
*
* It is probably not difficult to improve on this algorithm.
*
* This implementation assumes that the data is, like Landsat-7 FAST-L7A
* format, in raw binary files stored with each row contiguous. The pixel data
* type is defined in the pix typedef.
*
* This implementation doesn't do the linear regression.
*
* To the best of my knowledge, this program is free of any copyright claim.
* I, Kragen Sitaker, wrote it on 2000-03-31; I am an employee of Veridian ERIM
* International, which is the reason I have Landsat 7 data to play with, but
* resolution enhancement of Landsat 7 data is not part of my job. I just want
* a nice Landsat poster to hang on my wall.
*
* I hereby relinquish any copyright I may hold in it; since I am the only
* holder of copyright in it, that puts it in the public domain.
*/
/* hmm, there's a problem. It turns out it's not the edges of pixels that
* coincide; it's their centers. And therein lies the rub; pan width is not
* twice multi width, but (pan width - 1) is twice (multi width - 1).
*
* The cheap-ass thing to do is to do all but the last row and column; this
* turns out to be possible, but painful, and wrong.
*
* If I were rich, I could restore the invariants by doubling everything on
* disk. Unfortunately, that would require 800 megabytes for each band, and I
* don't have *that* much disk space. It would introduce funny morphological
* features, too.
*
* The right thing to do, following the paradigm of the original algorithm, is
* to find the normalized values from both of the multispectral pixels that
* overlap a given panchromatic pixel. Then average them.
*
* All in all, this is a pain in the ass. It removes some of the original
* algorithm's nicer properties. It requires that we buffer up two lines of
* multispectral data.
*
* So I'm sending out this broken program, which doesn't actually work on real
* Landsat data.
*
* It appears to work on some test data: the 2x2 "multispectral
* band" "ABCD" with the 4x4 "panchromatic band" "5455655554546565" produces
* the "resolution-enhanced band" "A@BBBABBCBDCDCED", which is fine. A@BA
* matches 5465 and averages to A; BBBB matches 5555 and averages to B; CBDC
* matches 5465 and averages to C; and DCED matches 5465 and averages to D.
* This may make more sense if you're looking at an ASCII chart and thinking
* about the numerical values of these characters.
*/
#include <stdio.h>
#include <stdlib.h>
/* we'll buffer one line of multispectral pixels and one multispectral line's
* worth of panchromatic pixels. */
typedef unsigned char pix;
/* no globals in this function, just too many args */
void compute_avgs(pix *multi_in, double *pan_factor, pix *pan, int multiwidth,
int width_ratio, int height_ratio) {
int i, j, k;
for (i = 0; i != multiwidth; i++) {
double sum = 0.0;
for (j = 0; j != height_ratio; j++) {
for (k = 0; k != width_ratio; k++) {
sum += pan[
i * width_ratio + k +
j * multiwidth * width_ratio
];
}
}
pan_factor[i] = multi_in[i]
/
(sum / height_ratio / width_ratio);
/* this would be where you would do linear regression */
}
}
/* again, no globals, but too damn many args. Assumes either truncation toward
* negative infinity or unsigned input values; the latter is true of Landsat.
*/
void compute_row(double *pan_factor, pix *pan, pix *multi_out, int multiwidth,
int width_ratio, int panrow) {
int i, j;
for (i = 0; i != multiwidth; i++) {
for (j = 0; j != width_ratio; j++) {
int offset = i * width_ratio + j;
multi_out[offset] = (int)
(pan_factor[i] *
pan[offset + panrow * multiwidth * width_ratio]
+ 0.5);
}
}
}
pix *multi_in, *pan, *multi_out;
double *pan_factor;
void usage(char *argv0) {
fprintf(stderr, "Usage: %s multiwidth multiheight panwidth "
"panheight multi_in_file \n"
"\tpan_in_file multi_out_file\n", argv0);
exit(2);
}
/* this function is too big and has too many variables */
int main(int argc, char **argv) {
int multiwidth, multiheight, panwidth, panheight;
int width_ratio, height_ratio;
int row;
FILE *multi_in_file, *pan_file, *multi_out_file;
if (argc != 8) {
usage(argv[0]);
}
#define num_arg(name, index) \
if (!(name = atoi(argv[index]))) { \
fprintf(stderr, "%s: bad arg %s '%s' (should be numeric and nonzero)\n", \
argv[0], #name, argv[index]); \
usage(argv[0]); \
}
num_arg(multiwidth, 1)
num_arg(multiheight, 2)
num_arg(panwidth, 3)
num_arg(panheight, 4)
#define filename_arg(name, index, mode) \
if (!(name = fopen(argv[index], mode))) { \
fprintf(stderr, "%s: %s: ", argv[0], #name); \
perror(argv[index]); \
usage(argv[0]); \
}
filename_arg(multi_in_file, 5, "rb")
filename_arg(pan_file, 6, "rb")
filename_arg(multi_out_file, 7, "wb")
/* note assumption of coincident pixel boundaries */
width_ratio = panwidth / multiwidth;
height_ratio = panheight / multiheight;
if (panwidth != width_ratio * multiwidth) {
fprintf(stderr, "%s: %d is not a multiple of %d, "
"donc Dieu existe --- repondez!\n",
/* --- Euler */
argv[0], panwidth, multiwidth);
return 1;
}
if (panheight != height_ratio * multiheight) {
fprintf(stderr, "%s: %d is not a multiple of %d, "
"donc Dieu existe --- repondez!\n",
argv[0], panheight, multiheight);
return 1;
}
if (!width_ratio || !height_ratio) {
fprintf(stderr, "%s: panchromatic band must be at least "
"as large in both dimensions\n as multispectral band\n",
argv[0]);
return 1;
}
multi_in = malloc(multiwidth * sizeof(multi_in[0]));
pan_factor = malloc(multiwidth * sizeof(pan_factor[0]));
pan = malloc(height_ratio * panwidth * sizeof(pan[0]));
multi_out = malloc(panwidth * sizeof(multi_out[0]));
if (!(multi_in && pan_factor && pan && multi_out)) {
fprintf(stderr, "%s: out of memory", argv[0]);
return 2;
}
for (row = 0; row != multiheight; row++) {
int panrow;
int size;
size = multiwidth * sizeof(multi_in[0]);
if (!fread(multi_in, size, 1, multi_in_file)) {
fprintf(stderr, "Error reading %d bytes from %s:"
" EOF at row %d\n",
size, argv[5], row);
return 1;
}
size = height_ratio * panwidth * sizeof(pan[0]);
if (!fread(pan, size, 1, pan_file)) {
fprintf(stderr, "Error reading %d bytes from %s:"
" EOF around row %d\n",
size, argv[6], row * height_ratio);
return 1;
}
compute_avgs(multi_in, pan_factor, pan,
multiwidth, width_ratio, height_ratio);
for (panrow = 0; panrow != height_ratio; panrow++) {
compute_row(pan_factor, pan, multi_out, multiwidth,
width_ratio, panrow);
if (!fwrite(multi_out, panwidth * sizeof(multi_out[0]),
1, multi_out_file)) {
fprintf(stderr, "Error writing ");
perror(argv[7]);
return 1;
}
}
}
if (fread(multi_in, 1, 1, multi_in_file)) {
fprintf(stderr, "Warning: %s too big (trailing ignored)\n",
argv[5]);
}
if (fread(multi_in, 1, 1, pan_file)) {
fprintf(stderr, "Warning: %s too big (trailing ignored)\n",
argv[6]);
}
fclose(multi_out_file);
free(multi_in);
free(pan_factor);
free(pan);
free(multi_out);
return 0;
}