-
Notifications
You must be signed in to change notification settings - Fork 1
ImportingExporting
Adrian Quintana edited this page Dec 11, 2017
·
1 revision
Spider and Xmipp file formats are the same for volumes, images and document files. In this document you will find some code in C, Matlab and instructions for loading Spider images with imagej.
You can use the following code to read and write Spider images:
/* Spider read/write routines ---------------------------------------------- */
#include <stdlib.h>
#include <cstdio>
#include <math.h>
#include <sys/stat.h>
typedef struct
{
/* 1 */
float fNslice; // NUMBER OF SLICES (PLANES) IN VOLUME
// (=1 FOR AN IMAGE) FOR NEW LONG LABEL
// FORMAT THE VALUE OF NSLICE STORED IN
// THE FILE IS NEGATIVE.
/* 2 */
float fNrow; // NUMBER OF ROWS PER SLICE (Y)
/* 3 */
float fNrec; // TOTAL NUMBER OF RECORDS (SEE NOTE #3).
/* 4 */
float fNlabel; // AUXILIARY NUMBER TO COMPUTE TOTAL NUMBER OF RECS
/* 5 */
float fIform; // FILE TYPE SPECIFIER.
// +3 FOR A 3-D FILE (FLOAT)
// +1 FOR A 2-D IMAGE (FLOAT)
// -1 FOR A 2-D FOURIER TRANSFORM
// -3 FOR A 3-D FOURIER TRANSFORM
// -5 FOR A NEW 2-D FOURIER TRANSFORM
// -7 FOR A NEW 3-D FOURIER TRANSFORM
// +8 FOR A 2-D EIGHT BIT IMAGE FILE
// +9 FOR A 2-D INT IMAGE FILE
// 10 FOR A 3-D INT IMAGE FILE
// 11 FOR A 2-D EIGHT BIT COLOR IMAGE FILE
/* 6 */
float fImami; // MAXIMUM/MINIMUM FLAG. IS SET AT 0 WHEN THE
// FILE IS CREATED, AND AT 1 WHEN THE MAXIMUM AND
// MINIMUM HAVE BEEN COMPUTED, AND HAVE BEEN STORED
// INTO THIS LABEL RECORD (SEE FOLLOWING WORDS)
/* 7 */
float fFmax; // MAXIMUM VALUE
/* 8 */
float fFmin; // MINIMUM VALUE
/* 9 */
float fAv; // AVERAGE VALUE
/* 10*/
float fSig; // STANDARD DEVIATION. A VALUE OF -1. INDICATES
// THAT SIG HAS NOT BEEN COMPUTED PREVIOUSLY.
/* 11*/
float fIhist; // FLAG INDICATING IF THE HISTOGRAM HAS BE
// COMPUTED. NOT USED IN 3D FILES!
/* 12*/
float fNcol; // NUMBER OF PIXELS PER LINE (Columns X)
/* 13*/
float fLabrec; // NUMBER OF LABEL RECORDS IN FILE HEADER
/* 14*/
float fIangle; // FLAG THAT TILT ANGLES HAVE BEEN FILLED
/* 15*/
float fPhi; // EULER: ROTATIONAL ANGLE
/* 16*/
float fTheta; // EULER: TILT ANGLE
/* 17*/
float fPsi; // EULER: PSI = TILT ANGLE
/* 18*/
float fXoff; // X TRANSLATION
/* 19*/
float fYoff; // Y TRANSLATION
/* 20*/
float fZoff; // Z TRANSLATION
/* 21*/
float fScale; // SCALE
/* 22*/
float fLabbyt; // TOTAL NUMBER OF BYTES IN LABEL
/* 23*/
float fLenbyt; // RECORD LENGTH IN BYTES
char fNada[24]; // this is a spider incongruence
/* 30*/
float fFlag; // THAT ANGLES ARE SET. 1 = ONE ADDITIONAL
// ROTATION IS PRESENT, 2 = ADDITIONAL ROTATION
// THAT PRECEEDS THE ROTATION THAT WAS STORED IN
// 15 FOR DETAILS SEE MANUAL CHAPTER VOCEUL.MAN
/* 31*/
float fPhi1;
/* 32*/
float fTheta1;
/* 33*/
float fPsi1;
/* 34*/
float fPhi2;
/* 35*/
float fTheta2;
/* 36*/
float fPsi2;
double fGeo_matrix[3][3]; // x9 = 72 bytes: Geometric info
float fAngle1; // angle info
float fr1;
float fr2; // lift up cosine mask parameters
/** Fraga 23/05/97 For Radon transforms **/
float RTflag; // 1=RT, 2=FFT(RT)
float Astart;
float Aend;
float Ainc;
float Rsigma; // 4*7 = 28 bytes
float Tstart;
float Tend;
float Tinc; // 4*3 = 12, 12+28 = 40B
/** Sjors Scheres 17/12/04 **/
float Weight; // For Maximum-Likelihood refinement
float Flip; // 0=no flipping operation (false), 1=flipping (true)
char fNada2[576]; // empty 700-76-40=624-40-8= 576 bytes
/*212-214*/
char szIDat[12]; // LOGICAL * 1 ARRAY DIMENSIONED 10, CONTAINING
// THE DATE OF CREATION (10 CHARS)
/*215-216*/
char szITim[8]; // LOGICAL * 1 ARRAY DIMENSIONED 8, CONTAINING
// THE TIME OF CREATION (8 CHARS)
/*217-256*/
char szITit[160]; // LOGICAL * 1 ARRAY DIMENSIONED 160
}
SpiderHeader;
size_t FREAD(void *dest, size_t size, size_t nitems, FILE *&fp, int reverse)
{
size_t retval;
if (!reverse)
retval = fread(dest, size, nitems, fp);
else
{
char *ptr = (char *)dest;
int end = 0;
retval = 0;
for (int n = 0; n < nitems; n++)
{
for (int i = size - 1; i >= 0; i--)
{
if (fread(ptr + i, 1, 1, fp) != 1)
{
end = 1;
break;
}
}
if (end)
break;
else
retval++;
ptr += size;
}
}
return retval;
}
/* This function reads a Spider volume. Call it as
int dim[3];
float *data=NULL;
readSpiderFile("myfile.vol",'V',dim,&data);
readSpiderFile("myfile.img",'I',dim,&data);
Code errors:
0 - OK
1 - Cannot open file
2 - File is not a Spider volume
3 - Problem when computing the size of the file
4 - filetype is not 'V' or 'I'
*/
int readSpiderFile(const char *filename, char filetype,
int dim[], float **data)
{
FILE* fp=NULL;
union {
unsigned char c[4];
float f;
} file_type;
int fileReversed=0, machineReversed=0, reversed=0;
SpiderHeader header;
unsigned long usfNcol, usfNrow, usfNslice, usfHeader, size, tmpSize,
volSize, n, currentPosition;
struct stat info;
float *ptr;
/* Set dimensions to 0, just in case it cannot be read */
dim[0]=0;
dim[1]=0;
dim[2]=0;
/* Check that the filetype is correct */
if (filetype!='V' and filetype!='I')
return 4;
/* Open file */
if ((fp = fopen(filename, "rb")) == NULL)
return 1;
/* Check that the input file is really a Spider volume */
currentPosition=ftell(fp);
/* Check file type */
fseek(fp, currentPosition+16, SEEK_SET);
for (int i = 0; i < 4; i++)
fread(&(file_type.c[i]), sizeof(unsigned char), 1, fp);
fseek(fp, currentPosition+0, SEEK_SET);
switch (filetype)
{
case 'V':
if (file_type.c[0]== 0 && file_type.c[1]== 0 &&
file_type.c[2]== 64 && file_type.c[3]==64)
{
fileReversed=0;
}
else if (file_type.c[0]==64 && file_type.c[1]==64 &&
file_type.c[2]== 0 && file_type.c[3]== 0)
{
fileReversed=1;
}
else
{
fclose(fp);
return 2;
}
break;
case 'I':
if (file_type.c[0]== 0 && file_type.c[1]== 0 &&
file_type.c[2]==128 && file_type.c[3]==63)
{
fileReversed=0;
}
else if (file_type.c[0]==63 && file_type.c[1]==128 &&
file_type.c[2]== 0 && file_type.c[3]== 0)
{
fileReversed=1;
}
else
{
fclose(fp);
return 2;
}
break;
}
/* Check machine type */
file_type.f = 1;
if (file_type.c[0]==63 && file_type.c[1]==128 && file_type.c[2]==0 &&
file_type.c[3]==0)
machineReversed=1;
/* Read header */
reversed=fileReversed ^ machineReversed;
if (!reversed)
fread(&header, sizeof(SpiderHeader), 1, fp);
else
{
FREAD(&header, sizeof(float), 36, fp, true);
FREAD(&header.fGeo_matrix, sizeof(double), 9, fp, true);
FREAD(&header.fAngle1, sizeof(float), 13, fp, true);
FREAD(&header.fNada2, sizeof(char), 756, fp, true);
}
/* Compute file size, header size and volume dimensions */
usfNcol = (unsigned long) header.fNcol;
usfNrow = (unsigned long) header.fNrow;
usfNslice = (unsigned long) header.fNslice;
usfHeader = (unsigned long) header.fNcol *
(unsigned long) header.fLabrec * sizeof(float);
if (fstat(fileno(fp), &info))
{
fclose(fp);
return 3;
}
/* Readjust the number of rows in "aberrant" images*/
if (filetype=='I' || header.fIform == 1)
if (usfNcol*usfNrow*sizeof(float) == info.st_size)
{
usfNrow = (unsigned long)(--header.fNrow);
--header.fNrec;
}
/* Check that the file size is correct */
switch (filetype)
{
case 'I':
size = usfHeader + usfNcol * usfNrow * sizeof(float);
if ((size != info.st_size) || (header.fIform != 1))
{
fclose(fp);
return 2;
}
break;
case 'V':
size = usfHeader + usfNslice * usfNcol * usfNrow * sizeof(float);
if ((size != info.st_size) || (header.fIform != 3))
{
fclose(fp);
return 2;
}
break;
}
/* Read the extra filling header space */
tmpSize = (unsigned long)(header.fNcol * header.fLabrec * 4);
tmpSize -= sizeof(SpiderHeader);
for (unsigned int i = 0; i < tmpSize / 4; i++)
{
float tmp;
fread(&tmp, sizeof(float), 1, fp);
}
currentPosition=ftell(fp);
/* Fill the dimensions */
dim[0]=(int)header.fNcol;
dim[1]=(int)header.fNrow;
if (filetype=='V') dim[2]=(int)header.fNslice;
else dim[2]=1;
volSize = (unsigned long) dim[0]*dim[1]*dim[2];
/* Read the whole file */
if (*data!=NULL)
free(data);
*data=(float *) malloc(volSize*sizeof(float));
ptr=*data;
for (n=0; n<volSize; n++)
FREAD(ptr++, sizeof(float), 1, fp, reversed);
/* Close file */
fclose(fp);
return 0;
}
/* This function writes a Spider volume. Call it as
writeSpiderFile("myfile.vol",'V',dim,data);
writeSpiderFile("myfile.img",'I',dim,data);
Code errors:
0 - OK
1 - Cannot open file
*/
int writeSpiderFile(const char *filename, char filetype,
int dim[], float *data)
{
FILE* fp;
SpiderHeader header;
unsigned long volSize, tmpSize;
float tmp=0.0F;
int headrec;
/* Check that the filetype is correct */
if (filetype!='V' and filetype!='I')
return 4;
/* Open file for writing */
if ((fp = fopen(filename, "wb")) == NULL)
return 1;
/* Adjust header */
header.fNcol = dim[0];
header.fNrow = dim[1];
header.fNslice=dim[2];
if ((header.fNcol != 0) && (header.fNrow != 0))
{
header.fNlabel = (float)((int)(256 / header.fNcol + 1));
header.fLabrec = (float) ceil((float)(256 / (float)header.fNcol));
headrec = (int) 1024 / ((int)header.fNcol * 4);
if ((1024 % (int)header.fNcol != 0))
{
header.fNrec = header.fNrow + 1;
headrec++;
}
else
header.fNrec = header.fNrow;
header.fLabbyt = header.fNcol * header.fLabrec * 4;
header.fLenbyt = (float) header.fNcol * 4;
}
if (filetype=='V') header.fIform = 3;
else if (filetype=='I') header.fIform = 1;
header.fScale = 1;
for (unsigned i = 0; i < 3; i++)
for (unsigned j = 0; j < 3; j++)
if (i == j)
header.fGeo_matrix[i][j] = (double)1.0;
else
header.fGeo_matrix[i][j] = (double)0.0;
header.fSig = -1;
header.fImami = header.fFmax = header.fFmin = header.fAv =
header.fIhist = header.fIangle = header.fPhi = header.fTheta =
header.fPsi = header.fXoff = header.fYoff = header.fZoff =
header.fPhi1 = header.fTheta1 = header.fPsi1 =
header.fPhi2 = header.fTheta2 = header.fPsi2 = header.fAngle1 =
header.fr1 = header.fr2 = header.RTflag = header.Astart =
header.Aend = header.Ainc = header.Rsigma = header.Tstart =
header.Tend = header.Tinc = header.Weight = header.Flip = 0.0F;
header.szIDat[0] = header.szITim[0] = header.szITit[0] =
header.fNada2[0] = header.fNada[0] = 0;
// Write the header and the necessary empty space
fwrite(&header, sizeof(SpiderHeader), 1, fp);
tmpSize = (header.fNcol * header.fLabrec * 4);
tmpSize -= sizeof(SpiderHeader);
for (unsigned int i = 0; i < tmpSize/4; i++)
fwrite(&tmp, sizeof(float), 1, fp);
// Now write the data
fwrite(data, sizeof(float), (unsigned long) dim[0]*dim[1]*dim[2], fp);
// Close file
fclose(fp);
return 0;
}
To try this code you may use the following one in C++ and comparing with the Xmipp own routines. You will find the volumes and images for trying attached to this wiki page.
/* Testing ----------------------------------------------------------------- */
#include <data/volume.h>
#include <data/image.h>
void testVolumeRead(const char * filename)
{
VolumeXmipp V(filename);
int dim[3];
float *data=NULL;
// Read the spider file
int result=readSpiderFile(filename,'V',dim,&data);
if (result==0)
{
if (data[0]!=V(0,0,0))
printf("The first value of %s does not coincide\n",filename);
if (data[XSIZE(V())*YSIZE(V())*ZSIZE(V())-1]!=
V(ZSIZE(V())-1,YSIZE(V())-1,XSIZE(V())-1))
printf("The last value of %s does not coincide\n",filename);
}
else
printf("There is a problem reading %s: %d\n",filename, result);
// Write it to make sure it is correct
result=writeSpiderFile(((std::string)filename+".io").c_str(),
'V',dim,data);
if (result!=0)
printf("There is a problem writing %s: %d\n",filename, result);
free(data);
}
void testImageRead(const char * filename)
{
ImageXmipp I(filename);
int dim[3];
float *data=NULL;
// Read the spider file
int result=readSpiderFile(filename,'I',dim,&data);
if (result==0)
{
if (data[0]!=I(0,0))
printf("The first value of %s does not coincide\n",filename);
if (data[XSIZE(I())*YSIZE(I())-1]!=
I(YSIZE(I())-1,XSIZE(I())-1))
printf("The last value of %s does not coincide\n",filename);
}
else
printf("There is a problem reading %s: %d\n",filename, result);
// Write it to make sure it is correct
result=writeSpiderFile(((std::string)filename+".io").c_str(),
'I',dim,data);
if (result!=0)
printf("There is a problem writing %s: %d\n",filename, result);
free(data);
}
int main()
{
testVolumeRead("1FFKnoise.vol");
testVolumeRead("1BRDnoise.vol");
testImageRead("g0ta000001noise.xmp");
testImageRead("g0ta000001noiseLarge.xmp");
return 0;
}
Use the following two functions to read images or volumes
function I=open_image(filename,endianness)
%
% Read a projection image in spider format.
% Return -1 on error
%
% Call as open_volume('myvolume.vol','ieee-le')
% Call as open_volume('myvolume.vol','ieee-be')
%
% Yoel Shkolnisky, November 2008
fid = fopen(filename, 'rb');
if (fid==-1)
I=-1;
return;
end
header=fread(fid, 256, 'float',0,endianness);
fNsam=header(12);
fLabrec=header(13);
header_size=fNsam*fLabrec;
% calculate the true header size
if header_size~=256
fseek(fid,0,'bof');
header=fread(fid, header_size, 'float',0,endianness);
end
ncol=header(12); %number of pixels per row
nrow=header(2); %number of rows
I = fread(fid, nrow*ncol, 'float',0,endianness);
I=reshape(I,nrow,ncol);
fclose(fid);
function V=open_volume(filename,endianness)
%
% Read a volume in spider format.
% Return -1 on error
%
% Call as open_volume('myvolume.vol','ieee-le')
% Call as open_volume('myvolume.vol','ieee-be')
%
% Carlos Oscar S. Sorzano, July 2009
fid = fopen(filename, 'rb');
if (fid==-1)
I=-1;
return;
end
header=fread(fid, 256, 'float', 0, endianness);
fNsam=header(12);
fLabrec=header(13);
header_size=fNsam*fLabrec;
% calculate the true header size
if header_size~=256
fseek(fid,0,'bof');
header=fread(fid, header_size, 'float', 0, endianness);
end
nslice=header(1); %number of slices
ncol=header(12); %number of pixels per row
nrow=header(2); %number of rows
V=fread(fid, nslice*nrow*ncol, 'float', 0, endianness);
V=reshape(V,nslice,nrow,ncol);
fclose(fid);
function vol=read_set_of_images(prefix,suffix)
%
% vol=read_spider(prefix)
% vol=read_spider(prefix,suffix)
%
% Read projections in spider format.
%
% Each projection is stored in a separate file. All filenames have the same
% prefix appended with a running number.
% The output is a stack of projections where vol(:,:,k) is the k'th
% projection.
% The default filename suffix is 'xmp'.
%
% Example:
% projs=read_spider('./dataset2/g1ta');
%
% Yoel Shkolnisky, November 2008
if nargin<2
suffix='xmp';
end;
proj_filenames=dir(strcat(prefix,'*','.',suffix));
fpath=fileparts(prefix);
K=length(proj_filenames);
if K==0
error('No images found');
end
for k=1:K;
proj=open_image(fullfile(fpath,proj_filenames(k).name));
if proj==-1
error('Problem with reading file');
end
if k==1
vol=zeros(size(proj,1),size(proj,2),K);
end
vol(:,:,k)=proj;
end
If they fail, you can use the following function but you need to know the actual header size and volume/image size
function I=myOpenSpider(filename,headerSize,dim,endianness)
fid = fopen(filename, 'rb', endianness);
fseek(fid, headerSize, 'bof');
aux=cumprod(dim);
I = fread(fid, aux(length(aux)), 'float', 0, endianness);
I=reshape(I,dim);
fclose(fid);
Call it as
I=open_image('image.xmp',1024,[64 64],'ieee-le');
or
V=open_image('volume.vol',1024,[64 64 64],'ieee-le');
You may also use the Spider read/write routines http://www.mathworks.com/matlabcentral/fileexchange/22902
For adding new file formats to ImageJ you will find a very good tutorial at http://www.mcdb.ucla.edu/research/hartenstein/acardona/imagej_programming_tutorials
In the particular case of Spider file formats you may simply copy the two .class files you will find attached to this wiki page into the plugins directory of ImageJ. The spider reader is done by the Institute Curie and can be found originally at http://u759.curie.u-psud.fr/softwaresu759.
-- Main.CoSS - 08 Aug 2008
-
/1BRDnoise.vol: Volume of the bacteriorhodopsin at a size of 64x64x64
-
/1FFKnoise.vol: Volume of the larg ribosome subunit at a size of 55x55x55
-
/g0ta000001noise.xmp: Image of the bacteriorhodopsin at a size of 64x64
-
/g0ta000001noiseLarge.xmp: The same image as before but at a size of 213x213
-
/SPIDER_Reader.class: Spider reader compiled class
-
/SPIDER_Reader.java: Spider reader java source code
-
/HandleExtraFileTypes.java: Import handler java source code
-
/Sel_Reader.class: Class for reading selfiles
-
/HandleExtraFileTypes.class: Import handler compiled class