calc_klt.c

JPEG2000实现的源码

/*****************************************************************************/
/* Copyright 2000, Eastman Kodak Company                                     */
/* All rights reserved                                                       */
/* File: "calc_klt.c"                                                        */
/* Description: Main driver code for KLT calculation.                        */
/* Author: Austin Lan                                                        */
/* Affiliation: Eastman Kodak Company                                        */
/* Version: VM8.5                                                            */
/* Last Revised: 11 September, 2000                                          */
/*****************************************************************************/
#include <local_services.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <ifc.h>
#include <time.h>
#include <image_io.h>

#include "meschach_utils.h" /* For eigenvector calculation */

/* ========================================================================= */
/* -------------------------- Internal Functions --------------------------- */
/* ========================================================================= */

#define COPY_VOID_TO_FLOAT(dest, src, type, rows, cols) {\
   int elem = (rows * cols);\
   float *dest_p = dest[0];\
   type *src_p = ((type **)src)[0];\
   while (elem-- > 0) *dest_p++ = (float)*src_p++;\
}

/*****************************************************************************/
/* STATIC                      print_generic_usage                           */
/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* Function:  print_generic_usage                                            */
/* Author: David Taubman                                                     */
/* Affiliation: Hewlett-Packard and                                          */
/*              The University of New South Wales, Australia                 */
/* Version: VM5.2                                                            */
/* Last Revised: 10 September, 1999                                          */
/*****************************************************************************/

static void
  print_generic_usage(char **usage)

 /* Formats and prints the usage information recovered from the `get_usage'
    function offered by any of the standard objects. */

{
  if (usage == NULL)
    return;
  for (; *usage != NULL; usage += 2)
    {
      local_printf(4,70,"%s",usage[0]);
      if (usage[1] != NULL)
        local_printf(7,70,"%s",usage[1]);
    }
}

/*****************************************************************************/
/* STATIC                      construct_reader                              */
/*****************************************************************************/

static image_reader_ref construct_reader(cmdl_ref cmdl)
{
   int k, num_args, have_pgx, have_mcraw;
   char **params;
   extern image_reader_ref create_mcraw_image_reader(void);
   extern image_reader_ref create_simple_image_reader(void);
   image_reader_ref reader;
   
   if ((num_args = cmdl->extract(cmdl,"-i",-1,&params)) <= 0) {
      return NULL;
   }

   have_pgx = have_mcraw = 0;
   for (k=0; k < num_args; k++) {
      if ((strstr(params[k], ".img") != NULL) ||
	  (strstr(params[k], ".nrm") != NULL) ||
	  (strstr(params[k], ".bil") != NULL) ||
	  (strstr(params[k], ".raw") != NULL) ||
	  (strstr(params[k], ".bsq") != NULL) ||
	  (strstr(params[k], ".bip") != NULL))
	 have_mcraw = 1;
      else
	 have_pgx = 1;
   }
   if ((have_pgx + have_mcraw) > 1) {
      local_error("Unable to support different input file formats on the "
		  "same command line...");
   }

   if (have_mcraw) {
      reader = create_mcraw_image_reader();
   } else {
      reader = create_simple_image_reader();
   }

   return reader;
}

/*****************************************************************************/
/* STATIC                        allocate_1d                                 */
/*****************************************************************************/

static void allocate_1d(void **ptr, int dim1, int byte_depth)
{
   void *mem_p;

   if ((mem_p = (void *)local_malloc("KLT", (dim1 * byte_depth))) == NULL) {
      local_error("allocate_1d():  calloc failure.");
   }
   memset((void *)mem_p, 0, (dim1 * byte_depth));

   *ptr = mem_p;
}

/*****************************************************************************/
/* STATIC                        allocate_2d                                 */
/*****************************************************************************/

static void allocate_2d(void ***ptr, int dim1, int dim2, int byte_depth)
{
   int i;
   void **mem_p;

   if ((mem_p = (void **)local_malloc("KLT", (dim1 * sizeof(void *)))) == NULL) {
      local_error("allocate_2d():  calloc failure (1).");
   }
   if ((mem_p[0] = (void *)local_malloc("KLT", (dim1 * dim2 * byte_depth))) == NULL) {
      local_error("allocate_2d():  calloc failure (2).");
   }
   for (i=1; i<dim1; i++)
      mem_p[i] = (char *)mem_p[0] + (i * dim2 * byte_depth);
   memset((void *)mem_p[0], 0, (dim1 * dim2 * byte_depth));

   *ptr = mem_p;
}

/*****************************************************************************/
/* STATIC                        deallocate_1d                               */
/*****************************************************************************/

static void deallocate_1d(void *ptr)
{
   if (ptr != NULL) {
      local_free(ptr);
   }
}

/*****************************************************************************/
/* STATIC                        deallocate_2d                               */
/*****************************************************************************/

static void deallocate_2d(void **ptr)
{
   if (ptr != NULL) {
      local_free(ptr[0]);
      local_free(ptr);
   }
}

/*****************************************************************************/
/* STATIC                calculate_cov_cor_mean_bil                          */
/*****************************************************************************/

static void calculate_cov_cor_mean_bil(image_reader_ref reader, 
				       int rows, int cols, int bands,
				       float **cov, float **cor, float *mean)
{
   int l,p,b1,b2,elem;
   void **in_page;
   float **line_slice;
   double **Sxy, *Sx;
   clock_t start;

   allocate_2d((void ***)&in_page, bands, cols, sizeof(ifc_int));
   allocate_2d((void ***)&line_slice, bands, cols, sizeof(float));
   allocate_2d((void ***)&Sxy, bands, bands, sizeof(double));
   allocate_1d((void **)&Sx, bands, sizeof(double));
   
   start = clock();

   local_printf(0,70, "\nCalculating covariance matrix ...");
   for (l=0; l<rows; l++) {

#undef printf
      printf("Line %d/%d\r", l+1, rows); fflush(stdout);
#define printf      __you_must_use_local_printf_instead__

      for (b1=0; b1<bands; b1++) {
	 reader->pull_line(reader, (ifc_int *)in_page[b1], b1, cols);
      }

      COPY_VOID_TO_FLOAT(line_slice, in_page, ifc_int, bands, cols);

      /* First band (b1 = 0) calculations include Sx and Sxy */
      {
	 float *y = *line_slice;
	 double *Sxy_p = *Sxy;
	 double *Sx_p = Sx;
	 for (b2=0; b2<bands; b2++) {
	    float *x = *line_slice;
	    float cur_Sx = 0.0;
	    float cur_Sxy = 0.0;
	    for (p=0; p<cols; p++) {
	       float val = *y++;
	       cur_Sx += val;
	       cur_Sxy += (*x++ * val);
	    }
	    *Sx_p++ += cur_Sx;
	    *Sxy_p++ += cur_Sxy;
	 }
      }

      /* The rest of the band calculations (1 <= b1 < bands) are for Sxy only */
      for (b1=1; b1<bands; b1++) {
	 float *y = line_slice[b1];
	 double *Sxy_p = &Sxy[b1][b1];
	 for (b2=b1; b2<bands; b2++) {
	    float *x = line_slice[b1];
	    float cur_Sxy = 0.0;
	    for (p=0; p<cols; p++) {
	       cur_Sxy += (*x++ * *y++);
	    }
	    *Sxy_p++ += cur_Sxy;
	 }
      }
   }
   local_printf(0,70, "\n");

   local_printf(0,70, "\nCovariance matrix calculated in ... %.5f seconds.\n", 
		(float)(clock() - start)/(float)CLOCKS_PER_SEC);

   elem = (rows * cols);