[Insight-users] Dicom Loader Error
Abayiz
abayiz at yahoo.com
Mon Jun 13 04:23:15 EDT 2011
Hello ITK users,
I am using the following code to load dicom series, but the problem is, it is not loading the files in the original order. The loaded files are in wrong sequence. Which part could be incorrect?
This is really very important to me, and sincerely hope to receive your help soon.
Many thanks in advance!
The code is as follows:
*********************************
BaseData* DICOMLoader::LoadData(ContentItem* contentItem)
{
Uint32 min = 4294967295;
Uint32 max = 0;
vector<string> files;
if(contentItem == NULL)
{
return NULL;
}
if(contentItem->HasChildItem())
{
ContentItem* item = contentItem->GetFirstItem();
while(item != NULL)
{
files.push_back(item->GetFSPath());
item = contentItem->GetNextItem();
}
}
else
{
files.push_back(contentItem->GetFSPath());
}
sort(files.begin(), files.end());
DcmFileFormat *dfile = new DcmFileFormat();
long int val_long;
int width, height, depth = 1;
OFString str_modality;
ModalityType modality;
BaseData *result = NULL;
Uint16 *dicom_buf;
unsigned long nObject;
DicomImage *di = NULL;
double m_min =4294967295,m_max = 0;
map<Sint32, double> sliceDistances;
EP_Interpretation epi;
OFString info;
E_TransferSyntax xfer;
vtkImageData *ImageData = NULL;
vtkLookupTable *lookup = NULL;
for(int i= 0; i < files.size(); i++)
{
if(dfile->loadFile((files[i]).c_str()).good())
{
DcmDataset *ds = dfile->getDataset();
xfer = ds->getOriginalXfer();
if(xfer >= EXS_JPEGProcess1TransferSyntax && EXS_JPEGProcess14SV1TransferSyntax >= xfer )
{
DJDecoderRegistration::registerCodecs();
}
di = new DicomImage(dfile, xfer, 0UL, 0UL, 0UL);
if(!di)
{
continue;
}
if(!i)
{
if(ds->findAndGetOFString(DCM_Modality,info).bad())
{
cout<<"Unable to get element: DCM_Modality\n";
}
else
{
if ((strcasecmp(info.c_str(), "CT") == 0) ||
(strcasecmp(info.c_str(), "DF") == 0) ||
(strcasecmp(info.c_str(), "DS") == 0) ||
(strcasecmp(info.c_str(), "DX") == 0) ||
(strcasecmp(info.c_str(), "MG") == 0) ||
(strcasecmp(info.c_str(), "PX") == 0) ||
(strcasecmp(info.c_str(), "RF") == 0) ||
(strcasecmp(info.c_str(), "RG") == 0) ||
(strcasecmp(info.c_str(), "XA") == 0))
{
modality = M_CT;
result = new CTData();
}
else if ((strcasecmp(info.c_str(), "MA") == 0) ||
(strcasecmp(info.c_str(), "MR") == 0) ||
(strcasecmp(info.c_str(), "MS") == 0))
{
modality = M_MRI;
result = new MRData();
}
else if ((strcasecmp(info.c_str(), "ST") == 0) ||
(strcasecmp(info.c_str(), "NM") == 0))
{
modality = M_SPECT;
result = new NMData();
}
else if ((strcasecmp(info.c_str(), "PT") == 0))
{
modality = M_PET;
result = new CTData();
}
else if ((strcasecmp(info.c_str(), "CR") == 0))
{
modality = M_CR;
result = new CRData();
}
else if ((strcasecmp(info.c_str(), "US") == 0))
{
modality = M_US;
result = new USData();
}
else
{
modality = UNKNOWN;
result = new BaseData();
}
}
info.clear();
if(ds->findAndGetOFString(DCM_PatientsName,info).bad())
{
cout<<"Unable to get element: DCM_PatientsName\n";
}
else
{
result->SetPatientName(info.c_str());
cout<<"Patient Ifo: "<<info<<endl;
}
info.clear();
if(ds->findAndGetOFString(DCM_PatientID,info).bad())
{
cout<<"Unable to get element: DCM_PatientID\n";
}
else
{
result->SetPatientID(info.c_str());
cout<<"Patient Ifo: "<<info<<endl;
}
info.clear();
if(ds->findAndGetOFString(DCM_PatientsSex,info).bad())
{
cout<<"Unable to get element: DCM_PatientsSex\n";
}
else
{
result->SetPatientSex(info.c_str());
cout<<"Patient Ifo: "<<info<<endl;
}
info.clear();
epi = di->getPhotometricInterpretation();
if(ds->findAndGetLongInt(DCM_Columns, val_long).bad())
{
cout<<"Unable to get element: DCM_Columns...\n";
val_long = 0;
}
width = val_long;
if(ds->findAndGetLongInt(DCM_Rows, val_long).bad())
{
cout<<"Unable to get element: DCM_Rows...\n";
val_long = 0;
}
height = val_long;
depth = di->getFrameCount() == 1 ? files.size() : di->getFrameCount();
if(di->getFrameCount() > 1)
{
cout<<"This is a multiframe instance ...\n";
LoadMultiFrame(ds, di, result, sliceDistances);
}
else
{
cout<<"This is a single frame instance ...\n";
LoadSingleFrame(ds, di, result, sliceDistances, i);
}
const DiPixel *dmp = di->getInterData();
switch(dmp->getRepresentation())
{
case EPR_Uint8:
result->GetScalar()->SetDataType(VTK_UNSIGNED_CHAR);
break;
case EPR_Sint8:
result->GetScalar()->SetDataType(VTK_CHAR);
break;
case EPR_Uint16:
result->GetScalar()->SetDataType(VTK_UNSIGNED_SHORT);
break;
case EPR_Sint16:
result->GetScalar()->SetDataType(VTK_SHORT);
break;
case EPR_Uint32:
result->GetScalar()->SetDataType(VTK_UNSIGNED_INT);
break;
case EPR_Sint32:
result->GetScalar()->SetDataType(VTK_INT);
break;
default:
break;
}
result->GetScalar()->SetInitials ( width, height, depth, result->GetScalar()->GetSpacingPtr()[0],
result->GetScalar()->GetSpacingPtr()[1], result->GetScalar()->GetSpacingPtr()[2], 0, 0,0 );
result->GetScalar()->CreateData();
ImageData = result->GetScalar()->GetImageData();
}
if(di->getFrameCount() > 1)
{
cout<<"This is a multiframe instance ...\n";
LoadMultiFrame(ds, di, result, sliceDistances);
}
else
{
cout<<"This is a single frame instance ...\n";
LoadSingleFrame(ds, di, result, sliceDistances, i);
}
const DiPixel *dmp = di->getInterData();
const void *pixelData = dmp->getData();
unsigned int* out;
if(epi == EPI_RGB )
{
cout<<"getOutputDataSize(): "<<di->getOutputDataSize(17)<<endl;
out = (unsigned int *)di->getOutputData(17,0,0);
}
else
{
cout<<"getOutputDataSize(): "<<di->getOutputDataSize(di->getDepth())<<endl;
}
cout<<"REP: "<<dmp->getRepresentation()<<endl;
switch(dmp->getRepresentation())
{
case EPR_Uint8:
{
if(epi == EPI_RGB )
{
unsigned int * vtkData = (unsigned int*)ImageData->GetScalarPointer();
unsigned int k,j;
unsigned int a[3];
unsigned j_max = height * width;
unsigned k_max = 3 * j_max;
for(k= 0, j = 0; k<k_max && j < j_max;j++)
{
a[0] = out[k++] ;
if(a[0] < min)
min = a[0];
if(a[0] > max)
max = a[0];
a[1] = out[k++] ;
if(a[1] < min)
min = a[1];
if(a[1] > max)
max = a[1];
a[2] = out[k++] ;
if(a[2] < min)
min = a[2];
if(a[2] > max)
max = a[2];
lookup->SetTableValue(j, ((double)a[0])/max/*131070.0*/, ((double)a[1])/max/*131070.0*/, ((double)a[2])/max/*131070.0*/, 0.5);
vtkData[0] = j;
vtkData = vtkData ++;
}
ImageData->GetPointData()->GetScalars()->SetLookupTable(lookup);
cout<<"RGB_min: "<<min<<" RGB_max:"<<max<<endl;
}
else
{
for (int j = 0; j < height; j++){
memcpy((Uint8*)ImageData->GetScalarPointer() + i * height * width + j*width,
(Uint8*)pixelData+ (height-(j+1))*width,
/*height * */width * di->getFrameCount() * sizeof(Uint8));
}
}
}
break;
case EPR_Sint8:
{
for (int j = 0; j < height; j++){
memcpy((Sint8*)ImageData->GetScalarPointer() + i * height * width + j*width,
(Sint8*)pixelData+ (height-(j+1))*width,
/*height * */width * di->getFrameCount() * sizeof(Sint8));
}
}
break;
case EPR_Uint16:
{
for (int j = 0; j < height; j++){
memcpy((Uint16*)ImageData->GetScalarPointer() + i * height * width + j*width,
(Uint16*)pixelData+ (height-(j+1))*width,
/*height * */width * di->getFrameCount() * sizeof(Uint16));
}
}
break;
case EPR_Sint16:
{
for (int j = 0; j < height; j++){
memcpy((Sint16*)ImageData->GetScalarPointer() + i * height * width + j*width,
(Sint16*)pixelData+ (height-(j+1))*width,
/*height * */width * di->getFrameCount() * sizeof(Sint16));
}
}
break;
case EPR_Uint32:
{
for (int j = 0; j < height; j++){
memcpy((Uint32*)ImageData->GetScalarPointer() + i * height * width + j*width,
(Uint32*)pixelData+ (height-(j+1))*width,
/*height * */width * di->getFrameCount() * sizeof(Uint32));
}
}
break;
case EPR_Sint32:
{
for (int j = 0; j < height; j++){
memcpy((Sint32*)ImageData->GetScalarPointer() + i * height * width + j*width,
(Sint32*)pixelData+ (height-(j+1))*width,
/*height * */width * di->getFrameCount() * sizeof(Sint32));
}
}
break;
default:
break;
}
result->SetDataSet(i, ds);
ImageData->Modified();
delete di;
dfile->clear();
if(xfer >= EXS_JPEGProcess1TransferSyntax && EXS_JPEGProcess14SV1TransferSyntax >= xfer )
{
DJDecoderRegistration::cleanup();
}
}
else
{
cout<<"can not load file: "<<files[i]<<endl;
}
}
map<Sint32, double>::iterator iter = sliceDistances.begin();
double sliceDistances1, sliceDistances2;
if (sliceDistances.size() == 1) {
sliceDistances1 = sliceDistances2 = iter->second;
}
else
{
iter++;
sliceDistances1 = iter->second;
iter++;
sliceDistances2 = iter->second;
}
printf("Real spacing[2] is %lf ....\n\n\n", sliceDistances1 -sliceDistances2);
if(sliceDistances1 != sliceDistances2)
{
if(sliceDistances1 > sliceDistances2)
{
result->GetScalar()->SetSpacing(result->GetScalar()->GetSpacingPtr()[0], result->GetScalar()->GetSpacingPtr()[1], sliceDistances1 -sliceDistances2);
ImageData->SetSpacing(result->GetScalar()->GetSpacingPtr()[0], result->GetScalar()->GetSpacingPtr()[1], sliceDistances1 -sliceDistances2);
}
else
{
result->GetScalar()->SetSpacing(result->GetScalar()->GetSpacingPtr()[0], result->GetScalar()->GetSpacingPtr()[1], sliceDistances2 -sliceDistances1);
ImageData->SetSpacing(result->GetScalar()->GetSpacingPtr()[0], result->GetScalar()->GetSpacingPtr()[1], sliceDistances2 -sliceDistances1);
}
}
ImageData->GetPipelineInformation()->Clear();
ImageData->UpdateInformation();
ImageData->UpdateData();
ImageData->Update();
int dimensions[3];
double spacing[3];
double origin[3];
ImageData->GetDimensions(dimensions);
ImageData->GetSpacing(spacing);
ImageData->GetOrigin(origin);
unifyOrientation(result);
ImageData = result->GetScalar()->GetImageData();
ImageData->SetDimensions(dimensions);
ImageData->SetSpacing(spacing);
ImageData->SetOrigin(origin);
ImageData->Update();
result->GetScalar()->GetImageData()->GetDimensions();
return result;
}
**********************
With regards,
Abayiz
More information about the Insight-users
mailing list