image.cc

来自「用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵」· CC 代码 · 共 1,083 行 · 第 1/3 页

// Copyright (C) 2002 Charless C. Fowlkes <[e-mail]>
// Copyright (C) 2002 David R. Martin <[e-mail]>
//
// 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., 59 Temple Place - Suite 330, Boston, MA
// 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html.

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <values.h>
#include <assert.h>
extern "C"
{
#include <jpeglib.h>
#include <jerror.h>
}

#include "util.hh"
#include "image.hh"
#include "sort.hh"


namespace Util
{
  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // read in a jpeg file into an ImageStack. 
  // if the jpeg file is grayscale then the resulting ImageStack only has 1 layer
  // otherwise it has 3 layers corresponding to the RGB colorspace.
  //

  bool readJpegFile(const char *filename, ImageStack& im)
  {
      assert(filename != NULL);

      int bytesPerPixel;
      int height;
      int width;
      unsigned char* imbuf;
      bool isGray;

      FILE* file = fopen(filename, "r");
      if (!file)
      {
        return (false);
      }

      jpeg_decompress_struct cinfo;
      jpeg_error_mgr jerr;
      cinfo.err = jpeg_std_error (&jerr);
      jpeg_create_decompress (&cinfo);
      jpeg_stdio_src (&cinfo, file);
      jpeg_read_header (&cinfo, (boolean) true);

      if (cinfo.out_color_space == JCS_GRAYSCALE)
      {
          isGray = true;
          cinfo.output_components = 1;
          bytesPerPixel = 1;
      }
      else
      {
          isGray = false;
          cinfo.out_color_space = JCS_RGB;
          cinfo.output_components = 3;
          bytesPerPixel = 3;
      }
      jpeg_calc_output_dimensions (&cinfo);
      jpeg_start_decompress (&cinfo);

      height = cinfo.output_height;
      width = cinfo.output_width;
      imbuf = new unsigned char[width * height * bytesPerPixel];

      const int lineSize = width * bytesPerPixel;
      unsigned char* p = imbuf;

      while (cinfo.output_scanline < cinfo.output_height)
      {
        jpeg_read_scanlines (&cinfo, &p, 1);
        p += lineSize;
      }

      jpeg_finish_decompress (&cinfo);
      jpeg_destroy_decompress (&cinfo);
      fclose (file);

      //now fill in our image array data structure
      if (isGray)
      {
        im.resize(1,width,height);
        for (int i = 0; i < width * height; i++)
        {
          const float gray = (float) imbuf[i] / 255;
          const int y = i / width;
          const int x = i % width;
          im(0,x,y) = gray;
          // this fails on fp type errors:  assert (gray >= 0 && gray <= 1);
        }
      }
      else
      {
        im.resize(3,width,height);
        for (int i = 0; i < width * height; i++)
        {
          const int y = i / width;
          const int x = i % width;
          float r = (float) imbuf[3 * i] / 255;
          float g = (float) imbuf[3 * i + 1] / 255;
          float b = (float) imbuf[3 * i + 2] / 255;
          if ((r < 0) || (r > 1))
          {
            std::cerr << "r = " << r << std::endl;
          }
          if ((g < 0) || (g > 1))
          {
            std::cerr << "g = " << g << std::endl;
          }
          if ((b < 0) || (b > 1))
          {
            std::cerr << "b = " << r << std::endl;
          }
          assert (r >= 0 && r <= 1);
          assert (g >= 0 && g <= 1);
          assert (b >= 0 && b <= 1);
          im(RGB_R,x,y) = r;
          im(RGB_G,x,y) = g;
          im(RGB_B,x,y) = b;
        }
      }
      delete[] imbuf;

      return true;
  }

  //////////////////////////////////////////////////////////////////////////////////////////////////////

  //
  // used for writing out a single channel image with
  // the "jet" psuedo-color map
  //
  static int jetR (float v)
  {
    assert (v >= 0 && v <= 1);
    int i = (uint) rint (v * 255);
    return Util::minmax (0, (450 - 5 * abs (i - 196)), 255);
  }

  static int jetG (float v)
  {
    assert (v >= 0 && v <= 1);
    int i = (uint) rint (v * 255);
    return Util::minmax (0, (450 - 5 * abs (i - 128)), 255);
  }

  static int jetB (float v)
  {
      assert (v >= 0 && v <= 1);
      int i = (uint) rint (v * 255);
      return Util::minmax (0, (450 - 5 * abs (i - 64)), 255);
  }

  //
  // create a normalized version of this image.
  // first run thru and find the min and max values
  // and then create a new image whose values range 
  // over [0.0,1.0].  the constant image is set to 0.0
  //
  // used in the Jpeg reading and writing.
  //
  static void getNormalized(const Image& im, Image& normalized) 
  {
      int width = im.size(1);
      int height = im.size(0);
      normalized.resize(width,height);

      float max = im(0,0);
      float min = im(0,0);

      for (int i = 0; i < width; i++)
      {
        for (int j = 0; j < height; j++)
        {
          float val = im(i,j);
          max = Util::max(max, val);
          min = Util::min(min, val);
        }
      }

      if ((max - min) > 0)
      {
        normalized = (im - min) / (max - min);
      }
      else
      {
        normalized.init(0);
      }
  }

  //
  // write out a grayscale image to a jpeg file.
  // if normalize=true, then the range of the image is adjusted to use the full scale
  // if jet=true then the image is written in pseudocolor rather than grayscale
  //
  bool writeJpegFile(const Image& im, const char *filespec, const bool normalize, const bool jet)
  {
      assert(filespec != NULL);

      //normalized version of this image
      Image normim;
      if (normalize)
      {
        getNormalized(im,normim); 
      }
      else
      {
        normim = im;
      }

      struct jpeg_error_mgr jerr;
      struct jpeg_compress_struct cinfo;

      // Open the output file.
      FILE* file = Util::openOutputStrm (filespec);
      if (!file)
      {
        return (false);
      }

      // Set up the normal JPEG error handling.
      cinfo.err = jpeg_std_error (&jerr);

      // Init a JPEG compression object.
      jpeg_create_compress (&cinfo);

      // Specify source of data.
      jpeg_stdio_dest (&cinfo, file);

      // Specify compression parameters.
      cinfo.image_width = im.size(0);
      cinfo.image_height = im.size(1);
      cinfo.input_components = 3;
      cinfo.in_color_space = JCS_RGB;
      jpeg_set_defaults (&cinfo);

      // Start compression.
      jpeg_start_compress (&cinfo, TRUE);

      // Allocate space for one scanline.
      JSAMPLE* buf = new JSAMPLE[im.size(0) * 3];

      // Write data one scanline at a time.
      for (int y = 0; y < im.size(1); y++)
      {
        for (int x = 0; x < im.size(0); x++)
        {
          float v = normim(x,y);
          assert (v >= 0 && v <= 1);
          if (jet)
          {
              buf[3*x + 0] = jetR(v);
              buf[3*x + 1] = jetG(v);
              buf[3*x + 2] = jetB(v);
          }
          else
          {
              int g = (int)rint(255 * v);
              buf[3*x + 0] = g;
              buf[3*x + 1] = g;
              buf[3*x + 2] = g;
          }
        }
        int c = jpeg_write_scanlines (&cinfo, &buf, 1);
        if (c != 1)
        {
          fclose (file);
          jpeg_destroy_compress (&cinfo);
          delete[]buf;
          return false;
        }
      }

      // Clean up.
      jpeg_finish_compress (&cinfo);
      fclose (file);
      jpeg_destroy_compress (&cinfo);
      delete[]buf;
      return true;
  }

  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // convert an RGB imagestack into an 1976 CIE L*a*b* imagestack.
  //

  void rgb2lab (const ImageStack& rgb, ImageStack& lab)
  {
      assert(rgb.size(0) == 3);
      const int w = rgb.size(1);
      const int h = rgb.size(2);

      lab.resize(3,w,h);

      for (int i = 0; i < w; i++)
      {
        for (int j = 0; j < h; j++)
        {
            // RGB
            const float R = rgb(RGB_R,i,j);
            const float G = rgb(RGB_G,i,j);
            const float B = rgb(RGB_B,i,j);
            assert (R >= 0 && R <= 1);
            assert (G >= 0 && G <= 1);
            assert (B >= 0 && B <= 1);
            // RGB --> XYZ
            const float X = 0.412453 * R + 0.357580 * G + 0.180423 * B;
            const float Y = 0.212671 * R + 0.715160 * G + 0.072169 * B;
            const float Z = 0.019334 * R + 0.119193 * G + 0.950227 * B;
            // XYZ of D65 reference white (R=G=B=1).
            const float Xn = 0.950456;
            const float Yn = 1.000000;
            const float Zn = 1.088754;
            // XYZ --> 1976 CIE L*a*b*
            const float rX = X / Xn;
            const float rY = Y / Yn;
            const float rZ = Z / Zn;
            const float thresh = 0.008856;
#define f(t) (t > thresh) ? pow(t,1./3.) : (7.787 * t + 16./116.)
            const float fX = f(rX);
            const float fY = f(rY);
            const float fZ = f(rZ);
#undef f
            const float L = (rY > thresh) ? 116. * pow (rY, 1. / 3.) - 16. : 903.3 * rY;
            const float a = 500. * (fX - fY);
            const float b = 200. * (fY - fZ);
            assert (L >= 0 && L <= 100);
            assert (a >= -120 && a <= 120);
            assert (b >= -120 && b <= 120);
            lab(LAB_L,i,j) = L;
            lab(LAB_A,i,j) = a;
            lab(LAB_B,i,j) = b;
        }
      }
  }

  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // normalize an lab image so that values lie in [0,1]
  //