00001 #include "1B.h"
00002 #include <stdio.h>
00003 #include <fstream>
00004
00005 using namespace std;
00006
00007 vu11121::vu11121() {}
00008
00009 vu11121::vu11121(byte *data, dword XSize, dword YSize, dword ZSize)
00010 {
00011 m_Data = data;
00012 m_DataSize = XSize * YSize * ZSize;
00013 m_Dim1Size = XSize;
00014 m_Dim2Size = YSize;
00015 m_Dim3Size = ZSize;
00016 }
00017
00018
00019
00020
00021
00022 bool vu11121::read(void)
00023 {
00024 if (m_FileName.isEmpty()) return setError("No file name specified.");
00025
00026 FILE *file = fopen(m_FileName,"rb");
00027 if (file != NULL)
00028 {
00029 bool success = read(file);
00030 fclose(file);
00031 return success;
00032 }
00033 else
00034 return setError("Could not open the specified file.");
00035 }
00036
00037
00038
00039
00040
00041 bool vu11121::write(void)
00042 {
00043 if (m_FileName.isEmpty()) return setError("No file name specified.");
00044
00045 FILE *file = fopen(m_FileName,"wb");
00046 if (file != NULL)
00047 {
00048 bool success = write(file);
00049 fclose(file);
00050 return success;
00051 }
00052 else
00053 return setError("Could not open the specified file.");
00054 }
00055
00056
00057
00058
00059
00060 bool vu11121::read(FILE *file)
00061 {
00062 int ret = 0;
00063 dword size = 0;
00064 int len = 0;
00065
00066
00067 bool success = vu1112::read(file);
00068 if (!success) return false;
00069
00070
00071 char dataName[64];
00072 ret = fscanf(file,"SCALARS %s ",dataName);
00073 if (ret != 1) return setInvalidFormatError();
00074
00075 m_DataName = dataName;
00076
00077
00078 fscanf(file,"byte LOOKUP_TABLE default%n",&len);
00079 if ((len < 25) || (fgetc(file) != '\n')) return setInvalidFormatError();
00080
00081
00082 if (m_DataSize != (m_Dim1Size * m_Dim2Size * m_Dim3Size))
00083 return setInvalidFormatError();
00084
00085
00086 m_Data = new byte[m_DataSize];
00087
00088
00089 if (m_Binary)
00090 size = fread(m_Data,sizeof(byte),m_DataSize,file);
00091 else
00092 {
00093 ret = 1;
00094 dword i = 0;
00095 while ((ret > 0) && (i < m_DataSize))
00096 ret = fscanf(file," %c",&m_Data[i++]);
00097 size = i;
00098 }
00099
00100
00101 if (size != m_DataSize) return setInvalidFormatError();
00102
00103 return true;
00104 }
00105
00106
00107
00108
00109
00110 bool vu11121::write(FILE *file)
00111 {
00112 int ret;
00113 dword size = 0;
00114
00115 m_Binary = true;
00116
00117 bool success = vu1112::write(file);
00118 if (!success) return false;
00119
00120 fprintf(file,"SCALARS ");
00121
00122
00123 if (!m_DataName.isEmpty())
00124 fprintf(file,"%s ",m_DataName.c_str());
00125 else
00126 fprintf(file,"data ");
00127
00128 fprintf(file,"byte\nLOOKUP_TABLE default\n");
00129
00130
00131 if (m_Binary)
00132 size = fwrite(m_Data,sizeof(byte),m_DataSize,file);
00133 else
00134 {
00135 ret = 1;
00136 dword i = 0;
00137 while ((ret > 0) && (i < m_DataSize))
00138 ret = fprintf(file," %d",m_Data[i++]);
00139 size = i;
00140 }
00141
00142 if (size == m_DataSize)
00143 return true;
00144 else
00145 return setWriteError();
00146 }
00147
00148 void vu11121::generateLapWeightHistogram()
00149 {
00156 unsigned long laplace_hist[256];
00157 unsigned long surface_area[256];
00158
00159
00160
00161
00162 for (int i=0; i<256; i++)
00163 {
00164 laplace_hist[i] = 0;
00165 surface_area[i] = 0;
00166 }
00167
00168 for (unsigned int k=0; k<m_Dim3Size; k++)
00169 for (unsigned int j=0; j<m_Dim2Size; j++)
00170 for (unsigned int i=0; i<m_Dim1Size; i++)
00171 {
00172 unsigned int val = getDataValue(i,j,k);
00173
00174 laplace_hist[val] +=
00175 ((getDataValue(i, j , k) - getDataValue(i+2, j, k)) - (getDataValue(i-2, j , k) - getDataValue(i, j, k))) +
00176 ((getDataValue(i, j , k) - getDataValue(i, j+2, k)) - (getDataValue(i, j-2 , k) - getDataValue(i, j, k))) +
00177 ((getDataValue(i, j , k) - getDataValue(i, j, k+2)) - (getDataValue(i, j , k-2) - getDataValue(i, j, k)));
00178
00179 if (val > getDataValue(i+1,j,k))
00180 surface_area[val] += 1;
00181 if (val < getDataValue(i+1,j,k))
00182 surface_area[val] -= 1;
00183 if (val > getDataValue(i-1,j,k))
00184 surface_area[val] += 1;
00185 if (val < getDataValue(i-1,j,k))
00186 surface_area[val] -= 1;
00187
00188 if (val > getDataValue(i,j+1,k))
00189 surface_area[val] += 1;
00190 if (val < getDataValue(i,j+1,k))
00191 surface_area[val] -= 1;
00192 if (val > getDataValue(i,j-1,k))
00193 surface_area[val] += 1;
00194 if (val < getDataValue(i,j-1,k))
00195 surface_area[val] -= 1;
00196
00197 if (val > getDataValue(i,j,k+1))
00198 surface_area[val] += 1;
00199 if (val < getDataValue(i,j,k+1))
00200 surface_area[val] -= 1;
00201 if (val > getDataValue(i,j,k-1))
00202 surface_area[val] += 1;
00203 if (val < getDataValue(i,j,k-1))
00204 surface_area[val] -= 1;
00205 }
00206
00207 for (int i = 254; i>= 0; i--)
00208 {
00209 laplace_hist[i] += laplace_hist[i+1];
00210 surface_area[i] += surface_area[i+1];
00211 }
00212
00213 cout << "Writing laplacian weighted histogram to file..." << endl;
00214
00215 FILE *fp1, *fp2;
00216
00217 fp1=fopen("lp_hist.dat","w");
00218 fp2=fopen("lp_area_hist.dat","w");
00219
00220
00221 for (int i=0; i<256; i++)
00222 {
00223 fprintf(fp1, "%d %d\n", i, (int)laplace_hist[i]);
00224 if (surface_area[i] != 0)
00225 fprintf(fp2, "%d %f\n", i, double(laplace_hist[i])/double(surface_area[i]));
00226 }
00227
00228 fclose(fp1);
00229 fclose(fp2);
00230
00231 }
00232
00233 void vu11121::cropFrom(const vu11121& vol, word cube[6])
00234 {
00235 if(cube[1]>vol.m_Dim1Size) cube[1] = vol.m_Dim1Size;
00236 if(cube[3]>vol.m_Dim2Size) cube[3] = vol.m_Dim2Size;
00237 if(cube[5]>vol.m_Dim3Size) cube[5] = vol.m_Dim3Size;
00238 if(cube[0]<cube[1] &&
00239 cube[2]<cube[3] &&
00240 cube[4]<cube[5] )
00241 {
00242 operator=(vol);
00243 if(m_Data) {
00244 delete [] m_Data;
00245 m_Data = NULL;
00246 }
00247 m_Dim1Size = cube[1] - cube[0];
00248 m_Dim2Size = cube[3] - cube[2];
00249 m_Dim3Size = cube[5] - cube[4];
00250 m_DataSize = m_Dim1Size*m_Dim2Size*m_Dim3Size;
00251 m_Data = new byte[m_DataSize];
00252 dword i,j,k;
00253 dword srclayer = vol.m_Dim1Size*vol.m_Dim2Size;
00254 dword layer = m_Dim1Size*m_Dim2Size;
00255 dword offset = srclayer*cube[4] + vol.m_Dim1Size*cube[2] + cube[0];
00256 for(k=0; k<m_Dim3Size; k++)
00257 for(j=0; j<m_Dim2Size; j++)
00258 {
00259 const byte *src = &vol.m_Data[
00260 srclayer*k + vol.m_Dim1Size*j + offset];
00261 byte *dst = &m_Data[layer*k + m_Dim1Size*j];
00262 for(i=0; i<m_Dim1Size; i++, src++, dst++)
00263 {
00264 *dst = *src;
00265 }
00266 }
00267 }
00268 }
00269
00270 void vu11121::scaleFrom(const vu11121& vol, word nx, word ny, word nz)
00271 {
00272 operator=(vol);
00273
00274 m_Dim1Size = nx;
00275 m_Dim2Size = ny;
00276 m_Dim3Size = nz;
00277 if(m_Data) {
00278 delete [] m_Data;
00279 m_Data = NULL;
00280 }
00281 m_DataSize = m_Dim1Size*m_Dim2Size*m_Dim3Size;
00282 m_Data = new byte[m_DataSize];
00283 dword i,j,k;
00284 dword layer = m_Dim1Size*m_Dim2Size;
00285 for(k=0; k<m_Dim3Size; k++)
00286 for(j=0; j<m_Dim2Size; j++)
00287 {
00288 byte *dst = &m_Data[layer*k + m_Dim1Size*j];
00289 for(i=0; i<m_Dim1Size; i++, dst++)
00290 {
00291
00292
00293
00294
00295
00296
00297
00298
00299 *dst = vol.getDataValue(i*vol.m_Dim1Size/m_Dim1Size,
00300 j*vol.m_Dim2Size/m_Dim2Size,
00301 k*vol.m_Dim3Size/m_Dim3Size);
00302 }
00303 }
00304 }
00305
00306 bool readRAW(const vuString& fname, vu11121& vol)
00307 {
00308 ifstream vf(fname);
00309 if(!vf.good()) return false;
00310 char line[2048];
00311 vf.getline(line, 2047, 0x20);
00312 if(sscanf(line, "%lu", &vol.m_Dim1Size) != 1)
00313 return false;
00314 vf.getline(line, 2047, 0x20);
00315 if(sscanf(line, "%lu", &vol.m_Dim2Size) != 1)
00316 return false;
00317 vf.getline(line, 2047, 0x20);
00318 if(sscanf(line, "%lu", &vol.m_Dim3Size) != 1)
00319 return false;
00320
00321
00322
00323
00324
00325
00326 if(vol.m_Data) delete [] vol.m_Data;
00327 vol.m_DataSize = vol.m_Dim1Size*vol.m_Dim2Size*vol.m_Dim3Size;
00328 vol.m_Data = new byte[vol.m_DataSize];
00329
00330 vf.read((char *)vol.m_Data, (int)vol.m_DataSize);
00331 if(vf.gcount() != (long)vol.m_DataSize || vf.bad())
00332 {
00333 cout << "gcount=" << vf.gcount() << " m_DataSize=" << vol.m_DataSize
00334 << " bad=" << (vf.bad()?"on":"off")
00335 << "eof" << (vf.eof() ? "on" : "off") << endl;
00336 return false;
00337 }
00338
00339 return true;
00340 }
00341
00342 bool vu11121::createHistogram(vuHistogram& hist) const
00343 {
00344 if(hist.getType() == vuHistogram::TYPE_INTENSITY) {
00345 hist.reset();
00346 const byte *dat = m_Data;
00347 dword p;
00348 for(p=0; p<m_Dim1Size*m_Dim2Size*m_Dim3Size; p++, dat++)
00349 hist.recordIntensity(*dat);
00350 return true;
00351 } else return false;
00352 }
00353
00354 void vu11121::remap(const vuMap& map)
00355 {
00356 dword i;
00357 byte *dat = m_Data;
00358 for(i=m_Dim1Size*m_Dim2Size*m_Dim3Size; i>0; i--, dat++)
00359 *dat = map[*dat];
00360 }