X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;ds=sidebyside;f=libctsim%2Fimagefile.cpp;h=b6525e1fda69077f252b2994aede0cb997d01d21;hb=7ae47cb0ff0a16d1c36797576155263434cc73ff;hp=b86ff794ae9cd0242201c21e33184bca9c513e3f;hpb=f7d2b7144f32a7bd157b7689022e62944b82fcc1;p=ctsim.git diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index b86ff79..b6525e1 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: imagefile.cpp,v 1.23 2000/12/21 03:40:58 kevin Exp $ +** $Id: imagefile.cpp,v 1.27 2000/12/29 20:15:37 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -28,8 +28,8 @@ #include "ct.h" -F32Image::F32Image (int nx, int ny) -: Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32) +F32Image::F32Image (int nx, int ny, int dataType) +: Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType) { } @@ -38,10 +38,11 @@ F32Image::F32Image (void) { setPixelFormat (Array2dFile::PIXEL_FLOAT32); setPixelSize (sizeof(kfloat32)); + setDataType (Array2dFile::DATA_TYPE_REAL); } -F64Image::F64Image (int nx, int ny) -: Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64) +F64Image::F64Image (int nx, int ny, int dataType) +: Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType) { } @@ -50,6 +51,7 @@ F64Image::F64Image (void) { setPixelFormat (PIXEL_FLOAT64); setPixelSize (sizeof(kfloat64)); + setDataType (Array2dFile::DATA_TYPE_REAL); } void @@ -62,7 +64,7 @@ ImageFile::filterResponse (const char* const domainName, double bw, const char* for (int i = -hx; i <= hx; i++) { for (int j = -hy; j <= hy; j++) { - double r = sqrt (i * i + j * j); + double r = ::sqrt (i * i + j * j); v[i+hx][j+hy] = filter.response (r); } @@ -170,7 +172,7 @@ ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, } } - d = sqrt (sqErrorSum / sqDiffFromMeanSum); + d = ::sqrt (sqErrorSum / sqDiffFromMeanSum); r = absErrorSum / absValueSum; int hx = m_nx / 2; @@ -212,26 +214,46 @@ ImageFile::printStatistics (std::ostream& os) const { double min, max, mean, mode, median, stddev; - statistics (min, max, mean, mode, median, stddev); - + statistics (min, max, mean, mode, median, stddev); + if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) + os << "Real Component Statistics" << std::endl; + os << " min: " << min << std::endl; os << " max: " << max << std::endl; os << " mean: " << mean << std::endl; os << " mode: " << mode << std::endl; os << "median: " << median << std::endl; - os << "stddev: " << stddev << std::endl; + os << "stddev: " << stddev << std::endl; + + if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) { + statistics (getImaginaryArray(), min, max, mean, mode, median, stddev); + os << std::endl << "Imaginary Component Statistics" << std::endl; + os << " min: " << min << std::endl; + os << " max: " << max << std::endl; + os << " mean: " << mean << std::endl; + os << " mode: " << mode << std::endl; + os << "median: " << median << std::endl; + os << "stddev: " << stddev << std::endl; + } } void ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const -{ - int nx = m_nx; - int ny = m_ny; - ImageFileArrayConst v = getArray(); - - if (v == NULL || nx == 0 || ny == 0) - return; +{ + ImageFileArrayConst v = getArray(); + statistics (v, min, max, mean, mode, median, stddev); +} + + +void +ImageFile::statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const +{ + int nx = m_nx; + int ny = m_ny; + + if (v == NULL || nx == 0 || ny == 0) + return; std::vector vecImage; int iVec = 0; @@ -242,9 +264,8 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou } vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev); -} - - +} + void ImageFile::getMinMax (double& min, double& max) const { @@ -266,7 +287,402 @@ ImageFile::getMinMax (double& min, double& max) const } } } + +bool +ImageFile::convertRealToComplex () +{ + if (dataType() != Array2dFile::DATA_TYPE_REAL) + return false; + + if (! reallocRealToComplex()) + return false; + + ImageFileArray vImag = getImaginaryArray(); + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumn vCol = vImag[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + *vCol++ = 0; + } + + return true; +} + +bool +ImageFile::convertComplexToReal () +{ + if (dataType() != Array2dFile::DATA_TYPE_COMPLEX) + return false; + + ImageFileArray vReal = getArray(); + ImageFileArray vImag = getImaginaryArray(); + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumn vRealCol = vReal[ix]; + ImageFileColumn vImagCol = vImag[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) { + std::complex c (*vRealCol, *vImagCol); + *vRealCol++ = std::abs (c); + vImagCol++; + } + } + + return reallocComplexToReal(); +} + +bool +ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const +{ + if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArrayConst vRHS = rRHS.getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + *out++ = *in1++ - *in2++; + } + + return true; +} +bool +ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const +{ + if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArrayConst vRHS = rRHS.getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + *out++ = *in1++ + *in2++; + } + + return true; +} + +bool +ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const +{ + if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArrayConst vRHS = rRHS.getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + *out++ = *in1++ * *in2++; + } + + return true; +} + +bool +ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const +{ + if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArrayConst vRHS = rRHS.getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) { + if (*in2 != 0.) + *out++ = *in1++ / *in2++; + else + *out++ = 0; + } + } + + return true; +} + + +bool +ImageFile::invertPixelValues (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + *out++ = - *in++; + } + + return true; +} + +bool +ImageFile::sqrt (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + if (*in < 0) + *out++ = -::sqrt(-*in++); + else + *out++ = ::sqrt(*in++); + } + + return true; +} + +bool +ImageFile::log (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + if (*in <= 0) + *out++ = 0; + else + *out++ = ::log(*in++); + } + + return true; +} + +bool +ImageFile::exp (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) + *out++ = ::exp (*in++); + } + + return true; +} + +bool +ImageFile::inverseFourier (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + return true; +} + +bool +ImageFile::fourier (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + if (result.dataType() == Array2dFile::DATA_TYPE_REAL) { + if (! result.reallocRealToComplex ()) + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vRealResult = result.getArray(); + ImageFileArray vImagResult = result.getImaginaryArray(); + + unsigned int ix, iy; + double* pY = new double [m_ny]; + std::complex** complexOut = new std::complex* [m_nx]; + for (ix = 0; ix < m_nx; ix++) + complexOut[ix] = new std::complex [m_ny]; + + for (ix = 0; ix < m_nx; ix++) { + for (iy = 0; iy < m_ny; iy++) + pY[iy] = vLHS[ix][iy]; + ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD); + } + delete pY; + + std::complex* pX = new std::complex [m_nx]; + std::complex* complexOutCol = new std::complex [m_nx]; + for (iy = 0; iy < m_ny; iy++) { + for (ix = 0; ix < m_nx; ix++) + pX[ix] = complexOut[ix][iy]; + ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD); + for (ix = 0; ix < m_nx; ix++) + complexOut[ix][iy] = complexOutCol[ix]; + } + delete [] pX; + + for (ix = 0; ix < m_nx; ix++) + ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny); + + + for (iy = 0; iy < m_ny; iy++) { + for (ix = 0; ix < m_nx; ix++) + complexOutCol[ix] = complexOut[ix][iy]; + ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);; + for (ix = 0; ix < m_nx; ix++) + complexOut[ix][iy] = complexOutCol[ix]; + + } + delete [] complexOutCol; + + for (ix = 0; ix < m_nx; ix++) + for (iy = 0; iy < m_ny; iy++) { + vRealResult[ix][iy] = complexOut[ix][iy].real(); + vImagResult[ix][iy] = complexOut[ix][iy].imag(); + } + + for (ix = 0; ix < m_nx; ix++) + delete [] complexOut[ix]; + + delete [] complexOut; + + return true; +} + +bool +ImageFile::magnitude (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArray vReal = getArray(); + ImageFileArray vImag = NULL; + if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) + vImag = getImaginaryArray(); + + ImageFileArray vRealResult = result.getArray(); + if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) { + ImageFileArray vImagResult = result.getImaginaryArray(); + for (unsigned int ix = 0; ix < m_nx; ix++) + for (unsigned int iy = 0; iy < m_ny; iy++) + vImagResult[ix][iy] = 0; + } + + for (unsigned int ix = 0; ix < m_nx; ix++) + for (unsigned int iy = 0; iy < m_ny; iy++) { + if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) + vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]); + else + vRealResult[ix][iy] = vReal[ix][iy]; + } + + return true; +} + +bool +ImageFile::phase (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArray vReal = getArray(); + ImageFileArray vImag = NULL; + if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) + vImag = getImaginaryArray(); + + ImageFileArray vRealResult = result.getArray(); + if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) { + ImageFileArray vImagResult = result.getImaginaryArray(); + for (unsigned int ix = 0; ix < m_nx; ix++) + for (unsigned int iy = 0; iy < m_ny; iy++) + vImagResult[ix][iy] = 0; + } + + for (unsigned int ix = 0; ix < m_nx; ix++) + for (unsigned int iy = 0; iy < m_ny; iy++) { + if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) + vRealResult[ix][iy] = ::atan (vImag[ix][iy] / vReal[ix][iy]); + else + vRealResult[ix][iy] = vReal[ix][iy]; + } + + return true; +} + +bool +ImageFile::square (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (unsigned int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (unsigned int iy = 0; iy < m_ny; iy++) { + *out++ = *in * *in; + in++; + } + } + + return true; +} + + void ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax) {