X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fimagefile.cpp;h=84c66741ec86b223bdac6254ebd6b6bc6be0bc15;hp=d4c973cc3a4a1aeeb1f7be4054c9f9ade7f581ff;hb=67a6c34b5a6f38d34e8cbe66091853f453fd2d7a;hpb=fd1d136a94a6d20013f38d6a997bdfefad0f5e98 diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index d4c973c..84c6674 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.22 2000/12/20 20:08:48 kevin Exp $ +** $Id: imagefile.cpp,v 1.26 2000/12/29 19:30:08 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,37 +214,58 @@ 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 (nx * ny); + std::vector vecImage; + int iVec = 0; + vecImage.resize (nx * ny); for (int ix = 0; ix < nx; ix++) { for (int iy = 0; iy < ny; iy++) - vecImage.push_back (v[ix][iy]); + vecImage[iVec++] = v[ix][iy]; } - vectorNumericStatistics (vecImage, min, max, mean, mode, median, stddev); -} - - + vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev); +} + void ImageFile::getMinMax (double& min, double& max) const { @@ -264,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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumn vCol = vImag[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumn vRealCol = vReal[ix]; + ImageFileColumn vImagCol = vImag[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in1 = vLHS[ix]; + ImageFileColumnConst in2 = vRHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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(); + + 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 (int ix = 0; ix < m_nx; ix++) + for (int iy = 0; iy < m_ny; iy++) + vImagResult[ix][iy] = 0; + } + + for (int ix = 0; ix < m_nx; ix++) + for (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 (int ix = 0; ix < m_nx; ix++) + for (int iy = 0; iy < m_ny; iy++) + vImagResult[ix][iy] = 0; + } + + for (int ix = 0; ix < m_nx; ix++) + for (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 (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (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) {