From 86f7a69b37e26dc298cbe25ec2523073cbbecab6 Mon Sep 17 00:00:00 2001 From: Peter Stockings Date: Sat, 17 May 2025 09:21:16 +1000 Subject: [PATCH] Initial commit --- .static | 0 CTAccordion .js | 108 +++++++++++++ UploadImageComponent.js | 326 ++++++++++++++++++++++++++++++++++++++++ app.js | 3 + fbp.js | 92 ++++++++++++ index.html | 27 ++++ sinogram.js | 209 ++++++++++++++++++++++++++ 7 files changed, 765 insertions(+) create mode 100644 .static create mode 100644 CTAccordion .js create mode 100644 UploadImageComponent.js create mode 100644 app.js create mode 100644 fbp.js create mode 100644 index.html create mode 100644 sinogram.js diff --git a/.static b/.static new file mode 100644 index 0000000..e69de29 diff --git a/CTAccordion .js b/CTAccordion .js new file mode 100644 index 0000000..d0b173c --- /dev/null +++ b/CTAccordion .js @@ -0,0 +1,108 @@ +export const steps = [ + { + title: "1. Upload / Input Image", + content: `You upload or use a default grayscale image. This is treated like a 2D slice of a physical object. + + The input image is a 2D function: + $ f(x, y) $ + This represents how much X-rays are absorbed at each point.`, + }, + { + title: "2. Radon Transform (Generating the Sinogram)", + content: `We rotate a virtual X-ray beam around the image and compute line integrals at each angle — simulating how X-rays pass through. + +\\[ +R(\\theta, s) = \\int_{-\\infty}^{\\infty} f(s \\cos\\theta - t \\sin\\theta,\\ s \\sin\\theta + t \\cos\\theta)\\ dt +\\] + +Where:
+- $\\theta$ = projection angle
+- $s$ = offset from center + + ![Radon transform](https://upload.wikimedia.org/wikipedia/commons/9/93/Radon_transform_sinogram.gif)`, + }, + { + title: "3. Sinogram", + content: `You now have a 2D image where: + - X-axis = detector position + - Y-axis = angle + Each row is a projection. + + Math: + $\text{Sinogram} = \{ R(\theta_1, s), R(\theta_2, s), \dots, R(\theta_n, s) \}$ + + ![Sinogram](https://www.researchgate.net/profile/Samuel-Asante-2/publication/299856137/figure/fig3/AS:348226420002817@1460035056245/A-Shepp-Logan-Phantom-and-reconstructed-Image-Sinogram-a-Original-image-b-radon.png)`, + }, + { + title: "4. Optional: Apply Ramp Filter (FBP)", + content: `If enabled, we sharpen each projection before back-projecting by amplifying high-frequency content. + + Math: +\\[ +P(\\omega) = \\mathcal{F}[R(\\theta, s)] \\\\ +P'(\\omega) = |\\omega| \\cdot P(\\omega) \\\\ +\\text{Filtered Projection} = \\mathcal{F}^{-1}[P'(\\omega)] +\\] + + ![Ramp filter graph](https://www.researchgate.net/publication/346858231/figure/fig4/AS:967035467091970@1607570630558/Applying-Ramp-filter-to-a-sinogram-preserves-high-frequency-features-The-filter-is.png)`, + }, + { + title: "5. Back Projection", + content: `Each (filtered) projection is \"smeared\" back into the image space along its angle. + + \\[ +f'(x, y) = \\int_0^{\\pi} R'(\\theta,\\ x \\cos\\theta + y \\sin\\theta)\\ d\\theta +\\]`, + }, + { + title: "6. Final Reconstruction", + content: `After all angles are added up, you get a reconstructed image resembling the original. + + \[ + f'(x, y) \approx f(x, y) + \]`, + }, +]; + +export const StepAccordion = { + view(vnode) { + const index = vnode.attrs.index; + const step = steps[index]; + const expanded = vnode.state.expanded ?? false; + + return m("div", { class: "border-b border-gray-300 py-1 my-2" }, [ + m( + "button", + { + class: + "w-full text-left font-bold text-lg text-gray-800 focus:outline-none", + onclick: () => { + vnode.state.expanded = !vnode.state.expanded; + }, + }, + step.title + ), + vnode.state.expanded && + m( + "div", + { + class: "mt-2 text-gray-700 whitespace-pre-wrap text-sm", + onupdate: () => { + if (window.MathJax) window.MathJax.typesetPromise(); + }, + }, + m.trust(markdownToHTML(step.content)) + ), + ]); + }, +}; + +// You must provide a markdownToHTML() function or use a library like marked.js +function markdownToHTML(text) { + return text + .replace(/\n/g, "
") + .replace( + /!\[(.*?)\]\((.*?)\)/g, + '$1' + ); +} diff --git a/UploadImageComponent.js b/UploadImageComponent.js new file mode 100644 index 0000000..90b6f5e --- /dev/null +++ b/UploadImageComponent.js @@ -0,0 +1,326 @@ +import { + generateSinogram, + reconstructImageFromSinogram, + convertToGrayscale, +} from "./sinogram.js"; + +import { StepAccordion } from "./CTAccordion .js"; + +export const UploadImageComponent = { + hasLoadedInitialImage: false, + angleCount: 180, + imageUrl: + "https://upload.wikimedia.org/wikipedia/commons/e/e5/Shepp_logan.png", + sinogramUrl: null, + reconstructedUrl: null, + defaultImageUrl: + "https://upload.wikimedia.org/wikipedia/commons/e/e5/Shepp_logan.png", + reconstructionFrames: [], + currentFrameIndex: 0, + renderMode: "grayscale", // or "heatmap" + useFBP: true, + + drawAngleOverlay(theta) { + const canvas = this.overlayCanvas; + if (!canvas || !this.imageElement) return; + + const ctx = canvas.getContext("2d"); + const w = canvas.width; + const h = canvas.height; + const cx = w / 2; + const cy = h / 2; + const len = Math.max(w, h); + + ctx.clearRect(0, 0, w, h); + + const dx = len * Math.cos(theta); + const dy = len * Math.sin(theta); + + ctx.strokeStyle = "rgba(255,0,0,0.8)"; + ctx.lineWidth = 2; + + ctx.beginPath(); + ctx.moveTo(cx - dx, cy - dy); + ctx.lineTo(cx + dx, cy + dy); + ctx.stroke(); + }, + + isOverlayReady() { + return ( + this.overlayCanvas && + this.imageElement && + this.imageElement.complete && + this.imageElement.naturalWidth > 0 + ); + }, + + async loadAndProcess(url, isUploaded = false) { + this.imageUrl = url; + this.sinogramUrl = "loading"; + this.reconstructedUrl = null; + m.redraw(); + + this.imageUrl = url; + let finalUrl = url; + + if (isUploaded) { + finalUrl = await convertToGrayscale(url); + this.imageUrl = finalUrl; + } + + this.sinogramUrl = await generateSinogram( + finalUrl, + this.angleCount, + this.drawAngleOverlay.bind(this) + ); + m.redraw(); + + this.reconstructionFrames = []; + this.currentFrameIndex = 0; + + this.reconstructedUrl = await reconstructImageFromSinogram( + this.sinogramUrl, + undefined, + (angle, frameUrl) => { + this.reconstructionFrames.push(frameUrl); + this.currentFrameIndex = this.reconstructionFrames.length - 1; + this.reconstructedUrl = frameUrl; + m.redraw(); + }, + this.renderMode, + this.useFBP + ); + }, + + oninit() { + this.loadAndProcessDebounced = debounce((url) => { + this.loadAndProcess(url); + }, 300); + }, + + view() { + return m( + "div", + { + class: "flex flex-col items-center min-h-screen bg-gray-100 py-10 px-4", + }, + [ + // Header + m("header", { class: "mb-10 text-center" }, [ + m( + "h1", + { class: "text-4xl font-bold text-gray-800 mb-2" }, + "Sinogram Generator" + ), + m( + "p", + { class: "text-gray-600 text-lg" }, + "Upload a grayscale image to simulate CT scan projections" + ), + ]), + + m(StepAccordion, { index: 0 }), + + // Upload Box + m( + "div", + { + class: + "w-full max-w-lg border-4 border-dashed border-gray-400 bg-white rounded-xl p-6 text-center hover:bg-gray-50 cursor-pointer transition", + ondragover: (e) => e.preventDefault(), + ondrop: (e) => { + e.preventDefault(); + const file = e.dataTransfer.files[0]; + if (file && file.type.startsWith("image/")) { + const url = URL.createObjectURL(file); + this.loadAndProcess(url, true); + } + }, + onclick: () => document.getElementById("fileInput").click(), + }, + [ + m( + "p", + { class: "text-gray-500" }, + "Click or drag a grayscale image here" + ), + m("input", { + id: "fileInput", + type: "file", + class: "hidden", + accept: "image/*", + onchange: (e) => { + const file = e.target.files[0]; + if (file && file.type.startsWith("image/")) { + const url = URL.createObjectURL(file); + this.loadAndProcess(url); + } + }, + }), + ] + ), + + m(StepAccordion, { index: 1 }), + + // Image Preview + m("div", { class: "relative mt-6 w-full max-w-md" }, [ + m("img", { + src: this.imageUrl, + class: "rounded shadow max-w-full h-auto mx-auto", + onload: (e) => { + this.imageElement = e.target; + + // Only start once both image and canvas are ready + if (this.isOverlayReady() && !this.hasLoadedInitialImage) { + this.hasLoadedInitialImage = true; + this.loadAndProcess(this.imageUrl); + } + }, + }), + m("canvas", { + width: this.imageElement?.width || 0, + height: this.imageElement?.height || 0, + style: "position:absolute; top:0; left:0; pointer-events:none;", + oncreate: ({ dom }) => { + this.overlayCanvas = dom; + + // Trigger load if image was already ready + if (this.isOverlayReady() && !this.hasLoadedInitialImage) { + this.hasLoadedInitialImage = true; + this.loadAndProcess(this.imageUrl); + } + }, + }), + ]), + + // Angle Slider + m("div", { class: "mt-6 w-full max-w-md" }, [ + m( + "label", + { class: "block text-sm font-medium text-gray-700 mb-1" }, + `Number of angles: ${this.angleCount}` + ), + m("input", { + type: "range", + min: 5, + max: 360, + value: this.angleCount, + step: 1, + class: "w-full", + oninput: (e) => { + this.angleCount = parseInt(e.target.value, 10); + this.loadAndProcessDebounced(this.imageUrl); // reprocess with new angle count + }, + }), + ]), + + m(StepAccordion, { index: 2 }), + + // Sinogram + this.sinogramUrl && + m("div", { class: "mt-10 w-full max-w-md text-center" }, [ + m( + "h2", + { class: "text-xl font-semibold text-gray-700 mb-4" }, + "Generated Sinogram" + ), + this.sinogramUrl === "loading" + ? m("p", "Processing...") + : m("img", { + src: this.sinogramUrl, + alt: "Sinogram", + class: "rounded shadow max-w-full h-auto mx-auto", + }), + ]), + + m(StepAccordion, { index: 3 }), + + // Reconstructed + this.reconstructedUrl && + m("div", { class: "mt-10 w-full max-w-md text-center" }, [ + m( + "h2", + { class: "text-xl font-semibold text-gray-700 mb-4" }, + "Reconstructed Image (Back Projection)" + ), + m("img", { + src: this.reconstructionFrames[this.currentFrameIndex], + alt: "Reconstructed", + class: "rounded shadow max-w-full h-auto mx-auto", + }), + m("div", { class: "mt-6 w-full max-w-md text-center" }, [ + m( + "label", + { class: "text-sm text-gray-600 mr-2" }, + "Render style:" + ), + m( + "select", + { + value: this.renderMode, + onchange: (e) => { + this.renderMode = e.target.value; + this.loadAndProcess(this.imageUrl); // re-render using selected mode + }, + }, + [ + m("option", { value: "heatmap" }, "Heatmap"), + m("option", { value: "grayscale" }, "Grayscale"), + ] + ), + ]), + + m("div", { class: "mt-4 w-full max-w-md text-left" }, [ + m("label", [ + m("input", { + type: "checkbox", + checked: this.useFBP, + onchange: (e) => { + this.useFBP = e.target.checked; + this.loadAndProcess(this.imageUrl); // regenerate with or without FBP + }, + }), + m( + "span", + { class: "ml-2 text-gray-700" }, + "Use Filtered Back Projection (Ramp)" + ), + ]), + ]), + + this.reconstructionFrames.length > 1 && + m("div", { class: "mt-4" }, [ + m("input", { + type: "range", + min: 0, + max: this.reconstructionFrames.length - 1, + value: this.currentFrameIndex, + step: 1, + oninput: (e) => { + this.currentFrameIndex = +e.target.value; + this.reconstructedUrl = + this.reconstructionFrames[this.currentFrameIndex]; + }, + }), + m( + "p", + { class: "text-sm text-gray-500 mt-1" }, + `Angle ${this.currentFrameIndex + 1} / ${ + this.reconstructionFrames.length + }` + ), + ]), + m(StepAccordion, { index: 4 }), + ]), + ] + ); + }, +}; + +function debounce(fn, delay) { + let timeout; + return (...args) => { + clearTimeout(timeout); + timeout = setTimeout(() => fn(...args), delay); + }; +} diff --git a/app.js b/app.js new file mode 100644 index 0000000..a2eab36 --- /dev/null +++ b/app.js @@ -0,0 +1,3 @@ +import { UploadImageComponent } from "./UploadImageComponent.js"; + +m.mount(document.body, UploadImageComponent); diff --git a/fbp.js b/fbp.js new file mode 100644 index 0000000..e9412f9 --- /dev/null +++ b/fbp.js @@ -0,0 +1,92 @@ +// Helper: Next power of 2 +function nextPow2(n) { + return 1 << (32 - Math.clz32(n - 1)); +} + +// Basic Cooley-Tukey FFT (real/imag arrays) +function fft1D(re, im) { + const N = re.length; + if (N <= 1) return; + + // Bit-reversal permutation + const rev = new Uint32Array(N); + let logN = Math.log2(N); + for (let i = 0; i < N; i++) { + let x = i, + y = 0; + for (let j = 0; j < logN; j++) { + y <<= 1; + y |= x & 1; + x >>= 1; + } + rev[i] = y; + } + for (let i = 0; i < N; i++) { + if (i < rev[i]) { + [re[i], re[rev[i]]] = [re[rev[i]], re[i]]; + [im[i], im[rev[i]]] = [im[rev[i]], im[i]]; + } + } + + // FFT + for (let s = 1; s <= logN; s++) { + const m = 1 << s; + const m2 = m >> 1; + const wAngle = (-2 * Math.PI) / m; + const cosW = Math.cos(wAngle); + const sinW = Math.sin(wAngle); + for (let k = 0; k < N; k += m) { + let wr = 1, + wi = 0; + for (let j = 0; j < m2; j++) { + const tRe = wr * re[k + j + m2] - wi * im[k + j + m2]; + const tIm = wr * im[k + j + m2] + wi * re[k + j + m2]; + const uRe = re[k + j]; + const uIm = im[k + j]; + re[k + j] = uRe + tRe; + im[k + j] = uIm + tIm; + re[k + j + m2] = uRe - tRe; + im[k + j + m2] = uIm - tIm; + const tempWr = wr; + wr = wr * cosW - wi * sinW; + wi = tempWr * sinW + wi * cosW; + } + } + } +} + +// Inverse FFT +function ifft1D(re, im) { + // Conjugate + for (let i = 0; i < re.length; i++) im[i] = -im[i]; + fft1D(re, im); + // Normalize and re-conjugate + const N = re.length; + for (let i = 0; i < N; i++) { + re[i] /= N; + im[i] = -im[i] / N; + } +} + +// Apply ramp filter in frequency domain +export function applyRampFilter(projection) { + const N = nextPow2(projection.length); + const re = new Float32Array(N); + const im = new Float32Array(N); + for (let i = 0; i < projection.length; i++) { + re[i] = projection[i]; + } + + fft1D(re, im); + + for (let i = 0; i < N / 2; i++) { + const freq = i / N; + re[i] *= freq; + im[i] *= freq; + re[N - i - 1] *= freq; + im[N - i - 1] *= freq; + } + + ifft1D(re, im); + return Array.from(re.slice(0, projection.length)); +} diff --git a/index.html b/index.html new file mode 100644 index 0000000..3748e2e --- /dev/null +++ b/index.html @@ -0,0 +1,27 @@ + + + + + + + Sinogram App + + + + + + + + + + + + \ No newline at end of file diff --git a/sinogram.js b/sinogram.js new file mode 100644 index 0000000..5b86b9f --- /dev/null +++ b/sinogram.js @@ -0,0 +1,209 @@ +import { applyRampFilter } from "./fbp.js"; + +export async function generateSinogram( + imageUrl, + angles = 180, + drawAngleCallback = null +) { + const image = await loadImage(imageUrl); + const size = Math.max(image.width, image.height); + const projections = []; + + const canvas = Object.assign(document.createElement("canvas"), { + width: size, + height: size, + }); + const ctx = canvas.getContext("2d"); + + for (let angle = 0; angle < angles; angle++) { + const theta = (angle * Math.PI) / angles; + + // 🔁 Call visual overlay for this angle + if (drawAngleCallback) drawAngleCallback(theta); + + // (Optional: add delay for animation) + await new Promise((r) => setTimeout(r, 0.01)); + + // Clear canvas + ctx.clearRect(0, 0, size, size); + + // Transform and draw rotated image + ctx.save(); + ctx.translate(size / 2, size / 2); + ctx.rotate(theta); + ctx.drawImage(image, -image.width / 2, -image.height / 2); + ctx.restore(); + + // Read pixel data + const { data } = ctx.getImageData(0, 0, size, size); + + // Sum brightness vertically (simulate X-ray projection) + const projection = []; + for (let x = 0; x < size; x++) { + let sum = 0; + for (let y = 0; y < size; y++) { + const i = (y * size + x) * 4; + const gray = data[i]; // red channel (since grayscale) + sum += gray; + } + projection.push(sum / size); // normalize + } + projections.push(projection); + } + + // Create sinogram canvas + const sinogramCanvas = Object.assign(document.createElement("canvas"), { + width: size, + height: angles, + }); + const sinCtx = sinogramCanvas.getContext("2d"); + const imgData = sinCtx.createImageData(size, angles); + + for (let y = 0; y < angles; y++) { + for (let x = 0; x < size; x++) { + const val = projections[y][x]; + const i = (y * size + x) * 4; + imgData.data[i + 0] = val; + imgData.data[i + 1] = val; + imgData.data[i + 2] = val; + imgData.data[i + 3] = 255; + } + } + + sinCtx.putImageData(imgData, 0, 0); + return sinogramCanvas.toDataURL(); +} + +export async function reconstructImageFromSinogram( + sinogramUrl, + size = 256, + onFrame = null, + renderMode = "heatmap", + useFBP = true +) { + const sinogramImage = await loadImage(sinogramUrl); + const canvas = Object.assign(document.createElement("canvas"), { + width: sinogramImage.width, + height: sinogramImage.height, + }); + const ctx = canvas.getContext("2d"); + ctx.drawImage(sinogramImage, 0, 0); + const sinogramData = ctx.getImageData( + 0, + 0, + sinogramImage.width, + sinogramImage.height + ).data; + + size = sinogramImage.width; // match size to sinogram resolution + const outputCanvas = Object.assign(document.createElement("canvas"), { + width: size, + height: size, + }); + const outputCtx = outputCanvas.getContext("2d"); + const accum = new Float32Array(size * size); + const center = size / 2; + + const angles = sinogramImage.height; + const width = sinogramImage.width; + + for (let angle = 0; angle < angles; angle++) { + const theta = (angle * Math.PI) / angles; + + let projection = []; + for (let x = 0; x < width; x++) { + const i = (angle * width + x) * 4; + projection.push(sinogramData[i]); + } + if (useFBP) { + projection = applyRampFilter(projection); + } + + for (let y = 0; y < size; y++) { + for (let x = 0; x < size; x++) { + const x0 = x - center; + const y0 = center - y; // flip y + const s = Math.round( + x0 * Math.cos(theta) + y0 * Math.sin(theta) + width / 2 + ); + if (s >= 0 && s < width) { + accum[y * size + x] += projection[s]; + } + } + } + + if (onFrame) { + // normalize and draw current frame + let maxVal = 0; + for (let i = 0; i < accum.length; i++) { + if (accum[i] > maxVal) maxVal = accum[i]; + } + const imageData = outputCtx.createImageData(size, size); + for (let i = 0; i < accum.length; i++) { + let val = accum[i] / maxVal; + val = Math.min(1, Math.max(0, val)); + let r, g, b; + if (renderMode === "grayscale") { + const gray = Math.round(val * 255); + r = g = b = gray; + } else { + [r, g, b] = getHeatmapColor(val); + } + imageData.data[i * 4 + 0] = r; + imageData.data[i * 4 + 1] = g; + imageData.data[i * 4 + 2] = b; + imageData.data[i * 4 + 3] = 255; + } + outputCtx.putImageData(imageData, 0, 0); + await new Promise((r) => setTimeout(r, 1)); + onFrame(angle, outputCanvas.toDataURL()); + } + } + + return outputCanvas.toDataURL(); +} + +// Heatmap mapping: blue → green → yellow → red +function getHeatmapColor(value) { + const r = Math.min(255, Math.max(0, 255 * Math.min(1, 4 * (value - 0.75)))); + const g = Math.min(255, Math.max(0, 255 * (4 * Math.abs(value - 0.5) - 1))); + const b = Math.min(255, Math.max(0, 255 * (1 - 4 * value))); + return [r, g, b]; +} + +function loadImage(src) { + return new Promise((resolve) => { + const img = new Image(); + img.crossOrigin = "anonymous"; + img.onload = () => resolve(img); + img.src = src; + }); +} + +export async function convertToGrayscale(imageUrl) { + const image = await loadImage(imageUrl); + const canvas = Object.assign(document.createElement("canvas"), { + width: image.width, + height: image.height, + }); + const ctx = canvas.getContext("2d"); + + // Draw original image + ctx.drawImage(image, 0, 0); + + // Get pixel data + const imageData = ctx.getImageData(0, 0, canvas.width, canvas.height); + const data = imageData.data; + + // Convert to grayscale: set R, G, B to luminance + for (let i = 0; i < data.length; i += 4) { + const r = data[i]; + const g = data[i + 1]; + const b = data[i + 2]; + const luminance = 0.299 * r + 0.587 * g + 0.114 * b; + data[i] = data[i + 1] = data[i + 2] = luminance; + } + + ctx.putImageData(imageData, 0, 0); + return canvas.toDataURL(); +}