calc_klt.c

JPEG2000实现的源码

   for (b1=0; b1<bands; b1++) {
      mean[b1] = (float)(Sx[b1] / (double)elem);
      for (b2=b1; b2<bands; b2++) {
	 cov[b1][b2] = (float)((Sxy[b1][b2] - Sx[b1]*Sx[b2]/(double)elem) / (double)(elem - 1));
	 cov[b2][b1] = cov[b1][b2];
      }
   }

   for (b1=0; b1<bands; b1++) {
      for (b2=b1; b2<bands; b2++) {
	 double num = cov[b1][b2];
	 double den = sqrt(cov[b1][b1] * cov[b2][b2]);
	 cor[b1][b2] = (float)((den == 0.0) ? 0.0 : (num / den));
	 cor[b2][b1] = cor[b1][b2];
      }
   }

   deallocate_2d((void **)in_page);
   deallocate_2d((void **)line_slice);
   deallocate_2d((void **)Sxy);
   deallocate_1d((void *)Sx);
}

/*****************************************************************************/
/* STATIC                           sort_eigen                               */
/*****************************************************************************/

static void sort_eigen(double **vec, double *val, int n)
{
   int i,j;
   int max_ind;
   double max_val;
   
   for (i=0; i<(n-1); i++) {
      for (j=(i+1), max_ind=-1, max_val=val[i]; j<n; j++) {
	 if (val[j] >= max_val) {
	    max_val = val[j];
	    max_ind = j;
	 }
      }
      if (max_ind >= 0) {
	 val[max_ind] = val[i];
	 val[i] = max_val;
	 for (j=0; j<n; j++) {
	    max_val = vec[j][max_ind];
	    vec[j][max_ind] = vec[j][i];
	    vec[j][i] = max_val;
	 }
      }
   }
}

/*****************************************************************************/
/* STATIC                        calculate_eigen                             */
/*****************************************************************************/

static void calculate_eigen(float **cov, float **eigen_vecs, float *eigen_vals, 
			    int n, int transpose)
{
   int i,j;
   MAT *A, *Q;
   VEC *out;

   A = m_get(n, n);
   Q = m_get(n, n);
   out = v_get(n);

   for (i=0; i<n; i++) {
      for (j=0; j<n; j++) {
	 A->me[i][j] = cov[i][j];
      }
   }

   symmeig(A, Q, out);
   sort_eigen(Q->me, out->ve, n);

   if (transpose) {
      for (i=0; i<n; i++) {
	 eigen_vals[i] = (float)out->ve[i];
	 for (j=0; j<n; j++) {
	    eigen_vecs[j][i] = (float)Q->me[i][j];
	 }
      }
   } else {
      for (i=0; i<n; i++) {
	 eigen_vals[i] = (float)out->ve[i];
	 for (j=0; j<n; j++) {
	    eigen_vecs[i][j] = (float)Q->me[i][j];
	 }
      }
   }

   m_free(A);
   m_free(Q);
   v_free(out);
}

/*****************************************************************************/
/* STATIC                           usage                                    */
/*****************************************************************************/

static void usage(char *fname)
{
   image_reader_ref image;
   extern image_reader_ref create_simple_image_reader(void);
   extern image_reader_ref create_mcraw_image_reader(void);

   local_printf(0,70, "\nUsage: %s\n", fname);
   local_printf(0,70, "\n");
   
   image = create_simple_image_reader();
   print_generic_usage(image->get_usage(image));
   image->terminate(image);

   local_printf(0,70, "...OR...");
   
   image = create_mcraw_image_reader();
   print_generic_usage(image->get_usage(image));
   image->terminate(image);

   local_printf(4,70, "-o <name of file for forward KLT result>");
   local_printf(4,70, "-eigval <name of text file for eigenvalues> (default = stdout)");

   local_exit(0);
}   

/* ========================================================================= */
/* -------------------------- External Functions --------------------------- */
/* ========================================================================= */

/*****************************************************************************/
/* EXTERN                             main                                   */
/*****************************************************************************/

int main(int argc, char *argv[])
{
   cmdl_ref cmdl;
   extern cmdl_ref create_std_cmdl(void);
   
   int i,p;
   char **params;
   char *out_fname = NULL, *eigval_fname = NULL;
   image_reader_ref reader;
   int rows, cols, bands;
   float **cov, **cor, *mean, **eigen_vecs, *eigen_vals;
   canvas_dims dims = {0};
   FILE *fp;

   if (argc == 1) {
      usage(argv[0]);
   }

   /* Create `cmdl' object */

   cmdl = create_std_cmdl();
   cmdl->initialize(cmdl, argc, argv);
   
   if ((p = cmdl->extract(cmdl,"-o",-1,&params)) == 1) {
      out_fname = params[0];
   } else {
      local_error("The `-o' argument is invalid.");
   }
   if ((p = cmdl->extract(cmdl,"-eigval",-1,&params)) >= 0) {
      if (p == 1) {
	 eigval_fname = params[0];
      } else {
	 local_error("The `-eigval' argument is invalid.");
      }
   }

   reader = construct_reader(cmdl);
   reader->initialize(reader, &bands, cmdl);

   cmdl->terminate(cmdl, 1 /* check args */);
   cmdl = NULL;

   if (out_fname == NULL) {
      local_error("Insufficient inputs on command-line.");
   } else {
      if ((fp = fopen(out_fname, "wb")) == NULL) {
	 local_error("Unable to open output file, %s.", out_fname);
      }
   }
   if (bands <= 1) {
      local_error("Input does not represent a multi-component image.");
   }

   reader->get_dims(reader, 0, &dims);
   for (i=0; i<bands; i++) {
      canvas_dims cur_dims = {0};
      reader->get_dims(reader, i, &cur_dims);
      if ((cur_dims.rows != dims.rows) || (cur_dims.cols != dims.cols)) {
	 local_error("Component dimensions must be the same!");
      }
   }
   rows = dims.rows;
   cols = dims.cols;

   allocate_2d((void ***)&cov, bands, bands, sizeof(float));
   allocate_2d((void ***)&cor, bands, bands, sizeof(float));
   allocate_1d((void **)&mean, bands, sizeof(float));
   allocate_2d((void ***)&eigen_vecs, bands, bands, sizeof(float));
   allocate_1d((void **)&eigen_vals, bands, sizeof(float));

   calculate_cov_cor_mean_bil(reader, rows, cols, bands, cov, cor, mean);
   calculate_eigen(cov, eigen_vecs, eigen_vals, bands, 1 /* transpose matrix */);

   if (fwrite((void *)eigen_vecs[0], sizeof(float), (bands * bands), fp) !=
       (unsigned int)(bands * bands)) {
      local_error("Error writing to output file, %s.", out_fname);
   }
   fclose(fp);

   if (eigval_fname != NULL) {
      if ((fp = fopen(eigval_fname, "w")) == NULL) {
	 local_error("Unable to open eigenvalue file, %s.", eigval_fname);
      }
      for (i=0; i<bands; i++) {
	 fprintf(fp, "%20.4f\n", eigen_vals[i]);
      }
      fclose(fp);
   } else {
      local_printf(0,70, "\n\nEigenvalues:\n");
      for (i=0; i<bands; i++) {
	 local_printf(4,70, "%20.4f\n", eigen_vals[i]);
      }
      local_printf(0,70, "\n");
   }

   reader->terminate(reader);
   deallocate_2d((void **)cov);
   deallocate_2d((void **)cor);
   deallocate_1d((void *)mean);
   deallocate_2d((void **)eigen_vecs);
   deallocate_1d((void *)eigen_vals);

   return 0;
}