diff --git a/.gitignore b/.gitignore index d1d270a..bf6c02a 100644 --- a/.gitignore +++ b/.gitignore @@ -5,13 +5,11 @@ testData test_data/psf.dv test_data/GPUsirecon/* test_data/proc.dv -IVE lapack *dist fftw2 cmake_build -src/IVE _test.sh src/cudaSirecon/cudasireconConfig.h src/cudasireconConfig.h diff --git a/README.md b/README.md index e00e338..5451472 100644 --- a/README.md +++ b/README.md @@ -258,40 +258,3 @@ conda env create -f environment-windows.yml conda activate simbuild bld.bat # see bld.bat for cmake details ``` - -### Building with MRC support - -The IVE/Priism libraries (for MRC/DV support), are not distributed with this -source code and must be acquired seperately from UCSF. Place them in a folder -called `IVE` at the top level of the source folder (same folder as the -cudaSirecon folder). It should minimally have the following files and folders -(example shown for linux, use `.a` or `.lib` as necessary for osx or windows) - -then build with `cmake -DBUILD_MRC=ON ....` - -
- -required IVE folder structure - -``` -cudasirecon -├── ... -└── IVE/ - ├── darwin64/ - │ ├── INCLUDE/ - │ └── LIB/ - ├── linux64/ - │ ├── INCLUDE/ - │ │ ├── IM.h - │ │ ├── IMInclude.h - │ │ └── IWApiConstants.h - │ └── LIB/ - │ ├── libimlib.a - │ └── libive.a - └── win64/ - ├── INCLUDE/ - └── LIB/ -``` - -
- diff --git a/VStutio_project/VStutio_project.vcxproj b/VStutio_project/VStutio_project.vcxproj index 5c9b677..7a29bae 100644 --- a/VStutio_project/VStutio_project.vcxproj +++ b/VStutio_project/VStutio_project.vcxproj @@ -144,7 +144,7 @@ true true WIN32;NDEBUG;_CONSOLE;_LIB;%(PreprocessorDefinitions) - $(CUDA_PATH)/include;../Buffers;C:/libtiff;../IVE/win64/INCLUDE;C:/Boost/boost_1_55_0 + $(CUDA_PATH)/include;../Buffers;C:/libtiff;C:/Boost/boost_1_55_0 None @@ -152,8 +152,8 @@ true true true - C:\Boost\boost_1_55_0\lib;$(CUDA_PATH)\lib\x64;C:\libtiff;..\lapack\win64;..\IVE\win64\lib;$(SolutionDir) - cuda.lib;cudart.lib;cufft.lib;libtiff.lib;liblapack.lib;libblas.lib;libimlib.lib;libive.lib;Buffers.lib;%(AdditionalDependencies) + C:\Boost\boost_1_55_0\lib;$(CUDA_PATH)\lib\x64;C:\libtiff;..\lapack\win64;$(SolutionDir) + cuda.lib;cudart.lib;cufft.lib;libtiff.lib;liblapack.lib;libblas.lib;Buffers.lib;%(AdditionalDependencies) 64 diff --git a/bld.bat b/bld.bat index eab9062..702e216 100644 --- a/bld.bat +++ b/bld.bat @@ -1,19 +1,15 @@ -del /Q cmake_build +@REM rmdir /S /Q cmake_build -@REM if not exist src\IVE\ ( -@REM powershell -Command "Invoke-WebRequest https://www.dropbox.com/s/2twvw0go3dr3aim/IVE.zip -OutFile IVE.zip" -@REM powershell -Command "Expand-Archive IVE.zip -DestinationPath src\" -@REM del IVE.zip ) +call "C:\Program Files (x86)\Microsoft Visual Studio\2019\Community\VC\Auxiliary\Build\vcvars64.bat" -mkdir cmake_build +@REM mkdir cmake_build cd cmake_build -cmake -G Ninja ^ - -DBUILD_MRC=ON ^ - -DBUILD_OTF_VIEWER=OFF ^ - -DCMAKE_BUILD_TYPE=Release ^ - -DCMAKE_INSTALL_PREFIX="%CONDA_PREFIX%/Library" ^ - ../src +@REM cmake -G Ninja ^ +@REM -DBUILD_OTF_VIEWER=OFF ^ +@REM -DCMAKE_BUILD_TYPE=Release ^ +@REM -DCMAKE_INSTALL_PREFIX="%CONDA_PREFIX%/Library" ^ +@REM ../src ninja -ninja install +@REM ninja install diff --git a/build.sh b/build.sh index 9be1703..7ab0d64 100755 --- a/build.sh +++ b/build.sh @@ -1,23 +1,14 @@ rm -rf cmake_build -if [ ! -d "src/IVE" ] -then - wget src/ https://www.dropbox.com/s/2twvw0go3dr3aim/IVE.zip - unzip -o -d src IVE.zip - rm IVE.zip -fi - - mkdir cmake_build cd cmake_build CXXFLAGS="$CXXFLAGS -Wfatal-errors -Wno-deprecated-declarations" cmake ${CMAKE_ARGS} \ - -DBUILD_MRC=ON \ -DBUILD_OTF_VIEWER=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX \ ../src -make -j 2 -make install +make -j 4 +# make install diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index db5a883..d90598c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -9,14 +9,6 @@ # includes the nvcc compiler For windows, you will need to install the CUDA # toolkit manually # -# 1. MRC/DV support (optional) If you want to build with MRC support add the IVE -# libraries to the src directory with the following paths: For Linux: -# src/IVE/linux64/INCLUDE src/IVE/linux64/LIB For Windows: -# src/IVE/win64/INCLUDE src/IVE/win64/LIB -# -# If you DON'T want to build with MRC or do not have access to the IVE libraries -# please run cmake with the `-DBUILD_MRC=OFF` option -# # To build the otfviewer, use `-DBUILD_OTF_VIEWER=ON` # # 1. run build.sh (linux) or bld.bat (windows) inspect those scripts for more @@ -81,33 +73,6 @@ endif(WIN32) message(STATUS "CUDA_NVCC_FLAGS: " ${CUDA_NVCC_FLAGS}) -if(BUILD_MRC OR BUILD_OTF_VIEWER) - message(STATUS "** Building WITH mrc support **") - - if(WIN32) - set(PLATFORM win64) - else() - set(PLATFORM linux64) - endif(WIN32) - - set(PRIISM_LIB_PATH "${CMAKE_SOURCE_DIR}/IVE/${PLATFORM}/LIB") - set(PRIISM_INCLUDE_PATH "${CMAKE_SOURCE_DIR}/IVE/${PLATFORM}/INCLUDE") - include_directories(${PRIISM_INCLUDE_PATH}) - link_directories(${PRIISM_LIB_PATH}) - - find_library( - IMLIB - NAMES imlib libimlib - PATHS ${PRIISM_LIB_PATH} REQUIRED) - find_library( - IVELIB - NAMES ive libive - PATHS ${PRIISM_LIB_PATH} REQUIRED) - add_definitions(-DMRC) -else() - message(STATUS "** Building WITHOUT mrc support **") -endif() - add_subdirectory(Buffers) add_subdirectory(cudaSirecon) add_subdirectory(otf) diff --git a/src/cudaSirecon/CMakeLists.txt b/src/cudaSirecon/CMakeLists.txt index 0ce959a..09164c4 100644 --- a/src/cudaSirecon/CMakeLists.txt +++ b/src/cudaSirecon/CMakeLists.txt @@ -3,8 +3,6 @@ include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/../Buffers") find_package(LAPACK REQUIRED) -message(STATUS "IMLIB: " ${IMLIB}) -message(STATUS "IVELIB: " ${IVELIB}) message(STATUS "LEG_STDIO: " ${LEG_STDIO}) message(STATUS "LAPACK_FOUND: " ${LAPACK_FOUND}) @@ -23,8 +21,7 @@ target_link_libraries( ${BOOST_LIBRARIES} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_FILESYSTEM_LIBRARY} - ${IMLIB} - ${IVELIB}) +) # cudasirecon executable file diff --git a/src/cudaSirecon/SIM_reconstructor.hpp b/src/cudaSirecon/SIM_reconstructor.hpp index 973c057..9f3de70 100644 --- a/src/cudaSirecon/SIM_reconstructor.hpp +++ b/src/cudaSirecon/SIM_reconstructor.hpp @@ -105,9 +105,7 @@ class SIM_Reconstructor { int m_zoffset; po::options_description m_progopts; po::variables_map m_varsmap; - #ifdef MRC IW_MRC_HEADER m_in_out_header; - #endif bool m_OTFfile_valid; int m_argc; diff --git a/src/cudaSirecon/cudaSirecon.cpp b/src/cudaSirecon/cudaSirecon.cpp index b47e905..c96de5b 100644 --- a/src/cudaSirecon/cudaSirecon.cpp +++ b/src/cudaSirecon/cudaSirecon.cpp @@ -1,13 +1,10 @@ -#include "cudaSirecon.h" #include "cudaSireconImpl.h" +#include "mrc.h" #include "SIM_reconstructor.hpp" #include "cudasireconConfig.h" #include -#ifdef MRC -#include "mrc.h" -#endif std::string version_number = std::to_string(cudasirecon_VERSION_MAJOR) + "." + std::to_string(cudasirecon_VERSION_MINOR) + "." + @@ -167,7 +164,6 @@ void getbg_and_slope(const char *corrfiles, float *background, int ixyz[3], mxyz[3], pixeltype; /* variables for IMRdHdr call */ float min, max, mean; /* variables for IMRdHdr call */ -#ifdef MRC IW_MRC_HEADER header; if (IMOpen(cstream_no, corrfiles, "ro")) { fprintf(stderr, "File %s does not exist\n", corrfiles); @@ -188,7 +184,6 @@ void getbg_and_slope(const char *corrfiles, float *background, IMRdSec(cstream_no, slope); IMClose(cstream_no); -#endif } void findModulationVectorsAndPhasesForAllDirections( @@ -271,7 +266,6 @@ void findModulationVectorsAndPhasesForAllDirections( // TO-DO std::cout << "No unmixed raw images were saved in TIFF mode\n"; } - #ifdef MRC else if (!params->bTIFF && params->bSaveSeparated) { CPUBuffer tmp((*rawImages)[0].getSize()); for (int phase = 0; phase < params->nphases; ++ phase) { @@ -284,7 +278,6 @@ void findModulationVectorsAndPhasesForAllDirections( } continue; // skip the k0 search and modamp fitting } - #endif /* After separation and FFT, the std::vector rawImages, now referred to as @@ -337,7 +330,6 @@ void findModulationVectorsAndPhasesForAllDirections( true); // ovlp0 shares buffer with tmp0 (hence "true" in the final parameter) ovlp0.save_tiff(params->fileOverlaps.c_str()); } - #ifdef MRC else { CPUBuffer tmp0(data->overlap0.getSize()); data->overlap0.set(&tmp0, 0, tmp0.getSize(), 0); @@ -350,7 +342,6 @@ void findModulationVectorsAndPhasesForAllDirections( IMWrSec(overlaps_stream_no, ol1Ptr + z * imgParams.nx * imgParams.ny); } } - #endif } printf("Initial guess by findk0() of k0[direction %d] = (%f,%f) pixels\n", direction, data->k0[direction].x/dkx, data->k0[direction].y/dky); @@ -429,7 +420,6 @@ void findModulationVectorsAndPhasesForAllDirections( // To-DO std::cout << "No overlaps were saved in TIFF mode"; } - #ifdef MRC else if (!params->bTIFF && order == 1 && params->bSaveOverlaps) {// output the overlaps // output the overlaps CPUBuffer tmp0(data->overlap0.getSize()); @@ -443,7 +433,6 @@ void findModulationVectorsAndPhasesForAllDirections( IMWrSec(overlaps_stream_no, ol1Ptr + z * imgParams.nx * imgParams.ny); } } - #endif } } /* if(searchforvector) ... else ... */ @@ -856,10 +845,6 @@ SIM_Reconstructor::SIM_Reconstructor(int argc, char **argv) // Suppress "unknown field" warnings TIFFSetWarningHandler(NULL); m_imgParams.ntimes = m_all_matching_files.size(); - } else { - #ifndef MRC - throw std::runtime_error("No tiff files found, and program was not compiled with MRC support."); - #endif } m_OTFfile_valid = false; @@ -867,10 +852,8 @@ SIM_Reconstructor::SIM_Reconstructor(int argc, char **argv) // deviceMemoryUsage(); setup(); - #ifdef MRC if (!m_myParams.bTIFF) ::setOutputHeader(m_myParams, m_imgParams, m_in_out_header); - #endif //! Load flat-field correction data bgAndSlope(m_myParams, m_imgParams, &m_reconData); } @@ -1074,7 +1057,6 @@ void SIM_Reconstructor::openFiles() { throw std::runtime_error( "Calling openFiles() but no input or output files were specified"); -#ifdef MRC if (!m_myParams.bTIFF && IMOpen(istream_no, m_myParams.ifiles.c_str(), "ro")) // TIFF files are not opened till loadAndRescaleImage() throw std::runtime_error("No matching input TIFF or MRC files not found"); @@ -1086,19 +1068,16 @@ void SIM_Reconstructor::openFiles() { std::cerr << "File " << m_myParams.ofiles << " can not be created.\n"; throw std::runtime_error("Error creating output file"); } -#endif if (m_myParams.bTIFF) { // will throw CImgIOException if file cannot be opened m_otf_tiff.assign(m_myParams.otffiles.c_str()); m_OTFfile_valid = true; } -#ifdef MRC else if (IMOpen(otfstream_no, m_myParams.otffiles.c_str(), "ro")) throw std::runtime_error("OTF file not found"); else m_OTFfile_valid = true; -#endif } void SIM_Reconstructor::setup(unsigned nx, unsigned ny, unsigned nImages, unsigned nChannels) @@ -1126,10 +1105,8 @@ void SIM_Reconstructor::setup() m_imgParams.nwaves = tiff0.spectrum(); // multi-color not supported yet m_imgParams.nz /= m_imgParams.nwaves * m_myParams.nphases * m_myParams.ndirs; } -#ifdef MRC else ::loadHeader(m_myParams, &m_imgParams, m_in_out_header); -#endif printf("nx_raw=%d, ny=%d, nz=%d\n", m_imgParams.nx_raw, m_imgParams.ny, m_imgParams.nz); @@ -1176,8 +1153,7 @@ void SIM_Reconstructor::setup_common() m_imgParams.nz0, m_imgParams.nwaves); printf("dxy=%f, dz=%f um\n", m_imgParams.dxy, m_imgParams.dz); - #ifdef MRC - // // Initialize headers for intermediate output files if requested + // Initialize headers for intermediate output files if requested if (!m_myParams.bTIFF) { if (m_myParams.bSaveAlignedRaw) { memcpy(&aligned_header, &m_in_out_header, sizeof(m_in_out_header)); @@ -1207,7 +1183,6 @@ void SIM_Reconstructor::setup_common() IMPutHdr(overlaps_stream_no, &overlaps_header); } } - #endif m_myParams.norders = 0; if (m_myParams.norders_output != 0) { @@ -1283,10 +1258,8 @@ int SIM_Reconstructor::processOneVolume() // deviceMemoryUsage(); - #ifdef MRC if (saveIntermediateDataForDebugging(m_myParams)) return 0; - #endif m_reconData.bigbuffer.resize((m_myParams.zoomfact * m_imgParams.nx) * (m_myParams.zoomfact * m_imgParams.ny) * (m_myParams.z_zoom * m_imgParams.nz0) * @@ -1339,7 +1312,6 @@ void SIM_Reconstructor::setFile(int it, int iw) if (m_myParams.bTIFF) { rawImage.assign(m_all_matching_files[it].c_str()); } - #ifdef MRC else { rawImage.assign(m_imgParams.nx_raw, m_imgParams.ny, m_imgParams.nz * m_myParams.nphases * m_myParams.ndirs); for (unsigned sec=0; sec using namespace cimg_library; -#ifdef MRC -#include // MRC file I/O routines -#endif - +#include "dvfile.h" // Block sizes for reduction kernels #define RED_BLOCK_SIZE_X 64 @@ -74,12 +71,10 @@ static const int aligned_stream_no = 10; static const int separated_stream_no = 11; static const int overlaps_stream_no = 12; -#ifdef MRC // static IW_MRC_HEADER header; static IW_MRC_HEADER aligned_header; static IW_MRC_HEADER sep_header; static IW_MRC_HEADER overlaps_header; -#endif struct myExtHeader { float timestamp; @@ -249,18 +244,14 @@ void SetDefaultParams(ReconParams *pParams); * */ // void setup(ReconParams* params, ImageParams* // imgParams, DriftParams* driftParams, ReconData* data); -#ifdef MRC void loadHeader(const ReconParams& params, ImageParams* imgParams, IW_MRC_HEADER &header); -#endif void allocateOTFs(ReconParams *pParams, int sizeOTF, std::vector > & otfs); void allocateImageBuffers(const ReconParams& params, const ImageParams& imgParams, ReconData* reconData); -#ifdef MRC void setOutputHeader(const ReconParams& myParams, const ImageParams& imgParams, IW_MRC_HEADER &header); -#endif void bgAndSlope(const ReconParams& myParams, const ImageParams& imgParams, ReconData* reconData); @@ -289,9 +280,7 @@ void deskewOneSection(CImg<> &rawSection, float* nxp2OutBuff, int z, int nz, // float *background, float backgroundExtra, float *slope, float inscale, // int bUsecorr); -#ifdef MRC int saveIntermediateDataForDebugging(const ReconParams& params); -#endif void matrix_transpose(float* mat, int nRows, int nCols); @@ -316,9 +305,7 @@ void writeResult(int it, int iw, const ReconParams& params, const ImageParams& imgParams, const ReconData& reconData); // This only works for MRC/DV files for now: -#ifdef MRC void saveCommandLineToHeader(int argc, char **argv, IW_MRC_HEADER &header); -#endif void dumpBands(std::vector* bands, int nx, int ny, int nz0); diff --git a/src/cudaSirecon/mrc.h b/src/cudaSirecon/mrc.h index b43e1e3..1b091f0 100644 --- a/src/cudaSirecon/mrc.h +++ b/src/cudaSirecon/mrc.h @@ -1,18 +1,6 @@ // MRC-specific functionality #include "cudaSirecon.h" -#include "IM.h" -#include "IMInclude.h" -#include "IWApiConstants.h" -#if defined(_MSC_VER) && (_MSC_VER >= 1900) -// for IMLIB with vs >2015 -// https://stackoverflow.com/questions/30412951/unresolved-external-symbol-imp-fprintf-and-imp-iob-func-sdl2 -FILE _iob[] = {*stdin, *stdout, *stderr}; -extern "C" FILE *__cdecl __iob_func(void) -{ - return _iob; -} -#endif // MRC header parser void loadHeader(const ReconParams& params, ImageParams* imgParams, diff --git a/src/dvfile.h b/src/dvfile.h new file mode 100644 index 0000000..b661f0a --- /dev/null +++ b/src/dvfile.h @@ -0,0 +1,547 @@ +#ifndef DVFILE_H +#define DVFILE_H + +#include +#include +#include +#include +#include +#include + +// Data Types +#define IW_AS_IS -1 +#define IW_BYTE 0 +#define IW_SHORT 1 +#define IW_FLOAT 2 +#define IW_COMPLEX_SHORT 3 +#define IW_COMPLEX 4 +#define IW_EMTOM 5 +#define IW_USHORT 6 +#define IW_LONG 7 + +// image sequence definitions +#define ZTW_SEQUENCE 0 /* "non-interleaved", by defintion */ +#define WZT_SEQUENCE 1 /* "interleaved", from R3D and others */ +#define ZWT_SEQUENCE 2 /* new sequence. Unsupported as of 11/97 */ + +// Possible values for the file type +#define IM_NORMAL_IMAGES 0 +#define IM_TILT_SERIES 1 +#define IM_STEREO_TILT_SERIES 2 +#define IM_AVERAGED_IMAGES 3 +#define IM_AVERAGED_STEREO_PAIRS 4 +#define IM_EM_TILT_SERIES 5 +#define IM_MULTIPOSITION 20 +#define IM_PUPIL_FUNCTION 8000 + +static const std::unordered_map pixelTypeSizes = { + {IW_BYTE, sizeof(uint8_t)}, {IW_SHORT, sizeof(int16_t)}, + {IW_FLOAT, sizeof(float)}, {IW_COMPLEX_SHORT, 2 * sizeof(int16_t)}, + {IW_COMPLEX, 2 * sizeof(float)}, {IW_EMTOM, sizeof(int16_t)}, + {IW_USHORT, sizeof(uint16_t)}, {IW_LONG, sizeof(int32_t)}}; + +typedef struct IW_MRC_Header { + int32_t nx, ny, nz; // nz=num_planes*num_waves*num_times + int32_t mode; // data type + int32_t nxst, nyst, nzst; // index of the first col/row/section + int32_t mx, my, mz; // number of intervals in x/y/z + float xlen, ylen, zlen; // pixel spacing for x/y/z + float alpha, beta, gamma; // cell angles + int32_t mapc, mapr, maps; // column/row/section axis + float amin, amax, amean; // min/max/mean intensity + int32_t ispg, inbsym; // space group number, number of bytes in extended header + int16_t nDVID, nblank; // ID value, unused + int32_t ntst; // starting time index (used for time series data) + char ibyte[24]; // 24 bytes of blank space + int16_t nint, nreal; // number of integers/floats in extended header per section + int16_t nres, nzfact; // number of sub-resolution data sets, reduction quotient for z axis + float min2, max2, min3, max3, min4, max4; // min/max intensity for 2nd, 3rd, 4th wavelengths + int16_t file_type, lens, n1, n2, v1, v2; // file type, lens ID, n1, n2, v1, v2 + float min5, max5; // min/max intensity for 5th wavelength + int16_t num_times; // number of time points + int16_t interleaved; // (0 = ZTW, 1 = WZT, 2 = ZWT) + float tilt_x, tilt_y, tilt_z; // x/y/z axis tilt angles + int16_t num_waves, iwav1, iwav2, iwav3, iwav4, iwav5; // number & values of wavelengths + float zorig, xorig, yorig; // z/x/y origin + int32_t nlab; // number of titles + char label[800]; + + std::string sequence_order() const { + switch (interleaved) { + case ZTW_SEQUENCE: return "ZTW"; + case WZT_SEQUENCE: return "WZT"; + case ZWT_SEQUENCE: return "ZWT"; + default: return "ZTW"; + } + } + + int num_planes() const { return nz / (num_waves ? num_waves : 1) / (num_times ? num_times : 1); } + + std::string image_type() const { + switch (file_type) { + case IM_NORMAL_IMAGES: return "NORMAL"; + case 100: return "NORMAL"; + case IM_TILT_SERIES: return "TILT_SERIES"; + case IM_STEREO_TILT_SERIES: return "STEREO_TILT_SERIES"; + case IM_AVERAGED_IMAGES: return "AVERAGED_IMAGES"; + case IM_AVERAGED_STEREO_PAIRS: return "AVERAGED_STEREO_PAIRS"; + case IM_EM_TILT_SERIES: return "EM_TILT_SERIES"; + case IM_MULTIPOSITION: return "MULTIPOSITION"; + case IM_PUPIL_FUNCTION: return "PUPIL_FUNCTION"; + default: return "UNKNOWN"; + } + } + + void print() { + std::cout << "Header:" << std::endl; + std::cout << " Dimensions: " << ny << "x" << nx << "x" << num_planes() << std::endl; + std::cout << " Number of wavelengths: " << num_waves << std::endl; + std::cout << " Number of time points: " << num_times << std::endl; + std::cout << " Pixel size: " << mode << std::endl; + std::cout << " Pixel spacing: " << xlen << "x" << ylen << "x" << zlen << std::endl; + std::cout << " mxyz: " << mx << "x" << my << "x" << mz << std::endl; + std::cout << " Min/Max/Mean: " << amin << "/" << amax << "/" << amean << std::endl; + std::cout << " Image type: " << image_type() << std::endl; + std::cout << " Sequence order: " << sequence_order() << std::endl; + } + +} IW_MRC_HEADER, *IW_MRC_HEADER_PTR; + +class DVFile { + private: + std::unique_ptr _file; + std::string _path; + bool _big_endian; + IW_MRC_Header hdr; + bool closed = true; + + // Private default constructor + DVFile() = default; + + public: + // Constructor for opening an existing file + DVFile(const std::string& path) { + _path = path; + _file = std::make_unique(path, std::ios::in | std::ios::out | std::ios::binary); + if (!_file->is_open()) { + throw std::runtime_error("Failed to open file"); + } + + // Determine byte order + _file->seekg(24 * 4); + char dvid[2]; + _file->read(dvid, 2); + if (dvid[0] == (char)0xA0 && dvid[1] == (char)0xC0) { + _big_endian = false; + } else if (dvid[0] == (char)0xC0 && dvid[1] == (char)0xA0) { + _big_endian = true; + } else { + throw std::runtime_error(path + " is not a recognized DV file."); + } + + // Read header + _file->seekg(0); + _file->read(reinterpret_cast(&hdr), sizeof(IW_MRC_Header)); + closed = false; + } + + static std::unique_ptr createNew(const std::string& path) { + std::unique_ptr dvfile(new DVFile()); // Use the private default constructor + dvfile->_path = path; + dvfile->_file = std::make_unique( + path, std::ios::in | std::ios::out | std::ios::binary | std::ios::trunc); + + if (!dvfile->_file->is_open()) { + throw std::runtime_error("Failed to create file"); + } + + dvfile->closed = false; + return dvfile; + } + + // Method to read the extended header + void readExtendedHeader() { + if (hdr.inbsym <= 0) { + throw std::runtime_error("No extended header present."); + } + + // Calculate the size of the extended header + int num_sections = hdr.nz; + size_t extended_header_size = (hdr.nint + hdr.nreal) * num_sections * sizeof(int32_t); + + // Allocate buffer for the extended header + size_t bufferSize = extended_header_size / sizeof(int32_t); + std::vector extended_header(bufferSize); + + // Read the extended header + _file->seekg(sizeof(IW_MRC_Header)); + _file->read(reinterpret_cast(extended_header.data()), extended_header_size); + + if (!_file->good()) { + throw std::runtime_error("Failed to read the extended header."); + } + + // Optionally, process the extended header based on your specific format + } + + /** + * @brief Calculate the offset of a section in the file. + * + * 0 (C/C++ macro is ZTW_SEQUENCE) + * This is the normal arrangement for processed images. Sometimes referred to as + * "non-interleaved". is = iz + nz * (it + nt * iw) Reversing that gives iz = is modulo nz, iw = + * is / (nz * nt), and it = (is / nz) modulo nt. 1 (C/C++ macro is WZT_SEQUENCE) This is the + * typical arrangement for images acquired from a microscope. Sometimes referred to as + * "interleaved" is = iw + nw * (iz + nz * it). Reversing that gives iz = (is / nw) modulo nz, iw + * = is modulo nw, and it = is / ( nw * nz). 2 (C/C++ macro is ZWT_SEQUENCE) Although not widely + * used, ZWT will find uses with certain processing algorithms. is = iz + nz * (iw + nw * it). + * Reversing that gives iz = is modulo nz, iw = (is / nz) modulo nw, and it = is / (nz * nw). + */ + int sectionOffset(int iz, int iw, int it) { + _validateZWT(iz, iw, it); + int is; + int nz = hdr.num_planes(); + int nt = hdr.num_times; + int nw = hdr.num_waves; + switch (hdr.interleaved) { + case WZT_SEQUENCE: is = iw + nw * (iz + nz * it); break; + case ZWT_SEQUENCE: is = iz + nz * (iw + nw * it); break; + case ZTW_SEQUENCE: + default: is = iz + nz * (it + nt * iw); break; + } + return is; + } + + // Set the read/write file pointer to the beginning of a section + void setCurrentZWT(int z, int w, int t) { + int offset = sizeof(IW_MRC_Header) + hdr.inbsym + sectionOffset(z, w, t) * _sectionSize(); + _file->seekg(offset); + _file->seekp(offset); + } + + void readSec(void* array) { + if (closed) { + throw std::runtime_error("Cannot read from closed file. Please reopen with .open()"); + } + _file->read(reinterpret_cast(array), _sectionSize()); + } + + void readSec(void* array, int t, int w, int z) { + setCurrentZWT(z, w, t); + readSec(array); + } + + // call setCurrentZWT to set the z/w/t plane before calling this + void writeSection(const void* array) { + if (closed) { + throw std::runtime_error("Cannot write to closed file. Please reopen with .open()"); + } + _file->write(reinterpret_cast(array), _sectionSize()); + } + + size_t getPixelSize() { return pixelTypeSizes.at(hdr.mode); } + + void open() { + if (closed) { + _file->open(_path, std::ios::binary); + if (!_file->is_open()) { + throw std::runtime_error("Failed to open file"); + } + closed = false; + } + } + + void close() { + if (!closed) { + _file->close(); + closed = true; + } + } + + ~DVFile() { close(); } + + std::string getPath() const { return _path; } + + IW_MRC_Header getHeader() const { return hdr; } + + void putHeader(const IW_MRC_Header& header) { + if (closed) { + throw std::runtime_error("Cannot write to closed file. Please reopen with .open()"); + } + _file->seekp(0); + _file->write(reinterpret_cast(&header), sizeof(IW_MRC_Header)); + hdr = header; + } + + bool isClosed() const { return closed; } + + private: + // Return the size of a frame/section in bytes + size_t _sectionSize() { return hdr.ny * hdr.nx * getPixelSize(); } + + void _validateZWT(int z, int w, int t) { + if (t >= hdr.num_times) { + throw std::runtime_error("Time index out of range"); + } + if (w >= hdr.num_waves) { + throw std::runtime_error("Wavelength index out of range"); + } + if (z >= hdr.num_planes()) { + throw std::runtime_error("Section index out of range"); + } + } +}; + +////////////////////////////////////////////////////////////////////////////// +// IVE API +////////////////////////////////////////////////////////////////////////////// + +// stream -> DVFile +static std::map> dvfile_map; + +inline DVFile& getDVFile(int istream) { + auto it = dvfile_map.find(istream); + if (it == dvfile_map.end()) { + throw std::runtime_error("Stream not found: " + std::to_string(istream)); + } + return *(it->second); +} + +/** + * @brief Open an image file and attach it to a stream. + * + * @param istream The input stream to be used for the operation. + * @param name The name of the file to be opened. + * @param attrib The file mode. + * - "ro" Opens an existing file or image window for reading only. + * - "new" Creates a file or image window and opens it for reading and writing. + * @return int 0 if successful and 1 if not. + */ +inline int IMOpen(int istream, const char* name, const char* attrib) { + // Check if the stream identifier is already in use and close it if necessary + if (dvfile_map.find(istream) != dvfile_map.end()) { + dvfile_map[istream]->close(); + dvfile_map.erase(istream); + std::cerr << "Warning: Reusing stream identifier " << istream << ". Previous stream closed." + << std::endl; + } + + if (std::string(attrib) == "ro") { + try { + dvfile_map[istream] = std::make_unique(name); + } catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << std::endl; + return 1; // Return non-zero to indicate failure + } + } else if (std::string(attrib) == "new") { + try { + dvfile_map[istream] = DVFile::createNew(name); + } catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << std::endl; + return 1; // Return non-zero to indicate failure + } + } else { + std::cerr << "Unknown file mode: " << attrib << std::endl; + return 1; // Return non-zero to indicate failure + } + return 0; // Return 0 to indicate success +} + +inline void IMClose(int istream) { + // call destructor of DVFile + dvfile_map.erase(istream); +} + +inline void IMGetHdr(int istream, IW_MRC_HEADER* header) { + *header = getDVFile(istream).getHeader(); +} + +inline void IMRdHdr(int istream, int ixyz[3], int mxyz[3], int* imode, float* min, float* max, + float* mean) { + IW_MRC_HEADER header; + IMGetHdr(istream, &header); + ixyz[0] = header.nx; + ixyz[1] = header.ny; + ixyz[2] = header.nz; + mxyz[0] = header.mx; + mxyz[1] = header.my; + mxyz[2] = header.mz; + *imode = header.mode; + *min = header.amin; + *max = header.amax; + *mean = header.amean; +} + +/** + * @brief Set the image conversion mode during read/write operations from image storage. + * + * By default in IVE, images that are read from image storage are converted to + * 4-byte floating-point data. Similarly, when images are written to image + * storage they are converted to the data type indicated by the image data type + * associated with the corresponding stream (see IMAlMode). The default in IVE + * is ConversionFlag=TRUE. + * We, however, don't ever convert the data type of the image data. So for now, + * this is a no-op. + * + * @param istream The input stream to be used for the operation. + * @param flag The flag indicating the type of operation to be performed. + */ +inline void IMAlCon(int istream, int flag) { + // if flag is 1, warn: + if (flag == 1) { + std::cerr << "Warning: IMAlCon is not implemented. ConversionFlag=TRUE is not supported." + << std::endl; + } +} + +/** + * @brief Change the image titles. + * + * @param istream The input stream to be used for the operation. + * @param titles The titles to be changed. Contains at least NumTitles title strings, each of + * which must contain exactly 80 characters. + * @param num_titles The number of titles to be changed. + */ +inline void IMAlLab(int istream, const char* labels, int nl) { + std::cerr << "Warning: IMAlLab is not implemented." << std::endl; +} + +/** + * @brief Enable or disable printing to standard output ("stdout"). + * + * Certain IM functions will print information to stdout if Format=TRUE, which is the default. To + * disable printing, set flag=FALSE. + * + * @param flag The flag indicating the type of operation to be performed. + */ +inline void IMAlPrt(int flag) { + if (flag == 1) { + std::cerr << "Warning: IMAlPrt is not implemented." << std::endl; + } +} + +/** + * @brief Position the read/write point at a particular Z, W, T section. + * + * If the stream points to a scratch window, IMPosnZWT can only change the destination wavelength. + * + * @param istream The input stream to be used for the operation. + * @param iz The Z section number. + * @param iw The wavelength number. + * @param it The time-point number. + * @return int 0 if successful and 1 if not. + */ +inline int IMPosnZWT(int istream, int iz, int iw, int it) { + DVFile& dvfile = getDVFile(istream); + + try { + dvfile.setCurrentZWT(iz, iw, it); + return 0; + } catch (const std::exception& e) { + std::cerr << e.what() << '\n'; + return 1; + } +} + +/** + * @brief Read the next section. + * + * Reads the next section into ImgBuffer and advances the file pointer to the + * section after that. The results are undefined if ImgBuffer does not have at + * least nx * ny elements or the file pointer does not point to the beginning of + * a section. + * In most cases, ImgBuffer will contain floating-point data. When + * image conversion is off, however, the data type of ImgBuffer should + * correspond to whatever data type is actually stored. For example, if the + * image data are stored as 16-bit integers, then ImgBuffer should point to a + * 16-bit buffer. See IMAlCon. + * + * @param istream The input stream to be used for the operation. + * @param ImgBuffer The image buffer to store the data. + */ +inline void IMRdSec(int istream, void* ImgBuffer) { + try { + getDVFile(istream).readSec(ImgBuffer); + } catch (const std::runtime_error& e) { + std::cerr << "Error reading section: " << e.what() << std::endl; + // Handle the error appropriately, e.g., by rethrowing or returning an error code + throw; // Rethrow the exception to propagate it further + } +} + +inline void IMWrSec(int istream, const void* array) { + try { + getDVFile(istream).writeSection(array); + } catch (const std::runtime_error& e) { + std::cerr << "Error writing section: " << e.what() << std::endl; + // Handle the error appropriately, e.g., by rethrowing or returning an error code + throw; // Rethrow the exception to propagate it further + } +} + +/** + * @brief Put an entire header into a stream. + * + * Header should point to a memory location that contains a complete image header. + * + * @param istream The input stream to be used for the operation. + * @param header The header to be saved in the stream. + */ +inline void IMPutHdr(int istream, const IW_MRC_HEADER* header) { + getDVFile(istream).putHeader(*header); +} + +/** + * @brief Write the image header to the storage device. + * + * Write image header associated with StreamNum to the storage device. Use this + * function to save the results of all IMAl functions. Header modifications are + * not saved until IMWrHdr is used! The contents of Title will be saved in the + * header according to the method indicated by ntflag. The minimum, maximum, and + * mean intensity - Min, Max, and Mean, respectively - of wavelength 0 are also + * saved in the header every time IMWrHdr is used. + * + * @param istream The input stream to be used for the operation. + * @param title The title to be saved in the header. + * @param ntflag The method to be used to save the title. + * - 0: use Title as the only title + * - 1: add Title to the end of the list + */ +inline void IMWrHdr(int istream, const char title[80], int ntflag, float dmin, float dmax, + float dmean) { + DVFile& dvfile = getDVFile(istream); + IW_MRC_HEADER header = dvfile.getHeader(); + header.amin = dmin; + header.amax = dmax; + header.amean = dmean; + if (ntflag == 0) { + // use Title as the only title + strncpy(header.label, title, 80); + } else if (ntflag == 1) { + // FIXME this is wrong. + // Append title to the end of the list + std::string new_title = title; + new_title += " "; + new_title += header.label; + strncpy(header.label, new_title.c_str(), 80); + } else { + throw std::runtime_error("Invalid ntflag: " + std::to_string(ntflag)); + } + + dvfile.putHeader(header); +} + +/** + * @brief Return extended header values for a particular Z section, wavelength, + * and time-point. + * + * The integer and floating-point values for the requested Z section (ZSecNum), + * wavelength (WaveNum), and time-point (TimeNum) are returned in IntValues and + * FloatValues, respectively. + * + */ +inline void IMRtExHdrZWT(int istream, int iz, int iw, int it, int ival[], float rval[]) { + std::cerr << "Warning: IMRtExHdrZWT is not implemented." << std::endl; +} + +#endif // DVFILE_H \ No newline at end of file diff --git a/src/otf/CMakeLists.txt b/src/otf/CMakeLists.txt index 09c26ac..039832e 100644 --- a/src/otf/CMakeLists.txt +++ b/src/otf/CMakeLists.txt @@ -1,10 +1,6 @@ include_directories("${CMAKE_SOURCE_DIR}/cudaSirecon") # for CImg.h add_executable(makeotf radialft.cpp) -if(BUILD_MRC OR BUILD_OTF_VIEWER) - target_link_libraries(makeotf ${IMLIB} ${IVELIB}) -endif() - # ######## START findFFTW: https://github.com/egpbos/findFFTW ###### configure_file(downloadFindFFTW.cmake.in findFFTW-download/CMakeLists.txt) @@ -52,7 +48,7 @@ if(BUILD_OTF_VIEWER) include_directories(${PRIISM_INCLUDE_PATH}) link_directories(${PRIISM_LIB_PATH}) add_executable(otfviewer otfviewer.cpp) - target_link_libraries(otfviewer ${IMLIB} ${IVELIB} ${TIFF_LIBRARIES}) + target_link_libraries(otfviewer ${TIFF_LIBRARIES}) if(WIN32) if(${MSVC_VERSION} GREATER 1800) diff --git a/src/otf/otfviewer.cpp b/src/otf/otfviewer.cpp index d4a9920..83dc609 100644 --- a/src/otf/otfviewer.cpp +++ b/src/otf/otfviewer.cpp @@ -9,18 +9,6 @@ using namespace cimg_library; #include #include -#if defined(_MSC_VER) && (_MSC_VER >= 1900) -// for IMLIB with vs >2015 -// https://stackoverflow.com/questions/30412951/unresolved-external-symbol-imp-fprintf-and-imp-iob-func-sdl2 -FILE _iob[] = {*stdin, *stdout, *stderr}; - -extern "C" FILE * __iob_func(void) -{ - return _iob; -} - -#endif - int main(int argc, char *argv[]) { @@ -64,7 +52,7 @@ int main(int argc, char *argv[]) int istream_no = 1; // Suppress IM header printout; somehow only in this mode would the rest of // IM calls go well on Windows. - IMAlPrt(0); + // IMAlPrt(0); // not needed with our dvfile implementation if (IMOpen(istream_no, argv[optind], "ro")) { std::cerr << "File " << argv[optind] << " cannot be opened.\n"; diff --git a/src/otf/radialft.cpp b/src/otf/radialft.cpp index ac6730f..3594b9f 100644 --- a/src/otf/radialft.cpp +++ b/src/otf/radialft.cpp @@ -22,9 +22,7 @@ #include using namespace cimg_library; -#ifdef MRC -#include -#endif +#include "dvfile.h" #include #include @@ -60,24 +58,11 @@ void rescale(std::complex *otfkxkz, int order, int nx, int nz, float *sca void outputdata(std::string &tiffile, std::vector *> &bands, int norders, int nx, int ny, int nz, float dkr, float dkz, int five_bands); -#ifdef MRC void outputdata(int ostream_no, IW_MRC_HEADER *header, std::vector *> &bands, int norders, int nx, int ny, int nz, float dkr, float dkz, int five_bands); void mrc_file_write(float *buffer, int nx, int ny, int nz, float rlen, float zlen, int mode, int iwave, const char *files); -#endif int commandline(int argc, char *argv[], int * twolens, int *rescale, float *beaddiam, float *k0angle, float *linespacing, int *five_bands, int *nphases, int *interpkr, int *leavekz, int *do_compen, int *I2M_inc, std::string &I2Mfiles, float *background, int *bBgInExtHdr, int *order0gen, std::string &order0files, int *conjugate, float *na, float *nimm, int *ifixz, int *ifixr, unsigned *wavelength, float *dr, float *dz, int *bCoherentBSIM, int *bForcedPIshift, std::string &ifiles, std::string &ofiles, int *ifilein, int *ofilein); -#if defined(_MSC_VER) && (_MSC_VER >= 1900) -// for IMLIB with vs >2015 -// https://stackoverflow.com/questions/30412951/unresolved-external-symbol-imp-fprintf-and-imp-iob-func-sdl2 -FILE _iob[] = {*stdin, *stdout, *stderr}; - -extern "C" FILE * __iob_func(void) -{ - return _iob; -} - -#endif int main(int argc, char **argv) { @@ -95,11 +80,7 @@ int main(int argc, char **argv) int bCoherentBSIM = 0; /* In coherent Bessel-SIM, do fixorigin() differently? */ int bForcedPIshift=0; - #ifdef MRC IW_MRC_HEADER header, otfheader; - IMAlPrt(0); - #endif - interpkr[0] = 0; interpkr[1] = 0; @@ -130,7 +111,6 @@ int main(int argc, char **argv) PSFtiff.assign(ifiles.c_str()); bTIFF = true; } - #ifdef MRC else { std::cout << "Not a TIFF file; now try MRC\n"; if (IMOpen(istream_no, ifiles.c_str(), "ro")) { @@ -138,18 +118,15 @@ int main(int argc, char **argv) return -1; } } - #endif if (!ofilein) { std::cout << "Output OTF file name: "; getline(std::cin, ofiles); } - #ifdef MRC if (!bTIFF) if (IMOpen(ostream_no, ofiles.c_str(), "new")) { std::cerr << "File " << ofiles << " cannot be created.\n"; return -1; } - #endif unsigned nx, nxExtra, ny, nz, nxy; if (bTIFF) { @@ -162,7 +139,6 @@ int main(int argc, char **argv) dz = user_dz; } - #ifdef MRC else { int ixyz[3], mxyz[3], pixeltype; float min, max, mean; @@ -176,7 +152,6 @@ int main(int argc, char **argv) dz = header.zlen; wavelength = header.iwav1; } - #endif if (!I2M_inc) nz /= nphases; @@ -219,12 +194,10 @@ int main(int argc, char **argv) for(int phase=0; phase