opencliper.lpi.tel.uva.es Git - OpenCLIPER/blob - src/KData.cpp

   1 /* Copyright (C) 2018 Federico Simmross Wattenberg,
   2  *                    Manuel Rodríguez Cayetano,
   3  *                    Javier Royuela del Val,
   4  *                    Elena Martín González,
   5  *                    Elisa Moya Sáez,
   6  *                    Marcos Martín Fernández and
   7  *                    Carlos Alberola López
   8  *
   9  * This file is part of OpenCLIPER.
  10  *
  11  * OpenCLIPER is free software; you can redistribute it and/or modify
  12  * it under the terms of the GNU General Public License as published by
  13  * the Free Software Foundation; version 3 of the License.
  14  *
  15  * OpenCLIPER is distributed in the hope that it will be useful, but
  16  * WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU General Public License
  21  * along with OpenCLIPER; If not, see <http://www.gnu.org/licenses/>.
  22  *
  23  *
  24  *  Contact:
  25  *
  26  *  Federico Simmross Wattenberg
  27  *  E.T.S.I. Telecomunicación
  28  *  Universidad de Valladolid
  29  *  Paseo de Belén 15
  30  *  47011 Valladolid, Spain.
  31  *  fedsim@tel.uva.es
  32  */
  33 /*
  34  * KData.cpp
  35  *
  36  *  Created on: 28 de oct. de 2016
  37  *      Author: manrod
  38  */
  39
  40 #include <OpenCLIPER/KData.hpp>
  41 #define ERRORPREFIX "OpenCLIPER::KData::"
  42
  43 namespace OpenCLIPER {
  44
  45 /**
  46  * @brief Constructor that creates an empty KData object.
  47  * @param[in] automaticStoreOnDevice sets the automatic copy of data from host memory to device memory feature if true
  48  * (default value: true)
  49  */
  50 KData::KData(bool automaticStoreOnDevice): Data(automaticStoreOnDevice) {
  51     // TODO Auto-generated constructor stub
  52 }
  53
  54 /**
  55  * @brief Constructor that creates a KData object from a vector of NDArrays.
  56  *
  57  * All parameters of type *& (reference to pointer) have move semantics: ownership of memory is moved from parameter to this
  58  * object, and parameter value is set to nullptr after method completion.
  59  * @param[in,out] pSensitivityMapsData pointer to a SensitivityMaps object (group of sensitivity maps for coils)
  60  * @param[in,out] pData pointer to vector of pointers to NDArray objects
  61  * @param[in,out] pCoord
  62  * @param[in] nCoils total number of coils available
  63  * @param[in] usedCoils number of coils used
  64  * @param[in,out] pDynDims pointer to vector of temporal dimensions
  65  * @param[in] trajectory trajectory used
  66  * @param[in,out] pDcf
  67  * @param[in,out] pDeltaK
  68  * @param[in] automaticStoreOnDevice sets the automatic copy of data from host memory to device memory feature if true
  69  * (default value: true)
  70  */
  71 KData::KData(SensitivityMapsData*& pSensitivityMapsData, vector<NDArray*>*& pData, vector<realType>*& pCoord, numCoilsType nCoils, set usedCoils,
  72              vector<dimIndexType>*& pDynDims, enum TrajType trajectory, vector<realType>*& pDcf, vector<realType>*& pDeltaK,
  73              bool automaticStoreOnDevice): Data(pData, automaticStoreOnDevice) {
  74     // TODO Auto-generated constructor stub
  75     /* Falta comprobar:
  76      * - tamaño coord sea nDims
  77      * - tamaño dcf sea nDims
  78      * - tamaño deltaK sea nDims
  79      * - tamaño usedCoils debe ser nCils
  80      * Problema: nDims hay que leerlo de un objeto de la clase NDArray,
  81      * pero KData tiene referencias a varios objetos de tipo NDArray
  82      * (vector data que pueden tener valores distintos de nDims)
  83      */
  84     this->pSensitivityMapsData.reset(pSensitivityMapsData);
  85     // delete(this->pSensitivityMapsData);
  86     // this->pSensitivityMapsData = pSensitivityMapsData;
  87     pSensitivityMapsData = nullptr;
  88     this->pCoord.reset(pCoord);
  89     pCoord = nullptr;
  90     this->nCoils = nCoils;
  91     this->usedCoils = usedCoils;
  92     this->trajectory = trajectory;
  93     this->pDcf.reset(pDcf);
  94     pDcf = nullptr;
  95     this->pDeltaK.reset(pDeltaK);
  96     pDeltaK = nullptr;
  97     this->pDynDims.reset(pDynDims);
  98     pDynDims = nullptr;
  99     dimIndexType totalNumberOfDynDims = getDynDimsTotalSize();
 100     if ((totalNumberOfDynDims * nCoils) != pData->size()) {
 101         throw std::invalid_argument(
 102             "sum of dynamic dimensions (DynDims) * number of coils must be equal to kData vector size");
 103     }
 104     internalSetData(pData);
 105 }
 106
 107 /**
 108  * @brief Constructor that creates an empty KData object but with spatial and temporal dimensions set.
 109  *
 110  * All parameters of type *& (reference to pointer) have move semantics: ownership of memory is moved from parameter to this
 111  * object, and parameter value is set to nullptr after method completion.
 112  * @param[in,out] pArraysDims pointer to a vector of pointers to vectors with the spatial dimensions of every image
 113  * @param[in] nCoils total number of coils available
 114  * @param[in,out] pDynDims pointer to vector of temporal dimensions
 115  * @param[in] automaticStoreOnDevice sets the automatic copy of data from host memory to device memory feature if true
 116  * (default value: true)
 117  */
 118 KData::KData(std::vector<std::vector<dimIndexType>*>*& pArraysDims, numCoilsType nCoils, vector<dimIndexType>*& pDynDims, bool automaticStoreOnDevice) :
 119     Data(pArraysDims, pDynDims, automaticStoreOnDevice) {
 120     // All dimensions for every NDArray must be valid
 121     //checkValidSpatialDimensions();
 122     this->nCoils = nCoils;
 123 }
 124
 125 /**
 126  * @brief Constructor that creates a KData object with data (optional), spatial and temporal dimensions got from
 127  * another KData object.
 128  * @param[in] sourceData KData object source of spatial and temporal dimensions
 129  * @param[in] copyData if true also data (not only dimensions) are copied from sourceData object to this object
 130  * @param[in] copySensitivityMaps if true also sensitivity maps (not only dimensions) are copied from sourceData object to this object
 131  * @param[in] copySamplingMasks if true also sampling masks (not only dimensions) are copied from sourceData object to this object
 132  * @param[in] automaticStoreOnDevice sets the automatic copy of data from host memory to device memory feature if true
 133  * (default value: true)
 134  */
 135 KData::KData(const KData* sourceData, bool copyData, bool copySensitivityMaps, bool copySamplingMasks, bool automaticStoreOnDevice) :
 136     Data ((const Data*) sourceData, copyData, automaticStoreOnDevice) {
 137     vector<realType>* pLocalCoord;
 138     if (sourceData->getCoord() == nullptr) {
 139         pLocalCoord = nullptr;
 140     } else {
 141         pLocalCoord = new vector<realType>(*(sourceData->getCoord()));
 142     }
 143     setCoord(pLocalCoord);
 144     setNCoils(sourceData->getNCoils());
 145     setUsedCoils(sourceData->getUsedCoils());
 146     setTrajectory(sourceData->getTrajectory());
 147     vector<realType>* pLocalDcf;
 148     if (sourceData->getDcf() == nullptr) {
 149         pLocalDcf = nullptr;
 150     } else {
 151         pLocalDcf = new vector<realType>(*(sourceData->getDcf()));
 152     }
 153     setDcf(pLocalDcf);
 154     vector<realType>* pLocalDeltaK;
 155     if (sourceData->getDeltaK() == nullptr) {
 156         pLocalDeltaK = nullptr;
 157     } else {
 158         pLocalDeltaK = new vector<realType>(*(sourceData->getDeltaK()));
 159     }
 160     setDeltaK(pLocalDeltaK);
 161     if (copySensitivityMaps == true) {
 162         SensitivityMapsData* sensitivityMapsData = new SensitivityMapsData(sourceData->getSensitivityMapsData(), true);
 163         setSensitivityMapsData(sensitivityMapsData);
 164     }
 165     if (copySamplingMasks == true) {
 166         SamplingMasksData* samplingMasksData = new SamplingMasksData(sourceData->getSamplingMasksData(), true);
 167         setSamplingMasksData(samplingMasksData);
 168     }
 169 }
 170
 171 /**
 172  * @brief Constructor tha loads data of a group of files in raw format (see @ref loadRawHostData) to NDArray objects
 173  * (every file contains one image and is stored into a NDArray object).
 174  *
 175  * Files of this group have names with the format
 176  * \<fileNamePrefix\>\<dims\>\<coilsFileNameSuffix\>\<coilNumber\>\<framesFileNameSuffix\>\<frameNumber\>\<fileNameExtension\>,
 177  * where \<dims\> is a string with the format \<width\>x\<height\>x\<depth\>, \<coilNumber\> is the identifier of the coil
 178  * used for image adquisition and \<frameNumber\> is the identifier of the time of image adquisition.
 179  * @param[in] dataFileNamePrefix fixed part of the file name before the variable part
 180  * @param[in,out] pArraysDims pointer to a vector of pointers to vectors with the dimensions of every image
 181  * @param[in] numCoils number of images per time frame (number of coils)
 182  * @param[in,out] pDynDims pointer to a vector with the temporal dimensions of every image of the sequence
 183  * @param[in] dataToLoad bit mask storing extra KData fields to be loaded (LOADSENSITIVITYMAPS for loading sensitivity maps and
 184  * LOADSAMPLINGMASKS for loading sampling masks)
 185  * @param[in] otherFieldsFileNamePrefixes name prefixes (after dataFileNamePrefix) for sensitivity maps and sampling masks
 186  * @param[in] coilsFileNameSuffix part of the file name before the coil number
 187  * @param[in] framesFileNameSuffix part of the file name before the frame number
 188  * @param[in] fileNameExtension extension for file names
 189  * @param[in] automaticStoreOnDevice sets the automatic copy of data from host memory to device memory feature if true
 190  * (default value: true)
 191  */
 192 KData::KData(const std::string dataFileNamePrefix, vector<vector< dimIndexType >*>*& pArraysDims, numCoilsType numCoils,
 193              vector <dimIndexType>*& pDynDims, uint dataToLoad,
 194              const vector<std::string> otherFieldsFileNamePrefixes, const std::string coilsFileNameSuffix,
 195              const std::string framesFileNameSuffix, const std::string fileNameExtension, bool automaticStoreOnDevice):
 196     Data(automaticStoreOnDevice) {
 197     loadRawHostData(dataFileNamePrefix, otherFieldsFileNamePrefixes, dataToLoad, pArraysDims, numCoils, pDynDims,
 198                     coilsFileNameSuffix, framesFileNameSuffix, fileNameExtension, automaticStoreOnDevice);
 199 }
 200
 201 /**
 202  * @brief Constructor that creates an KData object from a file in matlab format (containing one ore more variables).
 203  * @param[in] fileName name of the data file
 204  * @param[in] variableNames vector with names of matlab variables to be read from file (first element is the name of the
 205  * variable containing k-space image data, second is for the variable containing sensitivity maps data and third is
 206  * for the variable containing sampling masks data).
 207  * @param[in] automaticStoreOnDevice sets the automatic copy of data from host memory to device memory feature if true
 208  * (default value: true)
 209  */
 210 KData::KData(const std::string fileName, const vector<string> variableNames, bool automaticStoreOnDevice) :
 211     Data(automaticStoreOnDevice) {
 212     map<string, matvar_t*>*pMatlabVariablesMap;
 213     pMatlabVariablesMap = Data::readMatlabVariablesFromFile(fileName, variableNames);
 214     matvar_t *pKDataMatVar = nullptr, *pSensitivityMapsDataMatVar = nullptr, *pSamplingMasksDataMatVar = nullptr;
 215     pKDataMatVar = pMatlabVariablesMap->at(variableNames.at(0));
 216     if (variableNames.size() > 1) {
 217         pSensitivityMapsDataMatVar = pMatlabVariablesMap->at(variableNames.at(1));
 218     }
 219     if (variableNames.size() > 2) {
 220         pSamplingMasksDataMatVar = pMatlabVariablesMap->at(variableNames.at(2));
 221     }
 222     dimIndexType numOfSpatialDimensions, numOfNDArraysToBeRead;
 223
 224     // KData number of dimensions = number of spatial dimensions + 1 (coil dimension) + number of temporal dimensions
 225     // SensitivityMaps number of dimensions = number of spatial dimensions + 1 (coil dimension)
 226     // SamplimgMasks number of dimensions = number of spatial dimensions - 1 + number of temporal dimensions
 227
 228     // Special case: if number of coils is 1, sensitivity maps number of dimensions (spatial) is 2 (also for
 229     // KData and sampling masks spatial dimensions)
 230     if (pSensitivityMapsDataMatVar->rank == 2) {
 231          numOfSpatialDimensions = 2;
 232          setNCoils(1);
 233     } else {
 234         numOfSpatialDimensions = pSensitivityMapsDataMatVar->rank - 1;
 235         // Number of coils is value of first dimension after last spatial dimension
 236         // coil dimension is at numOfSpatialDimensions position in matlab dims array
 237         setNCoils(pKDataMatVar->dims[numOfSpatialDimensions]);
 238     }
 239     // Vector of temporal dimensions for Data object
 240     vector<dimIndexType>* pDynDims = new vector<dimIndexType>;
 241
 242     // Loop over temporal dimensions, first temporal dimension comes after number of coils
 243     for (dimIndexType i = numOfSpatialDimensions + 1; i < pKDataMatVar->rank; i++) {
 244         pDynDims->push_back(pKDataMatVar->dims[i]);
 245     }
 246     setDynDims(pDynDims);
 247     // Number of KData NDArrays to be read is number of coils * product of all temporal dimensions
 248     numOfNDArraysToBeRead = getNCoils() * getDynDimsTotalSize();
 249     ((Data*)(this))->loadMatlabHostData(pKDataMatVar, numOfSpatialDimensions, numOfNDArraysToBeRead);
 250
 251     try {
 252         // Create sensitivity maps and store them in class field
 253         SensitivityMapsData* pSensitivityMapsData = new SensitivityMapsData(pSensitivityMapsDataMatVar, automaticStoreOnDevice);
 254         setSensitivityMapsData(pSensitivityMapsData);
 255     } catch (std::invalid_argument &e) {
 256         CERR("Warning: no senstivity maps in KData\n");
 257     }
 258     try {
 259         // Create sampling masks and store them in class field (number of spatial dimensions of sampling mask data is
 260         // number of spatial dimensions of KData - 1
 261         SamplingMasksData* pSamplingMasksData = new SamplingMasksData(pSamplingMasksDataMatVar, automaticStoreOnDevice);
 262         setSamplingMasksData(pSamplingMasksData);
 263     } catch (std::invalid_argument &e) {
 264         CERR("Warning: no sampling masks in KData\n");
 265     }
 266 }
 267
 268 void KData::delApp(std::shared_ptr<CLapp> pCLapp) {
 269     Data::delApp(pCLapp);
 270     if (sensitivityMapsDataHandle != INVALIDDATAHANDLE) {
 271         pCLapp->getData(sensitivityMapsDataHandle)->delApp(pCLapp);
 272     }
 273     if (samplingMasksDataHandle != INVALIDDATAHANDLE) {
 274         pCLapp->getData(samplingMasksDataHandle)->delApp(pCLapp);
 275     }
 276 }
 277
 278 KData::~KData() {
 279 #ifdef KData_DEBUG
 280     CERR("~KData() begins..." << std::endl);
 281 #endif
 282     if (sensitivityMapsDataHandle != INVALIDDATAHANDLE) {
 283         sensitivityMapsDataHandle = INVALIDDATAHANDLE;
 284     }
 285     pSensitivityMapsData = nullptr;
 286
 287     pSensitivityMapsRMS = nullptr;
 288
 289     if (samplingMasksDataHandle != INVALIDDATAHANDLE) {
 290         samplingMasksDataHandle = INVALIDDATAHANDLE;
 291     }
 292     pSamplingMasksData = nullptr;
 293
 294     pCoord.reset(nullptr);
 295     pDcf.reset(nullptr);
 296     pDeltaK.reset(nullptr);
 297 #ifdef KData_DEBUG
 298     CERR("~KData() ends..." << std::endl);
 299 #endif
 300 }
 301
 302 /**
 303  * @brief Gets a pointer to one NDArray object (image data) from the vector of NDArray objects.
 304  * @param[in] dynIndexes vector with values for the temporal indexes used for indexing the NDArray objects vector
 305  * @param[in] coilId number of coil from indexing the NDArray objects vector
 306  * @return pointer to NDArray object at the specified position
 307  */
 308 const NDArray* KData::getDataAtDynPosAndCoilId(const vector<dimIndexType> dynIndexes, numCoilsType coilId) const {
 309     index1DType index1D;
 310     index1D = this->get1DIndexFromDynPos(dynIndexes);
 311     index1D = index1D + coilId * this->getDynDimsTotalSize();
 312     printf ("KData index1D: %d\n", index1D);
 313     return pData->at(index1D).get();
 314 }
 315
 316 /**
 317  * @brief Load data of a group of files to NDArray objects (every file contains one image and is stored into a NDArray object).
 318  *
 319  * Files of this group have names with the format
 320  * \<fileNamePrefix\>\<dims\>\<coilsFileNameSuffix\>\<coilNumber\>\<framesFileNameSuffix\>\<frameNumber\>\<fileNameExtension\>,
 321  * where \<dims\> is a string with the format \<width\>x\<height\>x\<depth\>, \<coilNumber\> is the identifier of the coil used
 322  * for image adquisition and \<frameNumber\> is the identifier of the time of image adquisition.
 323  * @param[in] dataFileNamePrefix fixed part of the file name before the variable part
 324  * @param[in] otherFieldsFileNamePrefixes name prefixes (after dataFileNamePrefix) for sensitivity maps and sampling masks
 325  * @param[in] dataToLoad bit mask storing extra KData fields to be loaded (LOADSENSITIVITYMAPS for loading sensitivity maps and
 326  * LOADSAMPLINGMASKS for loading sampling masks)
 327  * @param[in,out] pArraysDims pointer to a vector of pointers to vectors with the dimensions of every image
 328  * @param[in] numCoils number of images per time frame (number of coils)
 329  * @param[in,out] pDynDims pointer to a vector with the temporal dimensions of every image of the sequence
 330  * @param[in] coilsFileNameSuffix part of the file name before the coil number
 331  * @param[in] framesFileNameSuffix part of the file name before the frame number
 332  * @param[in] fileNameExtension extension for file names
 333  * @param[in] automaticStoreOnDevice sets the automatic copy of data from host memory to device memory feature if true
 334  * (default value: true)
 335  */
 336 void KData::loadRawHostData(const std::string dataFileNamePrefix, const vector<std::string> otherFieldsFileNamePrefixes,
 337                             uint dataToLoad, vector<vector< dimIndexType >*>*& pArraysDims, numCoilsType numCoils,
 338                             vector <dimIndexType>*& pDynDims, const std::string coilsFileNameSuffix,
 339                             const std::string framesFileNameSuffix, const std::string fileNameExtension,
 340                             bool automaticStoreOnDevice) {
 341     index1DType numFrames;
 342     if (pDynDims->size() == 0) {
 343         numFrames = 0;
 344     } else {
 345         numFrames = 1;
 346     }
 347     for (index1DType i = 0; i<pDynDims->size(); i++) {
 348         numFrames = numFrames * pDynDims->at(i);
 349     }
 350
 351     // Sensitivity maps loading
 352     if (dataToLoad & LOADSENSITIVITYMAPS) {
 353         // Spatial dimensions vector for sensitivity maps (all images of same size)
 354         vector<vector< dimIndexType >*>* pArraysDimsSensitivityMaps = new vector<vector< dimIndexType >*>();
 355         vector< dimIndexType >* pDims;
 356         //for (dimIndexType i = 0; i < pArraysDims->size(); i++) {
 357         // Size of pArraysDims == number of NDArrays == numCoils * getDynDimsTotalSize
 358         // spatial dimensions array for sampling masks only have
 359         // numCoil positions
 360         for (dimIndexType i = 0; i < numCoils; i++) {
 361             pDims = new vector<dimIndexType>(*(pArraysDims->at(i)));
 362             pArraysDimsSensitivityMaps->push_back(pDims);
 363         }
 364         // New Data object for storing SensitivityMaps data of input KData
 365         // Load SensitivityMaps data from raw files created from matlab
 366
 367         SensitivityMapsData* pSensitivityMapsData =
 368             new OpenCLIPER::SensitivityMapsData(dataFileNamePrefix + otherFieldsFileNamePrefixes.at(SENSITIVITYMAPSPREFIX), pArraysDimsSensitivityMaps,
 369                                            numCoils, coilsFileNameSuffix, fileNameExtension, automaticStoreOnDevice);
 370         setSensitivityMapsData(pSensitivityMapsData);
 371     }
 372
 373     // Sampling masks loading
 374     if (dataToLoad & LOADSAMPLINGMASKS) {
 375         // Spatial dimensions vector for sensitivity maps (all images of same size)
 376         vector<vector< dimIndexType >*>* pArraysDimsSamplingMasks = new vector<vector< dimIndexType >*>();
 377         vector< dimIndexType >* pDims;
 378         // for (dimIndexType i = 0; i < pArraysDims->size(); i++) {
 379         // Size of pArraysDims == number of NDArrays == numCoils * getDynDimsTotalSize
 380         // spatial dimensions array for sampling masks only have
 381         // getDynDimsTotalSize positions
 382         // Sampling masks has 1 spatial dimension (number of columns) for 2D-images
 383         for (dimIndexType i = 0; i < numFrames; i++) {
 384             pDims = new vector<dimIndexType>();
 385             pDims->push_back(pArraysDims->at(i)->at(0));
 386             pArraysDimsSamplingMasks->push_back(pDims);
 387         }
 388         // New Data object for storing SensitivityMaps data of input KData
 389         // Load SensitivityMaps data from raw files created from matlab
 390
 391         vector<dimIndexType>* pDynDimsSamplingMasks = new vector<dimIndexType>(*(pDynDims));
 392         SamplingMasksData* pSamplingMasksData =
 393             new OpenCLIPER::SamplingMasksData(dataFileNamePrefix + otherFieldsFileNamePrefixes.at(SAMPLINGMASKSPREFIX),
 394                                          pArraysDimsSamplingMasks, pDynDimsSamplingMasks, framesFileNameSuffix,
 395                                          fileNameExtension, automaticStoreOnDevice);
 396         setSamplingMasksData(pSamplingMasksData);
 397     }
 398
 399     // KData data loading
 400     this->setNCoils(numCoils);
 401     stringstream variableSuffixStream;
 402     std::vector<std::string> fileNameSuffixes;
 403     for (dimIndexType frame = 0; frame < numFrames; frame++) {
 404         for (numCoilsType coil = 0; coil < numCoils; coil++) {
 405             variableSuffixStream << coilsFileNameSuffix << setfill('0') << setw(2) << coil; // _coilCC
 406             variableSuffixStream << framesFileNameSuffix << setfill('0') << setw(2) << frame ; // _frameFF
 407             fileNameSuffixes.push_back(variableSuffixStream.str());
 408             variableSuffixStream.str(""); // emtpy string, we don't want to accumulate strings among iterations
 409         }
 410     }
 411     string dataFileNamePrefixWithDims = OpenCLIPER::Data::buildFileNamePrefix(dataFileNamePrefix, pArraysDims->at(0));
 412     (reinterpret_cast<Data*>(this))->loadRawHostData(dataFileNamePrefixWithDims, pArraysDims, pDynDims, fileNameSuffixes,
 413             fileNameExtension);
 414 #ifdef KData_DEBUG
 415     CERR("KData size: " << getData()->size() << std::endl);
 416 #endif
 417 }
 418
 419 /**
 420  * @brief Save data of NDArray objects to a group of files (every NDArray contains one image and is stored into a file).
 421  *
 422  * Files of this group have names with the format
 423  * \<fileNamePrefix\>\<dims\>\<coilsFileNameSuffix\>\<coilNumber\>\<framesFileNameSuffix\>\<frameNumber\>\<fileNameExtension\>,
 424  * where \<dims\> is a string with the format \<width\>x\<height\>x\<depth\>, \<coilNumber\> is the identifier of the coil
 425  * used for image adquisition and \<frameNumber\> is the identifier of the time of image adquisition.
 426  * @param[in] syncSource set the data source used for saving (IMAGE_ONLY or BUFFER_ONLY)
 427  * @param[in] fileNamePrefix fixed part of the file name before the variable part
 428  * @param[in] coilsFileNameSuffix part of the file name before the coil number
 429  * @param[in] framesFileNameSuffix part of the file name before the frame number
 430  * @param[in] fileNameExtension extension for file names
 431  */
 432 void KData::saveRawHostData(const SyncSource syncSource, const std::string fileNamePrefix,
 433                             const std::string coilsFileNameSuffix, const std::string framesFileNameSuffix,
 434                             const std::string fileNameExtension) {
 435     dimIndexType numFrames = getDynDimsTotalSize();
 436     numCoilsType numCoils = pData->size() / numFrames;
 437     stringstream variableSuffixStream;
 438     std::vector<std::string> fileNameSuffixes;
 439     for (dimIndexType frame = 0; frame < numFrames; frame++) {
 440         for (numCoilsType coil = 0; coil < numCoils; coil++) {
 441             variableSuffixStream << coilsFileNameSuffix << setfill('0') << setw(2) << coil; // _coilCC
 442             variableSuffixStream << framesFileNameSuffix << setfill('0') << setw(2) << frame ; // _frameFF
 443             fileNameSuffixes.push_back(variableSuffixStream.str());
 444             variableSuffixStream.str(""); // emtpy string, we don't want to accumulate strings among iterations
 445         }
 446     }
 447     (reinterpret_cast<Data*>(this))->saveRawHostData(syncSource, fileNamePrefix, fileNameSuffixes, fileNameExtension);
 448 #ifdef KData_DEBUG
 449     CERR("KData size: " << getData()->size() << std::endl);
 450 #endif
 451 }
 452
 453 /*
 454 void KData::save(const std::string fileName, const vector<string> variableNames, const SyncSource syncSource) {
 455 }
 456 */
 457
 458 void KData::setApp(std::shared_ptr<CLapp> pCLapp, SyncSource hostDeviceSync) {
 459     Data* pData = this;
 460     pData->Data::setApp(pCLapp, hostDeviceSync);
 461     if (pSensitivityMapsData != nullptr) {
 462         sensitivityMapsDataHandle = pCLapp->addData(pSensitivityMapsData, hostDeviceSync);
 463     }
 464     if (pSensitivityMapsRMS != nullptr) {
 465         pSensitivityMapsRMS->setApp(pCLapp, hostDeviceSync);
 466     }
 467     if (pSamplingMasksData != nullptr) {
 468         samplingMasksDataHandle = pCLapp->addData(pSamplingMasksData, SyncSource::BUFFER_ONLY);
 469     }
 470 }
 471
 472 /**
 473  * @brief Stores image/volume spatial and temporal dimensions in dataDims field.
 474  */
 475 void KData::calcDataDims() {
 476     Data::calcDataDims();
 477     pDataDimsVector->at(NumCoilsPos) = nCoils;
 478 }
 479
 480 /**
 481  * @brief Creates a new KData object with data following a pattern (odd negative numbers for real part, even negative numbers for imaginary part).
 482  *
 483  * @param[in] width width of the data (number of columns)
 484  * @param[in] height height of the data (number of rows)
 485  * @param[in] numCoils number of coils
 486  * @param[in] numFrames number of frames
 487  * @return a pointer to the new KData object
 488  */
 489 KData* KData::genTestKData(dimIndexType width, dimIndexType height, dimIndexType numFrames, numCoilsType numCoils) {
 490     CERR("width: " << width << std::endl << "height: " << height << std::endl);
 491     vector<dimIndexType> *pDimsInputImage, *pDimsSensitivityMap, *pDimsOutputImage, *pDimsSamplingMasks;
 492     vector <complexType>* pImageData;
 493     vector <complexType>* pSensitivityMap;
 494     vector <dimIndexType>* pSamplingMask;
 495     complexType elementImage, elementSensitivityMap;
 496     dimIndexType elementSamplingMask;
 497     std::vector<NDArray*>* pObjNDArraysImage = new std::vector<NDArray*>();
 498     std::vector<NDArray*>* pObjNDArraysSensitivityMap = new std::vector<NDArray*>();
 499     std::vector<NDArray*>* pObjNDArraysSamplingMasks = new std::vector<NDArray*>();
 500     std::vector<dimIndexType>* pInputDynDims = new std::vector<dimIndexType>({numFrames});
 501     std::vector<dimIndexType>* pSamplingMasksDynDims = new std::vector<dimIndexType>({numFrames});
 502     CERR("Creating sensitivity maps...\n");
 503     realType realElement = 1.0;
 504     for (numCoilsType coilId = 0; coilId < numCoils; coilId++) {
 505         pSensitivityMap = new vector<complexType>();
 506         for (index1DType i = 0; i < (width * height) * 2; i+=2) {
 507             elementSensitivityMap.real((-1.0)*(realElement)); // (-1, -3, -5, -7, ...)
 508             elementSensitivityMap.imag((-1.0)*(realElement+1)); // (-2, -4, -6, -8, ...)
 509             pSensitivityMap->push_back(elementSensitivityMap);
 510             realElement += 2.0;
 511         }
 512         pDimsSensitivityMap = new vector<dimIndexType>({width, height});
 513         NDArray* pObjNDArraySensitivityMap = NDArray::createNDArray<complexType>(pDimsSensitivityMap, pSensitivityMap);
 514         //new ConcreteNDArray<complexType>(pDimsSensitivityMap, pSensitivityMap);
 515         CERR("Sensitivity map coilId: " << coilId << std::endl);
 516         CERR(pObjNDArraySensitivityMap->hostDataToString("Sensitivity map") << std::endl);
 517         pObjNDArraysSensitivityMap->push_back(pObjNDArraySensitivityMap);
 518     }
 519     SensitivityMapsData* pSensitivityMapsData = new SensitivityMapsData(pObjNDArraysSensitivityMap, numCoils);
 520     CERR ("Creating sensitivity maps done.\n");
 521     CERR ("Creating KData...\n");
 522     realElement = 1.0;
 523     for (dimIndexType dynId = 0; dynId < pInputDynDims->at(0); dynId++) {
 524         for (numCoilsType coilId = 0; coilId < numCoils; coilId++) {
 525             pImageData = new vector<complexType>();
 526             for (index1DType i = 0; i < (width * height) * 2; i+=2) {
 527                 elementImage.real(realElement); // (1, 3, 5, 7, ...)
 528                 elementImage.imag(realElement+1); // (2, 4, 6, 8, ...)
 529                 pImageData->push_back(elementImage);
 530                 realElement += 2.0;
 531             }
 532             pDimsInputImage = new vector<dimIndexType>({width, height});
 533             pDimsOutputImage = new vector<dimIndexType>(*pDimsInputImage);
 534
 535             NDArray* pObjNDArrayImage = NDArray::createNDArray<complexType>(pDimsInputImage, pImageData);
 536             CERR("Initial data" << std::endl);
 537             CERR("coilId: " << coilId << "\tdynId: " << dynId << std::endl);
 538             CERR(pObjNDArrayImage->hostDataToString("NDArray:") << std::endl);
 539
 540             CERR("Objeto creado, ndims: " << (numberOfDimensionsType) pObjNDArrayImage->getNDims() << std::endl);
 541             // pObjNDArrays only needed for Data constructor but now is abstract => can't be instantiated
 542             // vector<NDArray*>* pObjNDArrays = new vector<NDArray*>;
 543             // pObjNDArrays->push_back(pObjNDArray);
 544             CERR("Creando objeto de la clase KData" << std::endl);
 545
 546             vector<realType>* pPixelSize = new vector<realType>({1});
 547             pObjNDArraysImage->push_back(pObjNDArrayImage);
 548         }
 549     }
 550     CERR ("Creating KData done.\n");
 551     CERR ("Creating sampling masks...\n");
 552     elementSamplingMask = 1;
 553     for (dimIndexType dynId = 0; dynId < pInputDynDims->at(0); dynId++) {
 554         pSamplingMask = new vector<dimIndexType>();
 555         for (index1DType i = 0; i < height; i++) {
 556             pSamplingMask->push_back(elementSamplingMask);
 557             elementSamplingMask = (elementSamplingMask + 1) % 2;
 558         }
 559         pDimsSamplingMasks = new vector<dimIndexType>({height});
 560
 561         NDArray* pObjNDArraySamplingMasks = NDArray::createNDArray<dimIndexType>(pDimsSamplingMasks, pSamplingMask);
 562         CERR("Initial data" << std::endl);
 563         CERR("dynId: " << dynId << std::endl);
 564         CERR(pObjNDArraySamplingMasks->hostDataToString("NDArray:") << std::endl);
 565
 566         CERR("Objeto creado, ndims: " << (numberOfDimensionsType) pObjNDArraySamplingMasks->getNDims() << std::endl);
 567         // pObjNDArrays only needed for Data constructor but now is abstract => can't be instantiated
 568         // vector<NDArray*>* pObjNDArrays = new vector<NDArray*>;
 569         // pObjNDArrays->push_back(pObjNDArray);
 570         pObjNDArraysSamplingMasks->push_back(pObjNDArraySamplingMasks);
 571     }
 572     SamplingMasksData* pSamplingMasksData = new SamplingMasksData(pObjNDArraysSamplingMasks, pSamplingMasksDynDims);
 573     CERR ("Creating sampling masks done.\n");
 574     pSamplingMasksData->externalToInternalFormat();
 575     for (unsigned int i = 0; i < pSamplingMasksData->getDynDimsTotalSize(); i ++) {
 576         CERR(pSamplingMasksData->getData()->at(i)->hostDataToString(string("Sampling mask (") + to_string(i) + string(")")));
 577     }
 578
 579     Data* pKData = new KData();
 580     pKData->setData(pObjNDArraysImage);
 581     pKData->setDynDims(pInputDynDims);
 582     (static_cast<KData*>(pKData))->setSensitivityMapsData(pSensitivityMapsData);
 583     (static_cast<KData*>(pKData))->setSamplingMasksData(pSamplingMasksData);
 584     (static_cast<KData*>(pKData))->setNCoils(numCoils);
 585     return (static_cast<KData*>(pKData));
 586 }
 587
 588 /**
 589  * @brief Saves data to a file in matlab format
 590  *
 591  * @param[in] fileName name of the file tha data will be saved to
 592  * @param[in] varNames names of the matlab variables stored in the file (for K-space data, sensitivity maps and
 593  * sampling masks, in this order)
 594  * @param[in] syncSource source of data: (OpenCL buffer or image)
 595  */
 596 void KData::matlabSave(string fileName, vector<string> varNames, SyncSource syncSource) {
 597     mat_t* matfp;
 598     vector<dimIndexType> matVarKDataDimsVector, matVarSensMapsDimsVector, matVarSampMasksDimsVector;
 599     matfp = Mat_CreateVer(fileName.c_str(),NULL,MAT_FT_DEFAULT);
 600     if ( NULL == matfp ) {
 601         throw invalid_argument(string("Error creating MAT file")  + fileName);
 602     }
 603     if (varNames.size() < 3) {
 604         throw invalid_argument("variable names for Data in k-space, sensitivity maps and sampling masks should be provided");
 605     }
 606     // Only image sequences with the same spatial dimensions are supported
 607     dimIndexType numMatVarKDataElements = pData->at(0)->size() * getNCoils() * getDynDimsTotalSize();
 608     if (NDARRAYWIDTH(pData->at(0)) == 0) {
 609         throw invalid_argument ("Invalid data size (width is 0)");
 610     }
 611
 612     MatVarInfo* pMatVarInfo = NDArray::newMatVarInfo(elementDataType, numMatVarKDataElements);
 613     vector<dimIndexType>* joinedDimensions = new vector<dimIndexType>(*(pData->at(0)->getDims()));
 614     joinedDimensions->push_back(getNCoils());
 615     joinedDimensions->insert(std::end(*joinedDimensions), std::begin(*getDynDims()), std::end(*getDynDims()));
 616     // Update matlab dimensions and rank with spatial dimensions (true for swapping width and depth dimensions)
 617     // Update matlab dimensions and rank with coils dimensions
 618     // Update matlab dimensions and rank with temporal dimensions
 619     pMatVarInfo->updateDimsAndRank(joinedDimensions,true);
 620     fillMatlabVarInfo(matfp, varNames.at(0), pMatVarInfo, syncSource);
 621     delete(pMatVarInfo);
 622     joinedDimensions->clear();
 623
 624     if (pSensitivityMapsData != nullptr) {
 625         dimIndexType numMatVarSensMapsElements = pSensitivityMapsData->getData()->at(0)->size() * getNCoils();
 626         pMatVarInfo = NDArray::newMatVarInfo(pSensitivityMapsData->getElementDataType(), numMatVarSensMapsElements);
 627         // Update matlab dimensions and rank with spatial dimensions (true for swapping width and depth dimensions)
 628         // Update matlab dimensions and rank with coils dimensions
 629         joinedDimensions->insert(std::end(*joinedDimensions), std::begin(*(pSensitivityMapsData->pData->at(0)->getDims())),
 630                                  std::end(*(pSensitivityMapsData->pData->at(0)->getDims())));
 631         joinedDimensions->push_back(pSensitivityMapsData->getNCoils());
 632         pMatVarInfo->updateDimsAndRank(joinedDimensions,true);
 633         pSensitivityMapsData->fillMatlabVarInfo(matfp, varNames.at(1), pMatVarInfo, syncSource);
 634         delete(pMatVarInfo);
 635         joinedDimensions->clear();
 636     } else {
 637         CERR("KData::matlabSave: KData without sensitivity maps\n");
 638     }
 639     if (pSamplingMasksData != nullptr) {
 640         pSamplingMasksData->internalToExternalFormat();
 641         dimIndexType numMatVarSampMasksElements = pSamplingMasksData->getData()->at(0)->size() * pSamplingMasksData->getDynDimsTotalSize();
 642         pMatVarInfo = NDArray::newMatVarInfo(pSamplingMasksData->getElementDataType(), numMatVarSampMasksElements);
 643         // First dimension of sampling masks is the number of image lines (its height), second and following
 644         // are the temporal dimensions;
 645         vector<dimIndexType>* pSamplingMasksDimensions = new vector<dimIndexType>();
 646         // Add spatial dimension: number of columns (width)
 647         pSamplingMasksDimensions->insert(pSamplingMasksDimensions->end(),
 648                                          pSamplingMasksData->pData->at(0)->getDims()->begin(),
 649                                          pSamplingMasksData->pData->at(0)->getDims()->end());
 650         /* Add number of rows (is 1): matlab needs at least two dimensions (and a sampling mask must be a row vector
 651          * with a number of columns equal to the image number of lines)
 652         */
 653         pSamplingMasksDimensions->push_back(1);
 654         // Add temporal dimensions
 655         pSamplingMasksDimensions->insert(pSamplingMasksDimensions->end(),
 656                                          pSamplingMasksData->getDynDims()->begin(),
 657                                          pSamplingMasksData->getDynDims()->end());
 658         // Update matlab dimensions and rank with sampling mask dimensions dimensions
 659         pMatVarInfo->updateDimsAndRank(pSamplingMasksDimensions, true);
 660         pSamplingMasksData->fillMatlabVarInfo(matfp, varNames.at(2), pMatVarInfo, syncSource);
 661         delete(pMatVarInfo);
 662         // Format must be reverted to internal after saving (sampling masks may be reused without beeing recreated)
 663         pSamplingMasksData->externalToInternalFormat();
 664     } else {
 665         CERR("KData::matlabSave: KData without sampling masks\n");
 666     }
 667
 668     Mat_Close(matfp);
 669 }
 670
 671
 672 } /* namespace OpenCLIPER */
 673