RE: Tiff issue
Valeri Onuchin <Valeri.Onoutchine <at> cern.ch>
2012-07-03 11:29:24 GMT
Ooops failed to send because of attachment size limit.
________________________________________
From: Valeri Onuchin
Sent: 03 July 2012 13:24
To: Esch, Shorty; roottalk (Mailing list for ROOT users.)
Subject: RE: Tiff issue
Hi Shorty,
the tiff format is not supported yet.
If you are working under windows
the attached codes demonstrate
- how to read a 16bit TIFF fie, convert it to "raw" format
- how to read the produced "raw" file in ROOT and
fill 2D hist with scaled red component.
HTH.
Regards. Valeriy
________________________________________
From: Esch, Shorty [ernst <at> lanl.gov]
Sent: 02 July 2012 19:02
To: roottalk (Mailing list for ROOT users.)
Subject: Tiff issue
Hi, has anybody a way to read a TIFF file grayscale into a 2d histogram in root?
Thanks,
Shorty
//$Id: doseHist.C 1092 2012-04-09 07:46:57Z onuchin $
// Author: Valeriy Onuchin 28.07.2011
/*************************************************************************
* *
* Copyright (C) 2012, Valeriy Onuchin *
* All rights reserved. *
* *
*************************************************************************/
#include <stdio.h>
#include "TSystem.h"
#include "TH2.h"
#include "TString.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TSpectrum2.h"
#include "TAxis.h"
#include "TFile.h"
//#include "LandscapeRed.h"
#include "PortraitRed.h"
#include "icru49.h"
double *dose_vs_LR;
//double *dose_vs_PR;
TH2D *dose2d = 0;
TH1D *dose1d = 0;
double *dose;
TCanvas *colz;
TCanvas *contz;
TCanvas *surf;
TCanvas *dose1dCanvas;
TCanvas *peakCanvas;
TH2 *dose2rebin = 0;
TH2 *dose2surf = 0;
TSpectrum2 *spectrum;
Float_t *xpeaks;
Float_t *ypeaks;
TFile *rootFile;
//______________________________________________________________________________
void doseHist(const char *file, const char *calibr = "PortraitRed.h")
{
// read RAW data file and create hists from this
int min, arrsz;
Bool_t land = false;
int ret = 0;//gROOT->LoadMacro(calibr);
if (ret) {
fprintf(stderr, "Failed to open calibr file: %s\n", calibr);
}
if (land) {
min = int(dose_vs_LR[0]);
arrsz = sizeof(dose_vs_LR)/sizeof(double) - 1;
} else {
min = int(dose_vs_PR[0]);
arrsz = sizeof(dose_vs_PR)/sizeof(double) - 1;
}
TString fname = gSystem->UnixPathName(file);
Bool_t tifile = kFALSE;
TString tmp_file = "";
TString rootfileName = gSystem->BaseName(fname);
// three possible extensions replace with .root
rootfileName.ReplaceAll(".tiff", ".root");
rootfileName.ReplaceAll(".tif", ".root");
rootfileName.ReplaceAll(".raw", ".root");
if (fname.EndsWith(".tif", TString::kIgnoreCase) ||
fname.EndsWith(".tiff", TString::kIgnoreCase)) {
tifile = kTRUE;
char *tmpname = gSystem->ConcatFileName(gSystem->TempDirectory(), gSystem->BaseName(fname));
TString tmp_file = tmpname;
tmp_file.ReplaceAll(".tif", ".raw");
delete tmpname;
TString tif2raw =gSystem->ExpandPathName("$(ROOTSYS)/macros/");
tif2raw += "tif2raw.exe -f " + fname + " -o " + tmp_file;
//printf("Convert tif file to raw: %s \n", tif2raw.Data());
gSystem->Exec(tif2raw);
fname = tmp_file;
}
//printf("%d\n", min);
FILE *fp = fopen(fname.Data(), "rb");
if (!fp) {
fprintf(stderr, "Failed to open file: %s", fname.Data());
return;
}
UInt_t w, h;
fread(&w, sizeof(w), 1, fp);
fread(&h, sizeof(h), 1, fp);
//printf("%d %d\n", w, h);
long sz = w*h*3;
unsigned short *rgb = new unsigned short[sz];
unsigned short r, g, b;
fread(rgb, sizeof(short), sz, fp);
fclose(fp);
// delete temporary file
if (tifile) {
gSystem->Unlink(tmp_file.Data());
}
fname = file;
TString title = "Dose distribution of ";
title += file;
TString name = "colz_";
name += fname;
// open ROOT file to save all hists
rootFile = new TFile(rootfileName, "RECREATE");
float dpi = 25.4/300;
dose2d = new TH2D(name, title, w, 0, w*dpi, h, 0, h*dpi);
dose2d->SetStats(0); //delete stats box
name = "1D_";
name += fname;
dose1d = new TH1D(name, title, 10000, 0, 50);
//printf("Start filling hitogram %s\n", dose2d->GetName());
Double_t value;
long j = 0;
long idx;
for (int y = 0; y < h; y++) {
for (int x = 0; x < w; x++) {
r = rgb[j];
//
idx = (r - min) > arrsz ? arrsz : r - min;
if (idx < 1) {
printf("x=%d y=%d r=%d\n", x, y, r);
}
idx = idx < 1 ? 1 : idx;
value = land ? dose_vs_LR[idx] : dose_vs_PR[idx];
g = rgb[++j];
b = rgb[++j];
dose2d->SetBinContent(x, h - y - 1, value);
dose1d->Fill(value);
j++;
}
}
delete [] rgb;
gStyle->SetCanvasPreferGL();
gStyle->SetOptStat(1001211);
title = "2D dose ";
title += fname;
colz = new TCanvas(fname + " colz", title + " colz in Gy", w*2-28, h*2-18);
colz->SetGridx();
colz->SetGridy();
colz->SetTickx();
colz->SetTicky();
//colz->SetCrosshair();
colz->ToggleEventStatus();
dose2d->SetTitle("Dose distribution of " + fname + " in Gy");
dose2d->Draw("colz");
dose2d->GetXaxis()->SetTitle("in mm");
dose2d->GetYaxis()->SetTitle("in mm");
colz->SetEditable(kFALSE);
int rebin = (int)300./25.4;
dose2rebin = dose2d->Rebin2D(10, 10, "rebinXY");
dose2rebin->SetStats(0); //delete stats box
dose2rebin->SetName("iso_" + fname);
dose2rebin->SetTitle("Isodose field of " + fname + " in cGy");
contz = new TCanvas(title + " contz", title + " contz", w*2-28, h*2-18);
contz->SetGridx();
contz->SetGridy();
contz->SetTickx();
contz->SetTicky();
//contz->SetCrosshair();
contz->ToggleEventStatus();
dose2rebin->Draw("cont1z");
dose2rebin->GetXaxis()->SetTitle("in mm");
dose2rebin->GetYaxis()->SetTitle("in mm");
//contz->SetEditable(kFALSE);
dose2surf = (TH2*)dose2rebin->Clone("dose2surf");
surf = new TCanvas(title + " surf", title + " surf in cGy");
dose2surf->Draw("glsurf2");
dose2surf->GetXaxis()->SetTitle("");
dose2surf->GetYaxis()->SetTitle("");
dose2surf->GetZaxis()->SetTitle("");
surf->SetGridx();
surf->SetGridy();
surf->SetTickx();
surf->SetTicky();
surf->ToggleEventStatus();
dose1dCanvas = new TCanvas(title + " 1D", "1D histogram of dose in Gy");
dose1dCanvas->SetGridx();
dose1dCanvas->SetGridy();
dose1dCanvas->SetTickx();
dose1dCanvas->SetTicky();
dose1dCanvas->SetLogy();
//dose1dCanvas->SetCrosshair();
dose1dCanvas->ToggleEventStatus();
dose1d->SetFillColor(2);
//dose1d->SetName("Dose_distribution_of_" + fname);
dose1d->SetTitle("Dose distribution of " + fname + " in Gy");
dose1d->Draw();
dose1d->GetXaxis()->SetTitle("in Gy");
//dose1dCanvas->SetEditable(kFALSE);
gStyle->SetCanvasPreferGL(0);
spectrum = new TSpectrum2(10);
TCanvas *tmp = new TCanvas();
Int_t nfound = spectrum->Search(dose2rebin, 4, "cont1z noMarkov");
delete tmp;
xpeaks = spectrum->GetPositionX();
ypeaks = spectrum->GetPositionY();
if (nfound >= 1) {
//peakCanvas = new TCanvas(name + " slices", name + " silces", 1024, 800);
//peakCanvas->SetCrosshair();
//peakCanvas->ToggleEventStatus();
//peakCanvas->Divide(2, nfound);
//peakCanvas->cd();
TH1D *xproj;
TH1D *yproj;
float idp = 300./25.4;
for (int i = 0; i < nfound; i++) {
//peakCanvas->cd(i+1);
printf("x = %f y = %f\n", xpeaks[i], ypeaks[i]);
TString prname = TString::Format("Y[%d] slice at X = ", i);
peakCanvas = new TCanvas(prname);
gPad->SetGridx();
gPad->SetGridy();
gPad->SetTickx();
gPad->SetTicky();
peakCanvas->ToggleEventStatus();
int xbin = int(xpeaks[i]*idp);
printf("i = %d xbin = %d\n" , i, xbin);
yproj = dose2d->ProjectionY(TString::Format("X = %f mm", xpeaks[i]), xbin, xbin, "d");
yproj->SetFillColor(9);
yproj->SetTitle(prname + TString::Format(" %f mm", xpeaks[i]));
//gPad->SetTitle(prname + TString::Format(" %f mm", xpeaks[i]));
//yproj->Draw();
yproj->GetXaxis()->SetTitle("in mm");
yproj->GetYaxis()->SetTitle("in Gy");
//peakCanvas->SetEditable(kFALSE);
//peakCanvas->cd(i+2);
prname = TString::Format("X[%d] slice at Y = ", i);
peakCanvas = new TCanvas(prname);
gPad->SetGridx();
gPad->SetGridy();
gPad->SetTickx();
gPad->SetTicky();
peakCanvas->ToggleEventStatus();
int ybin = int(ypeaks[i]*idp);
printf("i = %d ybin = %d\n" , i, ybin);
xproj = dose2d->ProjectionX(TString::Format("Y = %f mm", ypeaks[i]), ybin, ybin, "d");
xproj->SetFillColor(4);
xproj->SetTitle(prname + TString::Format(" %f mm", ypeaks[i]));
//gPad->SetTitle(prname + TString::Format(" %f mm", ypeaks[i]));
//xproj->Draw();
xproj->GetXaxis()->SetTitle("in mm");
xproj->GetYaxis()->SetTitle("in Gy");
//peakCanvas->SetEditable(kFALSE);
}
}
peakCanvas->Update();
rootFile->Write();
printf("Stop\n");
printf("\n\n==========================================================================\n");
double maxE = dose2d->GetMaximum();
double halfE = maxE*0.5;
TAxis *ax = dose1d->GetXaxis();
int b2 = ax->FindBin(maxE);
int b1 = ax->FindBin(halfE);
printf("Full integral = %E. Half integral = %E\n", dose2d->Integral(), dose1d->Integral(b1, b2));
int energy = 100; //100 MeV
//float dE1 = 7.2; // 7.2 MeV/cm at 100 MeV
//float dE2 = 9.57; // 9.57 MeV/cm
double iddp = (300./2.54)*(300./2.54); // (pixels per cantimeter)**2
double intensity = dose2d->Integral()/iddp/dEdx_ICRU49[energy]/10.*6.24e9; // dEdx_ICRU49 array -
stopping power of protons in water in Mev/mm
printf("Total Intensity = %E at proton Energy = %d (MeV)\n", intensity, energy);
}
//______________________________________________________________________________
void Intensity(int energy = 100)
{
// calculates total proton beam intensity based on total dose field integral
printf("\n\n==========================================================================\n");
double maxE = dose2d->GetMaximum();
double halfE = maxE*0.5;
TAxis *ax = dose1d->GetXaxis();
int b2 = ax->FindBin(maxE);
int b1 = ax->FindBin(halfE);
printf("Full integral = %E. Half integral = %E\n", dose2d->Integral(), dose1d->Integral(b1, b2));
double iddp = (300./2.54)*(300./2.54); // (pixels per cantimeter)**2
double intensity = dose2d->Integral()/iddp/dEdx_ICRU49[energy]/10.*6.24e9; // dEdx_ICRU49 array -
stopping power of protons in water in Mev/mm
printf("Total Intensity = %E at proton Energy = %d (MeV)\n", intensity, energy);
}
// $Id: tif2raw.cxx 945 2012-02-29 14:24:30Z onuchin $
// Author: Valeriy Onuchin 28.07.2011
/*************************************************************************
* *
* Copyright (C) 2011, Valeriy Onuchin *
* All rights reserved. *
* *
*************************************************************************/
#include <windows.h>
#include <cctype>
#include <stdio.h>
#include "vcclr.h" // for gcroot
using namespace System;
using namespace System::IO;
using namespace System::Drawing;
using namespace System::Windows::Media::Imaging;
#pragma warning(disable:4996) // CRT unsafe
#define LOCAL_DEBUG 1 // debugging
static char gUsageHelp[] = "\
Usage: %s -f infile\n\
Index:\n\
\"infile\" is a TIFF image file (must have \".tif\" or \".tiff\" extension)\n\
the output file will have the same name\n\
as \"infile\" but with \".raw\" extension.\n\n\
Example: tif2raw.exe -f scan111.tif\n\n";
static bool gTickIsOn = false;
//______________________________________________________________________________
static void tickThreadStub(void *)
{
//
while (1) {
::Sleep(1000);
#ifdef LOCAL_DEBUG
if (gTickIsOn) {
fprintf(stderr, ".");
}
#endif
}
}
////////////////////////////////////////////////////////////////////////////////
//______________________________________________________________________________
static char *in2out(char *infile)
{
// change file extension
char *out = 0;
out = strstr(infile, ".tif");
if (out) {
out[1] = 'r';
out[2] = 'a';
out[3] = 'w';
} else {
out = strstr(infile, ".tiff");
if (out) {
out[1] = 'r';
out[2] = 'a';
out[3] = 'w';
out[4] = '\x0';
} else {
fprintf(stderr, "Input file must have extension \".dat\"");
return 0;
}
}
return infile;
}
//______________________________________________________________________________
static char *basename(char *pathname)
{
static char *result;
char *filename = NULL;
DWORD len;
free(result);
len = 0;
len = ::GetFullPathNameA(pathname, len, result, &filename);
result = (char*)malloc(len);
if (result) {
if (::GetFullPathNameA(pathname, len, result, &filename) <= len) {
return filename;
}
}
return filename;
}
////////////////////////////////////////////////////////////////////////////////
//______________________________________________________________________________
int main(int argc, char *argv[])
{
//
char *optarg;
char rawfile[512];
char tiffile[512];
memset(rawfile, 0, sizeof(rawfile));
memset(tiffile, 0, sizeof(tiffile));
bool hasinput = false;
if ((argc > 1) && (!_stricmp(argv[1], "-h") || !_stricmp(argv[1], "/h") ||
!_stricmp(argv[1], "--help") ||
!_stricmp(argv[1], "-?") || !_stricmp(argv[1], "/?")) || (argc == 1)) {
fprintf(stderr, gUsageHelp, basename(argv[0]));
exit(0);
}
// check flag
if (!_stricmp(argv[1], "-f")) {
optarg = argv[2];
while (isspace(*optarg)) {
optarg++;
}
strncpy(tiffile, optarg, sizeof(tiffile) - 1);
tiffile[sizeof(tiffile) - 1] = '\x0';
hasinput = true;
}
if (!hasinput) {
fprintf(stderr, gUsageHelp, basename(argv[0]));
exit(0);
}
HANDLE tickThread = ::CreateThread(NULL, NULL,
(LPTHREAD_START_ROUTINE)tickThreadStub,
0, NULL, NULL);
#ifdef LOCAL_DEBUG
fprintf(stderr, "Reading TIFF file: %s\n", tiffile);
#endif
gTickIsOn = true;
FileStream^ stream;
try {
stream = gcnew FileStream(gcnew String(tiffile),
FileMode::Open,
FileAccess::Read,
FileShare::Read);
} catch (Exception^ e) {
Console::WriteLine("Failed to read file. Reason: {0}", e->Message);
throw;
}
Windows::Media::Imaging::TiffBitmapDecoder^ decoder;
try {
decoder = gcnew Windows::Media::Imaging::TiffBitmapDecoder(stream,
BitmapCreateOptions::PreservePixelFormat,
BitmapCacheOption::Default);
} catch (Exception^ e) {
Console::WriteLine("Failed to decode TIFF format. Reason: {0}", e->Message);
throw;
}
Windows::Media::Imaging::BitmapFrame^ bitmap = decoder->Frames[0];
int stride = bitmap->PixelWidth*6;
int picSize = bitmap->PixelWidth*bitmap->PixelHeight*3; // w*h*(rgb == 3 bytes)
array<short>^ pixels = gcnew array<short>(picSize);
if (bitmap->Format == Windows::Media::PixelFormats::Rgb48) {
bitmap->CopyPixels(pixels, stride, 0);
} else {
Console::WriteLine("Wrong image format: {0}", bitmap->Format);
}
stream->Close();
gTickIsOn = false;
if (tiffile[0]) {
char *out = in2out(tiffile);
if (!out) {
#ifdef LOCAL_DEBUG
fprintf(stderr, "Wrong output or input file name ... exiting\n");
#endif
exit(-1);
}
memcpy(rawfile, out, sizeof(rawfile));
}
#ifdef LOCAL_DEBUG
fprintf(stderr, "Writing RAW data to file: %s, size = %d bytes\n\t\tWidth = %d. Height = %d\n",
rawfile, (picSize + 2)*2, bitmap->PixelWidth, bitmap->PixelHeight);
#endif
FileStream^ fs = gcnew FileStream(gcnew String(rawfile), FileMode::Create);
BinaryWriter^ writer = gcnew BinaryWriter(fs);
writer->Write(bitmap->PixelWidth);
writer->Write(bitmap->PixelHeight);
for (int i = 0; i < picSize; i++) {
writer->Write(pixels[i]);
}
fs->Close();
#ifdef LOCAL_DEBUG
fprintf(stderr, "Done\n\n");
#endif
exit(0);
}
RSS Feed